| Table of contents | |
|---|---|
The Heidelberg group of Professor Rolf Rannacher, to which the three main authors of the deal.II library belonged during their PhD time and partly also afterwards, has been involved with adaptivity and error estimation for finite element discretizations since the mid-90ies. The main achievement is the development of error estimates for arbitrary functionals of the solution, and of optimal mesh refinement for its computation.
We will not discuss the derivation of these concepts in too great detail, but will implement the main ideas in the present example program. For a thorough introduction into the general idea, we refer to the seminal work of Becker and Rannacher [BR95],[BR96r], and the overview article of the same authors in Acta Numerica [BR01]; the first introduces the concept of error estimation and adaptivity for general functional output for the Laplace equation, while the second gives many examples of applications of these concepts to a large number of other, more complicated equations. For applications to individual types of equations, see also the publications by Becker [Bec95], [Bec98], Kanschat [Kan96], [FK97], Suttmeier [Sut96], [RS97], [RS98c], [RS99], Bangerth [BR99b], [Ban00w], [BR01a], [Ban02], and Hartmann [Har02], [HH01], [HH01b]. All of these works, from the original introduction by Becker and Rannacher to individual contributions to particular equations, have later been summarized in a book by Bangerth and Rannacher that covers all of these topics, see [BR03].
The basic idea is the following: in applications, one is not usually interested in the solution per se, but rather in certain aspects of it. For example, in simulations of flow problems, one may want to know the lift or drag of a body immersed in the fluid; it is this quantity that we want to know to best accuracy, and whether the rest of the solution of the describing equations is well resolved is not of primary interest. Likewise, in elasticity one might want to know about values of the stress at certain points to guess whether maximal load values of joints are safe, for example. Or, in radiative transfer problems, mean flux intensities are of interest.
In all the cases just listed, it is the evaluation of a functional
of the solution which we are interested in, rather than the values of
everywhere. Since the exact solution
is not available, but only its numerical approximation
, it is sensible to ask whether the computed value
is within certain limits of the exact value
, i.e. we want to bound the error with respect to this functional,
.
For simplicity of exposition, we henceforth assume that both the quantity of interest
, as well as the equation are linear, and we will in particular show the derivation for the Laplace equation with homogeneous Dirichlet boundary conditions, although the concept is much more general. For this general case, we refer to the references listed above. The goal is to obtain bounds on the error,
. For this, let us denote by
the solution of a dual problem, defined as follows:
where
is the bilinear form associated with the differential equation, and the test functions are chosen from the corresponding solution space. Then, taking as special test function
the error, we have that
and we can, by Galerkin orthogonality, rewrite this as
for all possible functions
from the discrete test space.
Concretely, for Laplace's equation, the error identity reads
For reasons that we will not explain, we do not want to use this formula as is, but rather split the scalar products into terms on all cells, and integrate by parts on each of them:
Next we use that
, and that
is a quantity that is continuous almost everywhere, so the terms involving
on one cell cancels with that on its neighbor, where the normal vector has the opposite sign. At the boundary of the domain, where there is no neighbor cell with which this term could cancel, the weight
can be chosen as zero, since
has zero boundary values, and
can be chosen to have the same.
Thus, we have
In a final step, note that when taking the normal derivative of
, we mean the value of this quantity as taken from this side of the cell (for the usual Lagrange elements, derivatives are not continuous across edges). We then rewrite the above formula by exchanging half of the edge integral of cell
with the neighbor cell
, to obtain
Using that for the normal vectors
holds, we define the jump of the normal derivative by
and get the final form after setting the discrete function
, which is by now still arbitrary, to the point interpolation of the dual solution,
:
With this, we have obtained an exact representation of the error of the finite element discretization with respect to arbitrary (linear) functionals
. Its structure is a weighted form of a residual estimator, as both
and
are cell and edge residuals that vanish on the exact solution, and
are weights indicating how important the residuals on a certain cell is for the evaluation of the given functional. Furthermore, it is a cell-wise quantity, so we can use it as a mesh refinement criterion. The question, is: how to evaluate it? After all, the evaluation requires knowledge of the dual solution
, which carries the information about the quantity we want to know to best accuracy.
In some, very special cases, this dual solution is known. For example, if the functional
is the point evaluation,
, then the dual solution has to satisfy
with the Dirac delta function on the right hand side, and the dual solution is the Green's function with respect to the point
. For simple geometries, this function is analytically known, and we could insert it into the error representation formula.
However, we do not want to restrict ourselves to such special cases. Rather, we will compute the dual solution numerically, and approximate
by some numerically obtained
. We note that it is not sufficient to compute this approximation
using the same method as used for the primal solution
, since then
, and the overall error estimate would be zero. Rather, the approximation
has to be from a larger space than the primal finite element space. There are various ways to obtain such an approximation (see the cited literature), and we will choose to compute it with a higher order finite element space. While this is certainly not the most efficient way, it is simple since we already have all we need to do that in place, and it also allows for simple experimenting. For more efficient methods, again refer to the given literature, in particular [BR95], [BR03].
With this, we end the discussion of the mathematical side of this program and turn to the actual implementation.
The step-14 example program builds heavily on the techniques already used in the step-13 program. Its implementation of the dual weighted residual error estimator explained above is done by deriving a second class, properly called DualSolver, from the Solver base class, and having a class (WeightedResidual) that joins the two again and controls the solution of the primal and dual problem, and then uses both to compute the error indicator for mesh refinement.
The program continues the modular concept of the previous example, by implementing the dual functional, describing quantity of interest, by an abstract base class, and providing two different functionals which implement this interface. Adding a different quantity of interest is thus simple.
One of the more fundamental differences is the handling of data. A common case is that you develop a program that solves a certain equation, and test it with different right hand sides, different domains, different coefficients and boundary values, etc. Usually, these have to match, so that exact solutions are known, or that their combination makes sense at all.
We demonstrate a way how this can be achieved in a simple, yet very flexible way. We will put everything that belongs to a certain setup into one class, and provide a little C++ mortar around it, so that entire setups (domains, coefficients, right hand sides, etc.) can be exchanged by only changing something in one place.
Going this way a little further, we have also centralized all the other parameters that describe how the program is to work in one place, such as the order of the finite element, the maximal number of degrees of freedom, the evaluation objects that shall be executed on the computed solutions, and so on. This allows for simpler configuration of the program, and we will show in a later program how to use a library class that can handle setting these parameters by reading an input file. The general aim is to reduce the places within a program where one may have to look when wanting to change some parameter, as it has turned out in practice that one forgets where they are as programs grow. Furthermore, putting all options describing what the program does in a certain run into a file (that can be stored with the results) helps repeatability of results more than if the various flags were set somewhere in the program, where their exact values are forgotten after the next change to this place.
Unfortunately, the program has become rather long. While this admittedly reduces its usefulness as an example program, we think that it is a very good starting point for development of a program for other kinds of problems, involving different equations than the Laplace equation treated here. Furthermore, it shows everything that we can show you about our way of a posteriori error estimation, and its structure should make it simple for you to adjust this method to other problems, other functionals, other geometries, coefficients, etc.
The author believes that the present program is his masterpiece among the example programs, regarding the mathematical complexity, as well as the simplicity to add extensions. If you use this program as a basis for your own programs, we would kindly like to ask you to state this fact and the name of the author of the example program, Wolfgang Bangerth, in publications that arise from that, of your program consists in a considerable part of the example program.
Start out with well known things...
#include <base/quadrature_lib.h> #include <base/function.h> #include <base/logstream.h> #include <base/thread_management.h> #include <lac/vector.h> #include <lac/full_matrix.h> #include <lac/sparse_matrix.h> #include <lac/solver_cg.h> #include <lac/precondition.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/grid_out.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/grid_refinement.h> #include <dofs/dof_handler.h> #include <dofs/dof_constraints.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <fe/fe_tools.h> #include <numerics/vectors.h> #include <numerics/matrices.h> #include <numerics/data_out.h> #include <numerics/error_estimator.h> #include <iostream> #include <fstream> #include <list> #include <algorithm> #include <numeric> #include <sstream>
The last step is as in all previous programs:
using namespace dealii;
As mentioned in the introduction, significant parts of the program have simply been taken over from the step-13 example program. We therefore only comment on those things that are new.
First, the framework for evaluation of solutions is unchanged, i.e. the base class is the same, and the class to evaluate the solution at a grid point is unchanged:
namespace Evaluation
{
template <int dim> class EvaluationBase { public: virtual ~EvaluationBase (); void set_refinement_cycle (const unsigned int refinement_cycle); virtual void operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const = 0; protected: unsigned int refinement_cycle; }; template <int dim> EvaluationBase<dim>::~EvaluationBase () {} template <int dim> void EvaluationBase<dim>::set_refinement_cycle (const unsigned int step) { refinement_cycle = step; }
template <int dim> class PointValueEvaluation : public EvaluationBase<dim> { public: PointValueEvaluation (const Point<dim> &evaluation_point); virtual void operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const; DeclException1 (ExcEvaluationPointNotFound, Point<dim>, << "The evaluation point " << arg1 << " was not found among the vertices of the present grid."); private: const Point<dim> evaluation_point; }; template <int dim> PointValueEvaluation<dim>:: PointValueEvaluation (const Point<dim> &evaluation_point) : evaluation_point (evaluation_point) {} template <int dim> void PointValueEvaluation<dim>:: operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const { double point_value = 1e20; typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); bool evaluation_point_found = false; for (; (cell!=endc) && !evaluation_point_found; ++cell) for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex) if (cell->vertex(vertex).distance (evaluation_point) < cell->diameter() * 1e-8) { point_value = solution(cell->vertex_dof_index(vertex,0)); evaluation_point_found = true; break; }; AssertThrow (evaluation_point_found, ExcEvaluationPointNotFound(evaluation_point)); std::cout << " Point value=" << point_value << std::endl; }
Besides the class implementing the evaluation of the solution at one point, we here provide one which evaluates the gradient at a grid point. Since in general the gradient of a finite element function is not continuous at a vertex, we have to be a little bit more careful here. What we do is to loop over all cells, even if we have found the point already on one cell, and use the mean value of the gradient at the vertex taken from all adjacent cells.
Given the interface of the PointValueEvaluation class, the declaration of this class provides little surprise, and neither does the constructor:
template <int dim> class PointXDerivativeEvaluation : public EvaluationBase<dim> { public: PointXDerivativeEvaluation (const Point<dim> &evaluation_point); virtual void operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const; DeclException1 (ExcEvaluationPointNotFound, Point<dim>, << "The evaluation point " << arg1 << " was not found among the vertices of the present grid."); private: const Point<dim> evaluation_point; }; template <int dim> PointXDerivativeEvaluation<dim>:: PointXDerivativeEvaluation (const Point<dim> &evaluation_point) : evaluation_point (evaluation_point) {}
The more interesting things happen inside the function doing the actual evaluation:
template <int dim> void PointXDerivativeEvaluation<dim>:: operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const {
This time initialize the return value with something useful, since we will have to add up a number of contributions and take the mean value afterwards...
double point_derivative = 0;
...then have some objects of which the meaning wil become clear below...
QTrapez<dim> vertex_quadrature; FEValues<dim> fe_values (dof_handler.get_fe(), vertex_quadrature, update_gradients | update_quadrature_points); std::vector<Tensor<1,dim> > solution_gradients (vertex_quadrature.size());
...and next loop over all cells and their vertices, and count how often the vertex has been found:
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); unsigned int evaluation_point_hits = 0; for (; cell!=endc; ++cell) for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex) if (cell->vertex(vertex) == evaluation_point) {
Things are now no more as simple, since we can't get the gradient of the finite element field as before, where we simply had to pick one degree of freedom at a vertex.
Rather, we have to evaluate the finite element field on this cell, and at a certain point. As you know, evaluating finite element fields at certain points is done through the FEValues class, so we use that. The question is: the FEValues object needs to be a given a quadrature formula and can then compute the values of finite element quantities at the quadrature points. Here, we don't want to do quadrature, we simply want to specify some points!
Nevertheless, the same way is chosen: use a special quadrature rule with points at the vertices, since these are what we are interested in. The appropriate rule is the trapezoidal rule, so that is the reason why we used that one above.
Thus: initialize the FEValues object on this cell,
fe_values.reinit (cell);
and extract the gradients of the solution vector at the vertices:
fe_values.get_function_grads (solution,
solution_gradients);
Now we have the gradients at all vertices, so pick out that one which belongs to the evaluation point (note that the order of vertices is not necessarily the same as that of the quadrature points):
unsigned int q_point = 0; for (; q_point<solution_gradients.size(); ++q_point) if (fe_values.quadrature_point(q_point) == evaluation_point) break;
Check that the evaluation point was indeed found,
Assert (q_point < solution_gradients.size(), ExcInternalError());
and if so take the x-derivative of the gradient there as the value which we are interested in, and increase the counter indicating how often we have added to that variable:
point_derivative += solution_gradients[q_point][0];
++evaluation_point_hits;
Finally break out of the innermost loop iterating over the vertices of the present cell, since if we have found the evaluation point at one vertex it cannot be at a following vertex as well:
break;
};
Now we have looped over all cells and vertices, so check whether the point was found:
AssertThrow (evaluation_point_hits > 0, ExcEvaluationPointNotFound(evaluation_point));
We have simply summed up the contributions of all adjacent cells, so we still have to compute the mean value. Once this is done, report the status:
point_derivative /= evaluation_point_hits;
std::cout << " Point x-derivative=" << point_derivative
<< std::endl;
}
Since this program has a more difficult structure (it computed a dual solution in addition to a primal one), writing out the solution is no more done by an evaluation object since we want to write both solutions at once into one file, and that requires some more information than available to the evaluation classes.
However, we also want to look at the grids generated. This again can be done with one such class. Its structure is analog to the SolutionOutput class of the previous example program, so we do not discuss it here in more detail. Furthermore, everything that is used here has already been used in previous example programs.
template <int dim> class GridOutput : public EvaluationBase<dim> { public: GridOutput (const std::string &output_name_base); virtual void operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &solution) const; private: const std::string output_name_base; }; template <int dim> GridOutput<dim>:: GridOutput (const std::string &output_name_base) : output_name_base (output_name_base) {} template <int dim> void GridOutput<dim>::operator () (const DoFHandler<dim> &dof_handler, const Vector<double> &/ *solution* /) const { std::ostringstream filename; filename << output_name_base << "-" << this->refinement_cycle << ".eps" << std::ends; std::ofstream out (filename.str().c_str()); GridOut().write_eps (dof_handler.get_tria(), out); } }
Next are the actual solver classes. Again, we discuss only the differences to the previous program.
namespace LaplaceSolver
{
Before everything else, forward-declare one class that we will have later, since we will want to make it a friend of some of the classes that follow, which requires the class to be known:
template <int dim> class WeightedResidual;
This class is almost unchanged, with the exception that it declares two more functions: output_solution will be used to generate output files from the actual solutions computed by derived classes, and the set_refinement_cycle function by which the testing framework sets the number of the refinement cycle to a local variable in this class; this number is later used to generate filenames for the solution output.
template <int dim> class Base { public: Base (Triangulation<dim> &coarse_grid); virtual ~Base (); virtual void solve_problem () = 0; virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const = 0; virtual void refine_grid () = 0; virtual unsigned int n_dofs () const = 0; virtual void set_refinement_cycle (const unsigned int cycle); virtual void output_solution () const = 0; protected: const SmartPointer<Triangulation<dim> > triangulation; unsigned int refinement_cycle; }; template <int dim> Base<dim>::Base (Triangulation<dim> &coarse_grid) : triangulation (&coarse_grid) {} template <int dim> Base<dim>::~Base () {} template <int dim> void Base<dim>::set_refinement_cycle (const unsigned int cycle) { refinement_cycle = cycle; }
Likewise, the Solver class is entirely unchanged and will thus not be discussed.
template <int dim> class Solver : public virtual Base<dim> { public: Solver (Triangulation<dim> &triangulation, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &boundary_values); virtual ~Solver (); virtual void solve_problem (); virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const; virtual unsigned int n_dofs () const; protected: const SmartPointer<const FiniteElement<dim> > fe; const SmartPointer<const Quadrature<dim> > quadrature; const SmartPointer<const Quadrature<dim-1> > face_quadrature; DoFHandler<dim> dof_handler; Vector<double> solution; const SmartPointer<const Function<dim> > boundary_values; virtual void assemble_rhs (Vector<double> &rhs) const = 0; private: struct LinearSystem { LinearSystem (const DoFHandler<dim> &dof_handler); void solve (Vector<double> &solution) const; ConstraintMatrix hanging_node_constraints; SparsityPattern sparsity_pattern; SparseMatrix<double> matrix; Vector<double> rhs; }; void assemble_linear_system (LinearSystem &linear_system); void assemble_matrix (LinearSystem &linear_system, const typename DoFHandler<dim>::active_cell_iterator &begin_cell, const typename DoFHandler<dim>::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const; }; template <int dim> Solver<dim>::Solver (Triangulation<dim> &triangulation, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &boundary_values) : Base<dim> (triangulation), fe (&fe), quadrature (&quadrature), face_quadrature (&face_quadrature), dof_handler (triangulation), boundary_values (&boundary_values) {} template <int dim> Solver<dim>::~Solver () { dof_handler.clear (); } template <int dim> void Solver<dim>::solve_problem () { dof_handler.distribute_dofs (*fe); solution.reinit (dof_handler.n_dofs()); LinearSystem linear_system (dof_handler); assemble_linear_system (linear_system); linear_system.solve (solution); } template <int dim> void Solver<dim>:: postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const { postprocessor (dof_handler, solution); } template <int dim> unsigned int Solver<dim>::n_dofs () const { return dof_handler.n_dofs(); } template <int dim> void Solver<dim>::assemble_linear_system (LinearSystem &linear_system) { typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator; const unsigned int n_threads = multithread_info.n_default_threads; std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (), dof_handler.end (), n_threads); Threads::ThreadMutex mutex; Threads::ThreadGroup<> threads; for (unsigned int thread=0; thread<n_threads; ++thread) threads += Threads::spawn (*this, &Solver<dim>::assemble_matrix) (linear_system, thread_ranges[thread].first, thread_ranges[thread].second, mutex); assemble_rhs (linear_system.rhs); linear_system.hanging_node_constraints.condense (linear_system.rhs); std::map<unsigned int,double> boundary_value_map; VectorTools::interpolate_boundary_values (dof_handler, 0, *boundary_values, boundary_value_map); threads.join_all (); linear_system.hanging_node_constraints.condense (linear_system.matrix); MatrixTools::apply_boundary_values (boundary_value_map, linear_system.matrix, solution, linear_system.rhs); } template <int dim> void Solver<dim>::assemble_matrix (LinearSystem &linear_system, const typename DoFHandler<dim>::active_cell_iterator &begin_cell, const typename DoFHandler<dim>::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const { FEValues<dim> fe_values (*fe, *quadrature, update_gradients | update_JxW_values); const unsigned int dofs_per_cell = fe->dofs_per_cell; const unsigned int n_q_points = quadrature->size(); FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell; cell!=end_cell; ++cell) { cell_matrix = 0; fe_values.reinit (cell); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) * fe_values.shape_grad(j,q_point) * fe_values.JxW(q_point)); cell->get_dof_indices (local_dof_indices); Threads::ThreadMutex::ScopedLock lock (mutex); for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=0; j<dofs_per_cell; ++j) linear_system.matrix.add (local_dof_indices[i], local_dof_indices[j], cell_matrix(i,j)); }; } template <int dim> Solver<dim>::LinearSystem:: LinearSystem (const DoFHandler<dim> &dof_handler) { hanging_node_constraints.clear (); void (*mhnc_p) (const DoFHandler<dim> &, ConstraintMatrix &) = &DoFTools::make_hanging_node_constraints; Threads::Thread<> mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints); sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); mhnc_thread.join (); hanging_node_constraints.close (); hanging_node_constraints.condense (sparsity_pattern); sparsity_pattern.compress(); matrix.reinit (sparsity_pattern); rhs.reinit (dof_handler.n_dofs()); } template <int dim> void Solver<dim>::LinearSystem::solve (Vector<double> &solution) const { SolverControl solver_control (5000, 1e-12); SolverCG<> cg (solver_control); PreconditionSSOR<> preconditioner; preconditioner.initialize(matrix, 1.2); cg.solve (matrix, solution, rhs, preconditioner); hanging_node_constraints.distribute (solution); }
The PrimalSolver class is also mostly unchanged except for overloading the functions solve_problem, n_dofs, and postprocess of the base class, and implementing the output_solution function. These overloaded functions do nothing particular besides calling the functions of the base class -- that seems superfluous, but works around a bug in a popular compiler which requires us to write such functions for the following scenario: Besides the PrimalSolver class, we will have a DualSolver, both derived from Solver. We will then have a final classes which derived from these two, which will then have two instances of the Solver class as its base classes. If we want, for example, the number of degrees of freedom of the primal solver, we would have to indicate this like so: PrimalSolver<dim>n_dofs(). However, the compiler does not accept this since the n_dofs function is actually from a base class of the PrimalSolver class, so we have to inject the name from the base to the derived class using these additional functions.
Regarding the implementation of the output_solution function, we keep the GlobalRefinement and RefinementKelly classes in this program, and they can then rely on the default implementation of this function which simply outputs the primal solution. The class implementing dual weighted error estimators will overload this function itself, to also output the dual solution.
Except for this, the class is unchanged with respect to the previous example.
template <int dim> class PrimalSolver : public Solver<dim> { public: PrimalSolver (Triangulation<dim> &triangulation, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values); virtual void solve_problem (); virtual unsigned int n_dofs () const; virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const; virtual void output_solution () const; protected: const SmartPointer<const Function<dim> > rhs_function; virtual void assemble_rhs (Vector<double> &rhs) const;
Now, in order to work around some problems in one of the compilers this library can be compiled with, we will have to declare a class that is actually derived from the present one, as a friend (strange as that seems). The full rationale will be explained below.
friend class WeightedResidual<dim>; }; template <int dim> PrimalSolver<dim>:: PrimalSolver (Triangulation<dim> &triangulation, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values) : Base<dim> (triangulation), Solver<dim> (triangulation, fe, quadrature, face_quadrature, boundary_values), rhs_function (&rhs_function) {} template <int dim> void PrimalSolver<dim>::solve_problem () { Solver<dim>::solve_problem (); } template <int dim> unsigned int PrimalSolver<dim>::n_dofs() const { return Solver<dim>::n_dofs(); } template <int dim> void PrimalSolver<dim>:: postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const { Solver<dim>::postprocess(postprocessor); } template <int dim> void PrimalSolver<dim>::output_solution () const { DataOut<dim> data_out; data_out.attach_dof_handler (this->dof_handler); data_out.add_data_vector (this->solution, "solution"); data_out.build_patches (); std::ostringstream filename; filename << "solution-" << this->refinement_cycle << ".gnuplot" << std::ends; std::ofstream out (filename.str().c_str()); data_out.write (out, DataOut<dim>::gnuplot); } template <int dim> void PrimalSolver<dim>:: assemble_rhs (Vector<double> &rhs) const { FEValues<dim> fe_values (*this->fe, *this->quadrature, update_values | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = this->fe->dofs_per_cell; const unsigned int n_q_points = this->quadrature->size(); Vector<double> cell_rhs (dofs_per_cell); std::vector<double> rhs_values (n_q_points); std::vector<unsigned int> local_dof_indices (dofs_per_cell); typename DoFHandler<dim>::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (; cell!=endc; ++cell) { cell_rhs = 0; fe_values.reinit (cell); rhs_function->value_list (fe_values.get_quadrature_points(), rhs_values); for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) cell_rhs(i) += (fe_values.shape_value(i,q_point) * rhs_values[q_point] * fe_values.JxW(q_point)); cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i<dofs_per_cell; ++i) rhs(local_dof_indices[i]) += cell_rhs(i); } }
For the following two classes, the same applies as for most of the above: the class is taken from the previous example as-is:
template <int dim> class RefinementGlobal : public PrimalSolver<dim> { public: RefinementGlobal (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values); virtual void refine_grid (); }; template <int dim> RefinementGlobal<dim>:: RefinementGlobal (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values) : Base<dim> (coarse_grid), PrimalSolver<dim> (coarse_grid, fe, quadrature, face_quadrature, rhs_function, boundary_values) {} template <int dim> void RefinementGlobal<dim>::refine_grid () { this->triangulation->refine_global (1); } template <int dim> class RefinementKelly : public PrimalSolver<dim> { public: RefinementKelly (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values); virtual void refine_grid (); }; template <int dim> RefinementKelly<dim>:: RefinementKelly (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values) : Base<dim> (coarse_grid), PrimalSolver<dim> (coarse_grid, fe, quadrature, face_quadrature, rhs_function, boundary_values) {} template <int dim> void RefinementKelly<dim>::refine_grid () { Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells()); KellyErrorEstimator<dim>::estimate (this->dof_handler, QGauss<dim-1>(3), typename FunctionMap<dim>::type(), this->solution, estimated_error_per_cell); GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation, estimated_error_per_cell, 0.3, 0.03); this->triangulation->execute_coarsening_and_refinement (); }
This class is a variant of the previous one, in that it allows to weight the refinement indicators we get from the library's Kelly indicator by some function. We include this class since the goal of this example program is to demonstrate automatic refinement criteria even for complex output quantities such as point values or stresses. If we did not solve a dual problem and compute the weights thereof, we would probably be tempted to give a hand-crafted weighting to the indicators to account for the fact that we are going to evaluate these quantities. This class accepts such a weighting function as argument to its constructor:
template <int dim> class RefinementWeightedKelly : public PrimalSolver<dim> { public: RefinementWeightedKelly (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values, const Function<dim> &weighting_function); virtual void refine_grid (); private: const SmartPointer<const Function<dim> > weighting_function; }; template <int dim> RefinementWeightedKelly<dim>:: RefinementWeightedKelly (Triangulation<dim> &coarse_grid, const FiniteElement<dim> &fe, const Quadrature<dim> &quadrature, const Quadrature<dim-1> &face_quadrature, const Function<dim> &rhs_function, const Function<dim> &boundary_values, const Function<dim> &weighting_function) : Base<dim> (coarse_grid), PrimalSolver<dim> (coarse_grid, fe, quadrature, face_quadrature, rhs_function, boundary_values), weighting_function (&weighting_function) {}
Now, here comes the main function, including the weighting:
template <int dim> void RefinementWeightedKelly<dim>::refine_grid () {
First compute some residual based error indicators for all cells by a method already implemented in the library. What exactly is computed can be read in the documentation of that class.
Vector<float> estimated_error (this->triangulation->n_active_cells()); KellyErrorEstimator<dim>::estimate (this->dof_handler, *this->face_quadrature, typename FunctionMap<dim>::type(), this->solution, estimated_error);
Now we are going to weight these indicators by the value of the function given to the constructor:
typename DoFHandler<dim>::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index) estimated_error(cell_index) *= weighting_function->value (cell->center()); GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation, estimated_error, 0.3, 0.03); this->triangulation->execute_coarsening_and_refinement (); } }
In this example program, we work with the same data sets as in the previous one, but as it may so happen that someone wants to run the program with different boundary values and right hand side functions, or on a different grid, we show a simple technique to do exactly that. For more clarity, we furthermore pack everything that has to do with equation data into a namespace of its own.
The underlying assumption is that this is a research program, and that there we often have a number of test cases that consist of a domain, a right hand side, boundary values, possibly a specified coefficient, and a number of other parameters. They often vary all at the same time when shifting from one example to another. To make handling such sets of problem description parameters simple is the goal of the following.
Basically, the idea is this: let us have a structure for each set of data, in which we pack everything that describes a test case: here, these are two subclasses, one called BoundaryValues for the boundary values of the exact solution, and one called RightHandSide, and then a way to generate the coarse grid. Since the solution of the previous example program looked like curved ridges, we use this name here for the enclosing class. Note that the names of the two inner classes have to be the same for all enclosing test case classes, and also that we have attached the dimension template argument to the enclosing class rather than to the inner ones, to make further processing simpler. (From a language viewpoint, a namespace would be better to encapsulate these inner classes, rather than a structure. However, namespaces cannot be given as template arguments, so we use a structure to allow a second object to select from within its given argument. The enclosing structure, of course, has no member variables apart from the classes it declares, and a static function to generate the coarse mesh; it will in general never be instantiated.)
The idea is then the following (this is the right time to also take a brief look at the code below): we can generate objects for boundary values and right hand side by simply giving the name of the outer class as a template argument to a class which we call here Data::SetUp, and it then creates objects for the inner classes. In this case, to get all that characterizes the curved ridge solution, we would simply generate an instance of Data::SetUp<Data::CurvedRidge>, and everything we need to know about the solution would be static member variables and functions of that object.
This approach might seem like overkill in this case, but will become very handy once a certain set up is not only characterized by Dirichlet boundary values and a right hand side function, but in addition by material properties, Neumann values, different boundary descriptors, etc. In that case, the SetUp class might consist of a dozen or more objects, and each descriptor class (like the CurvedRidges class below) would have to provide them. Then, you will be happy to be able to change from one set of data to another by only changing the template argument to the SetUp class at one place, rather than at many.
With this framework for different test cases, we are almost finished, but one thing remains: by now we can select statically, by changing one template argument, which data set to choose. In order to be able to do that dynamically, i.e. at run time, we need a base class. This we provide in the obvious way, see below, with virtual abstract functions. It forces us to introduce a second template parameter dim which we need for the base class (which could be avoided using some template magic, but we omit that), but that's all.
Adding new testcases is now simple, you don't have to touch the framework classes, only a structure like the CurvedRidges one is needed.
namespace Data
{
Based on the above description, the SetUpBase class then looks as follows. To allow using the SmartPointer class with this class, we derived from the Subscriptor class.
template <int dim> struct SetUpBase : public Subscriptor { virtual const Function<dim> & get_boundary_values () const = 0; virtual const Function<dim> & get_right_hand_side () const = 0; virtual void create_coarse_grid (Triangulation<dim> &coarse_grid) const = 0; };
And now for the derived class that takes the template argument as explained above. For some reason, C++ requires us to define a constructor (which maybe empty), as otherwise a warning is generated that some data is not initialized.
Here we pack the data elements into private variables, and allow access to them through the methods of the base class.
template <class Traits, int dim> struct SetUp : public SetUpBase<dim> { SetUp () {} virtual const Function<dim> & get_boundary_values () const; virtual const Function<dim> & get_right_hand_side () const; virtual void create_coarse_grid (Triangulation<dim> &coarse_grid) const; private: static const typename Traits::BoundaryValues boundary_values; static const typename Traits::RightHandSide right_hand_side; };
We have to provide definitions for the static member variables of the above class:
template <class Traits, int dim> const typename Traits::BoundaryValues SetUp<Traits,dim>::boundary_values; template <class Traits, int dim> const typename Traits::RightHandSide SetUp<Traits,dim>::right_hand_side;
And definitions of the member functions:
template <class Traits, int dim> const Function<dim> & SetUp<Traits,dim>::get_boundary_values () const { return boundary_values; } template <class Traits, int dim> const Function<dim> & SetUp<Traits,dim>::get_right_hand_side () const { return right_hand_side; } template <class Traits, int dim> void SetUp<Traits,dim>:: create_coarse_grid (Triangulation<dim> &coarse_grid) const { Traits::create_coarse_grid (coarse_grid); }
The class that is used to describe the boundary values and right hand side of the curved ridge problem already used in the step-13 example program is then like so:
template <int dim> struct CurvedRidges { class BoundaryValues : public Function<dim> { public: BoundaryValues () : Function<dim> () {} virtual double value (const Point<dim> &p, const unsigned int component) const; }; class RightHandSide : public Function<dim> { public: RightHandSide () : Function<dim> () {} virtual double value (const Point<dim> &p, const unsigned int component) const; }; static void create_coarse_grid (Triangulation<dim> &coarse_grid); }; template <int dim> double CurvedRidges<dim>::BoundaryValues:: value (const Point<dim> &p, const unsigned int / *component* /) const { double q = p(0); for (unsigned int i=1; i<dim; ++i) q += std::sin(10*p(i)+5*p(0)*p(0)); const double exponential = std::exp(q); return exponential; } template <int dim> double CurvedRidges<dim>::RightHandSide::value (const Point<dim> &p, const unsigned int / *component* /) const { double q = p(0); for (unsigned int i=1; i<dim; ++i) q += std::sin(10*p(i)+5*p(0)*p(0)); const double u = std::exp(q); double t1 = 1, t2 = 0, t3 = 0; for (unsigned int i=1; i<dim; ++i) { t1 += std::cos(10*p(i)+5*p(0)*p(0)) * 10 * p(0); t2 += 10*std::cos(10*p(i)+5*p(0)*p(0)) - 100*std::sin(10*p(i)+5*p(0)*p(0)) * p(0)*p(0); t3 += 100*std::cos(10*p(i)+5*p(0)*p(0))*std::cos(10*p(i)+5*p(0)*p(0)) - 100*std::sin(10*p(i)+5*p(0)*p(0)); }; t1 = t1*t1; return -u*(t1+t2+t3); } template <int dim> void CurvedRidges<dim>:: create_coarse_grid (Triangulation<dim> &coarse_grid) { GridGenerator::hyper_cube (coarse_grid, -1, 1); coarse_grid.refine_global (2); }
This example program was written while giving practical courses for a lecture on adaptive finite element methods and duality based error estimates. For these courses, we had one exercise, which required to solve the Laplace equation with constant right hand side on a square domain with a square hole in the center, and zero boundary values. Since the implementation of the properties of this problem is so particularly simple here, lets do it. As the number of the exercise was 2.3, we take the liberty to retain this name for the class as well.
template <int dim> struct Exercise_2_3 {
We need a class to denote the boundary values of the problem. In this case, this is simple: it's the zero function, so don't even declare a class, just a typedef:
typedef ZeroFunction<dim> BoundaryValues;
Second, a class that denotes the right hand side. Since they are constant, just subclass the corresponding class of the library and be done:
class RightHandSide : public ConstantFunction<dim> { public: RightHandSide () : ConstantFunction<dim> (1.) {} };
Finally a function to generate the coarse grid. This is somewhat more complicated here, see immediately below.
static void create_coarse_grid (Triangulation<dim> &coarse_grid); };
As stated above, the grid for this example is the square [-1,1]^2 with the square [-1/2,1/2]^2 as hole in it. We create the coarse grid as 4 times 4 cells with the middle four ones missing.
Of course, the example has an extension to 3d, but since this function cannot be written in a dimension independent way we choose not to implement this here, but rather only specialize the template for dim=2. If you compile the program for 3d, you'll get a message from the linker that this function is not implemented for 3d, and needs to be provided.
For the creation of this geometry, the library has no predefined method. In this case, the geometry is still simple enough to do the creation by hand, rather than using a mesh generator.
template <> void Exercise_2_3<2>:: create_coarse_grid (Triangulation<2> &coarse_grid) {
First define the space dimension, to allow those parts of the function that are actually dimension independent to use this variable. That makes it simpler if you later takes this as a starting point to implement the 3d version.
const unsigned int dim = 2;
Then have a list of vertices. Here, they are 24 (5 times 5, with the middle one omitted). It is probably best to draw a sketch here. Note that we leave the number of vertices open at first, but then let the compiler compute this number afterwards. This reduces the possibility of having the dimension to large and leaving the last ones uninitialized.
static const Point<2> vertices_1[] = { Point<2> (-1., -1.), Point<2> (-1./2, -1.), Point<2> (0., -1.), Point<2> (+1./2, -1.), Point<2> (+1, -1.), Point<2> (-1., -1./2.), Point<2> (-1./2, -1./2.), Point<2> (0., -1./2.), Point<2> (+1./2, -1./2.), Point<2> (+1, -1./2.), Point<2> (-1., 0.), Point<2> (-1./2, 0.), Point<2> (+1./2, 0.), Point<2> (+1, 0.), Point<2> (-1., 1./2.), Point<2> (-1./2, 1./2.), Point<2> (0., 1./2.), Point<2> (+1./2, 1./2.), Point<2> (+1, 1./2.), Point<2> (-1., 1.), Point<2> (-1./2, 1.), Point<2> (0., 1.), Point<2> (+1./2, 1.), Point<2> (+1, 1.) }; const unsigned int n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
From this static list of vertices, we generate an STL vector of the vertices, as this is the data type the library wants to see.
const std::vector<Point<dim> > vertices (&vertices_1[0],
&vertices_1[n_vertices]);
Next, we have to define the cells and the vertices they contain. Here, we have 8 vertices, but leave the number open and let it be computed afterwards:
static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell] = {{0, 1, 5, 6}, {1, 2, 6, 7}, {2, 3, 7, 8}, {3, 4, 8, 9}, {5, 6, 10, 11}, {8, 9, 12, 13}, {10, 11, 14, 15}, {12, 13, 17, 18}, {14, 15, 19, 20}, {15, 16, 20, 21}, {16, 17, 21, 22}, {17, 18, 22, 23}}; const unsigned int n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
Again, we generate a C++ vector type from this, but this time by looping over the cells (yes, this is boring). Additionally, we set the material indicator to zero for all the cells:
std::vector<CellData<dim> > cells (n_cells, CellData<dim>()); for (unsigned int i=0; i<n_cells; ++i) { for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j) cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; };
Finally pass all this information to the library to generate a triangulation. The last parameter may be used to pass information about non-zero boundary indicators at certain faces of the triangulation to the library, but we don't want that here, so we give an empty object:
coarse_grid.create_triangulation (vertices,
cells,
SubCellData());
And since we want that the evaluation point (3/4,3/4) in this example is a grid point, we refine once globally:
coarse_grid.refine_global (1); } }
As you have now read through this framework, you may be wondering why we have not chosen to implement the classes implementing a certain setup (like the CurvedRidges class) directly as classes derived from Data::SetUpBase. Indeed, we could have done very well so. The only reason is that then we would have to have member variables for the solution and right hand side classes in the CurvedRidges class, as well as member functions overloading the abstract functions of the base class giving access to these member variables. The SetUp class has the sole reason to relieve us from the need to reiterate these member variables and functions that would be necessary in all such classes. In some way, the template mechanism here only provides a way to have default implementations for a number of functions that depend on external quantities and can thus not be provided using normal virtual functions, at least not without the help of templates.
However, there might be good reasons to actually implement classes derived from Data::SetUpBase, for example if the solution or right hand side classes require constructors that take arguments, which the Data::SetUpBase class cannot provide. In that case, subclassing is a worthwhile strategy. Other possibilities for special cases are to derive from Data::SetUp<SomeSetUp> where SomeSetUp denotes a class, or even to explicitly specialize Data::SetUp<SomeSetUp>. The latter allows to transparently use the way the SetUp class is used for other set-ups, but with special actions taken for special arguments.
A final observation favoring the approach taken here is the following: we have found numerous times that when starting a project, the number of parameters (usually boundary values, right hand side, coarse grid, just as here) was small, and the number of test cases was small as well. One then starts out by handcoding them into a number of switch statements. Over time, projects grow, and so does the number of test cases. The number of switch statements grows with that, and their length as well, and one starts to find ways to consider impossible examples where domains, boundary values, and right hand sides do not fit together any more, and starts loosing the overview over the whole structure. Encapsulating everything belonging to a certain test case into a structure of its own has proven worthwhile for this, as it keeps everything that belongs to one test case in one place. Furthermore, it allows to put these things all in one or more files that are only devoted to test cases and their data, without having to bring their actual implementation into contact with the rest of the program.
As with the other components of the program, we put everything we need to describe dual functionals into a namespace of its own, and define an abstract base class that provides the interface the class solving the dual problem needs for its work.
We will then implement two such classes, for the evaluation of a point value and of the derivative of the solution at that point. For these functionals we already have the corresponding evaluation objects, so they are comlementary.
namespace DualFunctional
{
First start with the base class for dual functionals. Since for linear problems the characteristics of the dual problem play a role only in the right hand side, we only need to provide for a function that assembles the right hand side for a given discretization:
template <int dim> class DualFunctionalBase : public Subscriptor { public: virtual void assemble_rhs (const DoFHandler<dim> &dof_handler, Vector<double> &rhs) const = 0; };
As a first application, we consider the functional corresponding to the evaluation of the solution's value at a given point which again we assume to be a vertex. Apart from the constructor that takes and stores the evaluation point, this class consists only of the function that implements assembling the right hand side.
template <int dim> class PointValueEvaluation : public DualFunctionalBase<dim> { public: PointValueEvaluation (const Point<dim> &evaluation_point); virtual void assemble_rhs (const DoFHandler<dim> &dof_handler, Vector<double> &rhs) const; DeclException1 (ExcEvaluationPointNotFound, Point<dim>, << "The evaluation point " << arg1 << " was not found among the vertices of the present grid."); protected: const Point<dim> evaluation_point; }; template <int dim> PointValueEvaluation<dim>:: PointValueEvaluation (const Point<dim> &evaluation_point) : evaluation_point (evaluation_point) {}
As for doing the main purpose of the class, assembling the right hand side, let us first consider what is necessary: The right hand side of the dual problem is a vector of values J(phi_i), where J is the error functional, and phi_i is the i-th shape function. Here, J is the evaluation at the point x0, i.e. J(phi_i)=phi_i(x0).
Now, we have assumed that the evaluation point is a vertex. Thus, for the usual finite elements we might be using in this program, we can take for granted that at such a point exactly one shape function is nonzero, and in particular has the value one. Thus, we set the right hand side vector to all-zeros, then seek for the shape function associated with that point and set the corresponding value of the right hand side vector to one:
template <int dim> void PointValueEvaluation<dim>:: assemble_rhs (const DoFHandler<dim> &dof_handler, Vector<double> &rhs) const