solution.cc

Go to the documentation of this file.

 
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_generator.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
 
#include <fstream>
#include <iostream>
 
using namespace dealii;
 
 
template <int dim, int n_equations>
class Solution : public Function<dim>
{
public:
    Solution()
      : Function<dim>(n_equations)
    {}
 
    virtual void vector_value(const Point<dim>  &p,
                              Vector<double>    &values) const override;
 
    virtual void vector_gradient(const Point<dim>  &p,
                                 std::vector< Tensor< 1, dim, double > > & gradients ) const override;
 
private:
    int R = 1;
};
 
 
// Out-of-line implementation of Solution::vector_value() (see the
// class declaration for documentation). Several manufactured solutions
// are coded up here; all but one (test case 2) are commented out. To
// switch test cases, uncomment the desired block and comment out the
// others.
template <int dim, int n_equation>
void Solution<dim, n_equation>::vector_value(const Point<dim> &p,
                                             Vector<double> &values) const
{
    // ---- test case 1 (commented out) ---------------------------------
    //double pi = numbers::PI;
    //values(0) = std::sin(p[0]) * std::sin(p[1]);
    //values(1) = (- 2) * p[0] * p[1];
 
    // ---- test case 2 (currently active) ------------------------------
    values(0) = std::exp(p[0]) * std::cos(p[1]);
    values(1) = (- 2) * p[0] * p[1];
 
    // ---- test case 3 (commented out) ---------------------------------
    //double pi = numbers::PI;
    //values(0) = std::sin(pi*p[0]) * std::sin(pi*p[1]);
    //values(1) = (- 2) * p[0] * p[1];
}
 
// Out-of-line implementation of Solution::vector_gradient() (see the
// class declaration for documentation). Returns the analytic gradient
// of the currently active manufactured solution (test case 2). The
// function first resizes @p gradients to @p n_equations and then fills
// it component by component.
template <int dim, int n_equations>
void Solution<dim, n_equations>::vector_gradient(const Point<dim> &p,
                                                 std::vector< Tensor<1, dim, double > > & gradients ) const
{
    //double pi = numbers::PI;
    gradients.resize(n_equations); // Resize the gradients vector
    //std::vector< Tensor<1, dim, double> > return_value;
 
    // ---- test case 1 (commented out) ---------------------------------
    //double pi = numbers::PI;
    //gradients[0][0] = std::cos(p[0]) * std::sin(p[1]);
    //gradients[0][1] = std::sin(p[0]) * std::cos(p[1]);
    //gradients[1][0] = - 2 * p[1];
    //gradients[1][1] = - 2 * p[0];
 
    // ---- test case 2 (currently active) ------------------------------
    gradients[0][0] = std::exp(p[0]) * std::cos(p[1]);
    gradients[0][1] = std::exp(p[0]) * (-1) * std::sin(p[1]);
    gradients[1][0] = - 2 * p[1];
    gradients[1][1] = - 2 * p[0];
 
    // ---- test case 3 (commented out) ---------------------------------
    //gradients[0][0] =  pi*std::cos(pi*p[0]) * std::sin(pi*p[1]);
    //gradients[0][1] =  pi*std::sin(pi*p[0]) * std::cos(pi*p[1]);
    //gradients[1][0] = - 2 * p[1];
    //gradients[1][1] = - 2 * p[0];
 
    // ---- test case 4 (polar coordinates, commented out) --------------
    //const auto r = p.norm();
    //const auto theta = std::atan( p[1]/p[0] );
    //
    //gradients[0][0] = 0;
    //gradients[0][1] = 0;
    //gradients[1][0] = std::cos( theta / R);
    //gradients[1][1] = - r * (1/R) * std::sin( theta / R);
    //return return_value;
}