--- /dev/null
+/* $Id: step-14.cc 2008 2006-01-23 23:53:16Z wolf $ */
+/* Author: Wolfgang Bangerth, ETH Zurich, 2002 */
+
+/* $Id: step-14.cc 2008 2006-01-23 23:53:16Z wolf $ */
+/* Version: $Name$ */
+/* */
+/* Copyright (C) 2002, 2003, 2004, 2008 by the deal.II authors */
+/* */
+/* This file is subject to QPL and may not be distributed */
+/* without copyright and license information. Please refer */
+/* to the file deal.II/doc/license.html for the text and */
+/* further information on this license. */
+
+
+ // 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/compressed_sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/vector_memory.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 <grid/tria_boundary_lib.h>
+#include <grid/intergrid_map.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 <fe/mapping_q.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+#include <numerics/derivative_approximation.h>
+
+#include <iostream>
+#include <fstream>
+#include <list>
+#include <algorithm>
+#include <numeric>
+#include <stdlib.h>
+
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
+
+using namespace dealii;
+
+
+/**
+ * The number of iterations the program shall do. This is given as a
+ * parameter to the executable.
+ */
+unsigned int n_steps = 0;
+/**
+ * The present step.
+ */
+unsigned int step;
+
+ /**
+ * Declare the coefficient in front
+ * of the Laplace operator of the
+ * equation.
+ */
+template <int dim>
+class LaplaceCoefficient : public Function<dim>
+{
+ public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ return 1.+p*p;
+ }
+};
+
+
+ /**
+ * Same for the coefficient in front
+ * of the mass term.
+ */
+template <int dim>
+class MassCoefficient : public Function<dim>
+{
+ public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ double factor = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ factor *= std::sin(4*p(d));
+ return 1.+factor/2;
+ }
+};
+
+
+
+ // @sect3{Evaluating the solution}
+
+ // 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
+{
+ // @sect4{The EvaluationBase class}
+ 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;
+
+ 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;
+ }
+
+
+ // @sect4{The PointValueEvaluation class}
+ 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.");
+
+ 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;
+ }
+
+
+ // @sect4{The PointXDerivativeEvaluation class}
+
+ // 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.");
+
+ 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;
+ MappingQ<dim> mapping (4);
+ FEValues<dim> fe_values (mapping, dof_handler.get_fe(),
+ vertex_quadrature,
+ update_gradients | update_q_points);
+ std::vector<Tensor<1,dim> >
+ solution_gradients (vertex_quadrature.n_quadrature_points);
+
+ // ...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;
+ }
+
+
+
+ // @sect4{The GridOutput class}
+
+ // 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;
+
+ 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
+ {
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ostringstream filename;
+#else
+ std::ostrstream filename;
+#endif
+ filename << "spec2006-447.dealII/"
+ << output_name_base << "-"
+ << this->refinement_cycle
+ << ".eps"
+ << std::ends;
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ofstream out (filename.str().c_str());
+#else
+ std::ofstream out (filename.str());
+#endif
+
+ GridOut().write_eps (dof_handler.get_tria(), out);
+ }
+}
+
+
+ // @sect3{The Laplace solver classes}
+
+ // 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;
+
+
+ // @sect4{The Laplace solver base class}
+
+ // 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;
+
+
+ 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;
+ }
+
+
+ // @sect4{The Laplace Solver class}
+
+ // 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;
+
+
+ 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;
+
+
+ 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
+ {
+ MappingQ<dim> mapping (4);
+ FEValues<dim> fe_values (mapping, *fe, *quadrature,
+ UpdateFlags(update_gradients | update_values |
+ update_q_points |
+ update_JxW_values));
+
+ const unsigned int dofs_per_cell = fe->dofs_per_cell;
+ const unsigned int n_q_points = quadrature->n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ std::vector<double> laplace_coefficients (fe_values.n_quadrature_points);
+ std::vector<double> mass_coefficients (fe_values.n_quadrature_points);
+
+
+ for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+ cell!=end_cell; ++cell)
+ {
+ cell_matrix = 0;
+
+ fe_values.reinit (cell);
+
+ LaplaceCoefficient<dim>().value_list (fe_values.get_quadrature_points(),
+ laplace_coefficients);
+ MassCoefficient<dim>().value_list (fe_values.get_quadrature_points(),
+ mass_coefficients);
+
+
+ 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) *
+ laplace_coefficients[q_point]
+ +
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point) *
+ mass_coefficients[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);
+
+ // make sparsity pattern. since
+ // in 3d the usual way just blows
+ // the roof w.r.t. memory
+ // consumption, use the detour
+ // via a compressed sparsity
+ // pattern that we later copy
+ // over
+ CompressedSparsityPattern csp(dof_handler.n_dofs(),
+ dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, csp);
+
+ mhnc_thread.join ();
+ hanging_node_constraints.close ();
+ hanging_node_constraints.condense (csp);
+
+ sparsity_pattern.copy_from(csp);
+ 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 (solution.size(), 1e-6);
+ PrimitiveVectorMemory<> vector_memory;
+ SolverCG<> cg (solver_control, vector_memory);
+
+ PreconditionJacobi<> preconditioner;
+ preconditioner.initialize(matrix);
+
+ cg.solve (matrix, solution, rhs, preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+ }
+
+
+
+
+ // @sect4{The PrimalSolver class}
+
+ // 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;
+
+
+ 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 use some
+ // workarounds. This will
+ // require that we 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
+ {
+ abort ();
+ }
+
+
+
+ template <int dim>
+ void
+ PrimalSolver<dim>::
+ assemble_rhs (Vector<double> &rhs) const
+ {
+ MappingQ<dim> mapping (4);
+ FEValues<dim> fe_values (mapping, *this->fe, *this->quadrature,
+ UpdateFlags(update_values |
+ update_q_points |
+ update_JxW_values));
+
+ const unsigned int dofs_per_cell = this->fe->dofs_per_cell;
+ const unsigned int n_q_points = this->quadrature->n_quadrature_points;
+
+ 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);
+ }
+ }
+
+
+ // @sect4{The RefinementGlobal and RefinementKelly classes}
+
+ // 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,
+ QGauss3<dim-1>(),
+ typename FunctionMap<dim>::type(),
+ this->solution,
+ estimated_error_per_cell);
+ GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
+ estimated_error_per_cell,
+ 0.2, 0.02);
+ this->triangulation->execute_coarsening_and_refinement ();
+ }
+
+
+
+ // @sect4{The RefinementWeightedKelly class}
+
+ // 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 ();
+
+
+ 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.2, 0.02);
+ this->triangulation->execute_coarsening_and_refinement ();
+ }
+
+}
+
+
+ // @sect3{Equation data}
+ //
+ // 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
+{
+ // @sect4{The SetUpBase and SetUp classes}
+
+ // 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;
+
+
+ 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);
+ }
+
+
+ // @sect4{The CurvedRidges class}
+
+ // 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);
+ }
+
+
+ // @sect4{The Exercise_2_3 class}
+
+ // 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, 6,5},
+ {1, 2, 7, 6},
+ {2, 3, 8, 7},
+ {3, 4, 9, 8},
+ {5, 6, 11, 10},
+ {8, 9, 13, 12},
+ {10, 11, 15, 14},
+ {12, 13, 18, 17},
+ {14, 15, 20, 19},
+ {15, 16, 21, 20},
+ {16, 17, 22, 21},
+ {17, 18, 23, 22}};
+ 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);
+ }
+*/
+
+
+ template <>
+ void
+ Exercise_2_3<3>::
+ create_coarse_grid (Triangulation<3> &coarse_grid)
+ {
+ GridGenerator::hyper_ball (coarse_grid);
+ static HyperBallBoundary<3> boundary;
+ coarse_grid.set_boundary (0, boundary);
+ coarse_grid.refine_global (1);
+ }
+
+}
+
+ // @sect4{Discussion}
+ //
+ // 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.
+
+
+ // @sect3{Dual functionals}
+
+ // 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
+{
+ // @sect4{The DualFunctionalBase class}
+
+ // 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;
+ };
+
+
+ // @sect4{The PointValueEvaluation class}
+
+ // 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.");
+
+
+ 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
+ {
+ // So, first set everything to
+ // zeros...
+ rhs.reinit (dof_handler.n_dofs());
+
+ // ...then loop over cells and
+ // find the evaluation point
+ // among the vertices (or very
+ // close to a vertex, which may
+ // happen due to floating point
+ // round-off):
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int vertex=0;
+ vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ if (cell->vertex(vertex).distance(evaluation_point)
+ < cell->diameter()*1e-8)
+ {
+ // Ok, found, so set
+ // corresponding entry,
+ // and leave function
+ // since we are finished:
+ rhs(cell->vertex_dof_index(vertex,0)) = 1;
+ return;
+ };
+
+ // Finally, a sanity check: if we
+ // somehow got here, then we must
+ // have missed the evaluation
+ // point, so raise an exception
+ // unconditionally:
+ AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point));
+ }
+
+
+ // @sect4{The PointXDerivativeEvaluation class}
+
+ // As second application, we again
+ // consider the evaluation of the
+ // x-derivative of the solution at
+ // one point. Again, the
+ // declaration of the class, and
+ // the implementation of its
+ // constructor is not too
+ // interesting:
+ template <int dim>
+ class PointXDerivativeEvaluation : public DualFunctionalBase<dim>
+ {
+ public:
+ PointXDerivativeEvaluation (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.");
+
+
+ const Point<dim> evaluation_point;
+ };
+
+
+ template <int dim>
+ PointXDerivativeEvaluation<dim>::
+ PointXDerivativeEvaluation (const Point<dim> &evaluation_point)
+ :
+ evaluation_point (evaluation_point)
+ {}
+
+
+ // What is interesting is the
+ // implementation of this
+ // functional: here,
+ // J(phi_i)=d/dx phi_i(x0).
+ //
+ // We could, as in the
+ // implementation of the respective
+ // evaluation object take the
+ // average of the gradients of each
+ // shape function phi_i at this
+ // evaluation point. However, we
+ // take a slightly different
+ // approach: we simply take the
+ // average over all cells that
+ // surround this point. The
+ // question which cells
+ // ``surrounds'' the evaluation
+ // point is made dependent on the
+ // mesh width by including those
+ // cells for which the distance of
+ // the cell's midpoint to the
+ // evaluation point is less than
+ // the cell's diameter.
+ //
+ // Taking the average of the
+ // gradient over the area/volume of
+ // these cells leads to a dual
+ // solution which is very close to
+ // the one which would result from
+ // the point evaluation of the
+ // gradient. It is simple to
+ // justify theoretically that this
+ // does not change the method
+ // significantly.
+ template <int dim>
+ void
+ PointXDerivativeEvaluation<dim>::
+ assemble_rhs (const DoFHandler<dim> &dof_handler,
+ Vector<double> &rhs) const
+ {
+ // Again, first set all entries
+ // to zero:
+ rhs.reinit (dof_handler.n_dofs());
+
+ // Initialize a ``FEValues''
+ // object with a quadrature
+ // formula, have abbreviations
+ // for the number of quadrature
+ // points and shape functions...
+ QGauss4<dim> quadrature;
+ MappingQ<dim> mapping (4);
+ FEValues<dim> fe_values (mapping, dof_handler.get_fe(), quadrature,
+ update_gradients |
+ update_q_points |
+ update_JxW_values);
+ const unsigned int n_q_points = fe_values.n_quadrature_points;
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+
+ // ...and have two objects that
+ // are used to store the global
+ // indices of the degrees of
+ // freedom on a cell, and the
+ // values of the gradients of the
+ // shape functions at the
+ // quadrature points:
+ Vector<double> cell_rhs (dofs_per_cell);
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ // Finally have a variable in
+ // which we will sum up the
+ // area/volume of the cells over
+ // which we integrate, by
+ // integrating the unit functions
+ // on these cells:
+ double total_volume = 0;
+
+ // Then start the loop over all
+ // cells, and select those cells
+ // which are close enough to the
+ // evaluation point:
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ if (cell->center().distance(evaluation_point) <=
+ cell->diameter())
+ {
+ // If we have found such a
+ // cell, then initialize
+ // the ``FEValues'' object
+ // and integrate the
+ // x-component of the
+ // gradient of each shape
+ // function, as well as the
+ // unit function for the
+ // total area/volume.
+ fe_values.reinit (cell);
+ cell_rhs = 0;
+
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
+ fe_values.JxW (q);
+ total_volume += fe_values.JxW (q);
+ };
+
+ // If we have the local
+ // contributions,
+ // distribute them to the
+ // global vector:
+ 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);
+ };
+
+ // After we have looped over all
+ // cells, check whether we have
+ // found any at all, by making
+ // sure that their volume is
+ // non-zero. If not, then the
+ // results will be botched, as
+ // the right hand side should
+ // then still be zero, so throw
+ // an exception:
+ AssertThrow (total_volume > 0,
+ ExcEvaluationPointNotFound(evaluation_point));
+
+ // Finally, we have by now only
+ // integrated the gradients of
+ // the shape functions, not
+ // taking their mean value. We
+ // fix this by dividing by the
+ // measure of the volume over
+ // which we have integrated:
+ rhs.scale (1./total_volume);
+ }
+
+
+}
+
+
+ // @sect3{Extending the LaplaceSolver namespace}
+namespace LaplaceSolver
+{
+
+ // @sect4{The DualSolver class}
+
+ // In the same way as the
+ // ``PrimalSolver'' class above, we
+ // now implement a
+ // ``DualSolver''. It has all the
+ // same features, the only
+ // difference is that it does not
+ // take a function object denoting
+ // a right hand side object, but
+ // now takes a
+ // ``DualFunctionalBase'' object
+ // that will assemble the right
+ // hand side vector of the dual
+ // problem. The rest of the class
+ // is rather trivial.
+ //
+ // Since both primal and dual
+ // solver will use the same
+ // triangulation, but different
+ // discretizations, it now becomes
+ // clear why we have made the
+ // ``Base'' class a virtual one:
+ // since the final class will be
+ // derived from both
+ // ``PrimalSolver'' as well as
+ // ``DualSolver'', it would have
+ // two ``Base'' instances, would we
+ // not have marked the inheritance
+ // as virtual. Since in many
+ // applications the base class
+ // would store much more
+ // information than just the
+ // triangulation which needs to be
+ // shared between primal and dual
+ // solvers, we do not usually want
+ // to use two such base classes.
+ template <int dim>
+ class DualSolver : public Solver<dim>
+ {
+ public:
+ DualSolver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Quadrature<dim-1> &face_quadrature,
+ const DualFunctional::DualFunctionalBase<dim> &dual_functional);
+
+ virtual
+ void
+ solve_problem ();
+
+ virtual
+ unsigned int
+ n_dofs () const;
+
+ virtual
+ void
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
+
+
+ const SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;
+ virtual void assemble_rhs (Vector<double> &rhs) const;
+
+ static const ZeroFunction<dim> boundary_values;
+
+ // Same as above -- make a
+ // derived class a friend of
+ // this one:
+ friend class WeightedResidual<dim>;
+ };
+
+ template <int dim>
+ const ZeroFunction<dim> DualSolver<dim>::boundary_values;
+
+ template <int dim>
+ DualSolver<dim>::
+ DualSolver (Triangulation<dim> &triangulation,
+ const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Quadrature<dim-1> &face_quadrature,
+ const DualFunctional::DualFunctionalBase<dim> &dual_functional)
+ :
+ Base<dim> (triangulation),
+ Solver<dim> (triangulation, fe,
+ quadrature, face_quadrature,
+ boundary_values),
+ dual_functional (&dual_functional)
+ {}
+
+
+ template <int dim>
+ void
+ DualSolver<dim>::solve_problem ()
+ {
+ Solver<dim>::solve_problem ();
+ }
+
+
+
+ template <int dim>
+ unsigned int
+ DualSolver<dim>::n_dofs() const
+ {
+ return Solver<dim>::n_dofs();
+ }
+
+
+ template <int dim>
+ void
+ DualSolver<dim>::
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
+ {
+ Solver<dim>::postprocess(postprocessor);
+ }
+
+
+
+ template <int dim>
+ void
+ DualSolver<dim>::
+ assemble_rhs (Vector<double> &rhs) const
+ {
+ dual_functional->assemble_rhs (this->dof_handler, rhs);
+ }
+
+
+ // @sect4{The WeightedResidual class}
+
+ // Here finally comes the main
+ // class of this program, the one
+ // that implements the dual
+ // weighted residual error
+ // estimator. It joins the primal
+ // and dual solver classes to use
+ // them for the computation of
+ // primal and dual solutions, and
+ // implements the error
+ // representation formula for use
+ // as error estimate and mesh
+ // refinement.
+ //
+ // The first few of the functions
+ // of this class are mostly
+ // overriders of the respective
+ // functions of the base class:
+ template <int dim>
+ class WeightedResidual : public PrimalSolver<dim>,
+ public DualSolver<dim>
+ {
+ public:
+ WeightedResidual (Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &primal_fe,
+ const FiniteElement<dim> &dual_fe,
+ const Quadrature<dim> &quadrature,
+ const Quadrature<dim-1> &face_quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &boundary_values,
+ const DualFunctional::DualFunctionalBase<dim> &dual_functional);
+
+ virtual
+ void
+ solve_problem ();
+
+ virtual
+ void
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
+
+ virtual
+ unsigned int
+ n_dofs () const;
+
+ virtual void refine_grid ();
+
+ virtual
+ void
+ output_solution () const;
+
+
+ // In the private section, we
+ // have two functions that are
+ // used to call the
+ // ``solve_problem'' functions
+ // of the primal and dual base
+ // classes. These two functions
+ // will be called in parallel
+ // by the ``solve_problem''
+ // function of this class.
+ void solve_primal_problem ();
+ void solve_dual_problem ();
+ // Then declare abbreviations
+ // for active cell iterators,
+ // to avoid that we have to
+ // write this lengthy name
+ // over and over again:
+
+ typedef
+ typename DoFHandler<dim>::active_cell_iterator
+ active_cell_iterator;
+
+ // Next, declare a data type
+ // that we will us to store the
+ // contribution of faces to the
+ // error estimator. The idea is
+ // that we can compute the face
+ // terms from each of the two
+ // cells to this face, as they
+ // are the same when viewed
+ // from both sides. What we
+ // will do is to compute them
+ // only once, based on some
+ // rules explained below which
+ // of the two adjacent cells
+ // will be in charge to do
+ // so. We then store the
+ // contribution of each face in
+ // a map mapping faces to their
+ // values, and only collect the
+ // contributions for each cell
+ // by looping over the cells a
+ // second time and grabbing the
+ // values from the map.
+ //
+ // The data type of this map is
+ // declared here:
+ typedef
+ typename std::map<typename DoFHandler<dim>::face_iterator,double>
+ FaceIntegrals;
+
+ // In the computation of the
+ // error estimates on cells and
+ // faces, we need a number of
+ // helper objects, such as
+ // ``FEValues'' and
+ // ``FEFaceValues'' functions,
+ // but also temporary objects
+ // storing the values and
+ // gradients of primal and dual
+ // solutions, for
+ // example. These fields are
+ // needed in the three
+ // functions that do the
+ // integration on cells, and
+ // regular and irregular faces,
+ // respectively.
+ //
+ // There are three reasonable
+ // ways to provide these
+ // fields: first, as local
+ // variables in the function
+ // that needs them; second, as
+ // member variables of this
+ // class; third, as arguments
+ // passed to that function.
+ //
+ // These three alternatives all
+ // have drawbacks: the third
+ // that their number is not
+ // neglectable and would make
+ // calling these functions a
+ // lengthy enterprise. The
+ // second has the drawback that
+ // it disallows
+ // parallelization, since the
+ // threads that will compute
+ // the error estimate have to
+ // have their own copies of
+ // these variables each, so
+ // member variables of the
+ // enclosing class will not
+ // work. The first approach,
+ // although straightforward,
+ // has a subtle but important
+ // drawback: we will call these
+ // functions over and over
+ // again, many thousand times
+ // maybe; it has now turned out
+ // that allocating vectors and
+ // other objects that need
+ // memory from the heap is an
+ // expensive business in terms
+ // of run-time, since memory
+ // allocation is expensive when
+ // several threads are
+ // involved. In our experience,
+ // more than 20 per cent of the
+ // total run time of error
+ // estimation functions are due
+ // to memory allocation, if
+ // done on a per-call level. It
+ // is thus significantly better
+ // to allocate the memory only
+ // once, and recycle the
+ // objects as often as
+ // possible.
+ //
+ // What to do? Our answer is to
+ // use a variant of the third
+ // strategy, namely generating
+ // these variables once in the
+ // main function of each
+ // thread, and passing them
+ // down to the functions that
+ // do the actual work. To avoid
+ // that we have to give these
+ // functions a dozen or so
+ // arguments, we pack all these
+ // variables into two
+ // structures, one which is
+ // used for the computations on
+ // cells, the other doing them
+ // on the faces. Instead of
+ // many individual objects, we
+ // will then only pass one such
+ // object to these functions,
+ // making their calling
+ // sequence simpler.
+ struct CellData
+ {
+ MappingQ<dim> mapping;
+ FEValues<dim> fe_values;
+ const SmartPointer<const Function<dim> > right_hand_side;
+
+ std::vector<double> cell_residual;
+ std::vector<double> rhs_values;
+ std::vector<double> dual_weights;
+ typename std::vector<Tensor<2,dim> > cell_grad_grads;
+ CellData (const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &right_hand_side);
+ };
+
+ struct FaceData
+ {
+ MappingQ<dim> mapping;
+ FEFaceValues<dim> fe_face_values_cell;
+ FEFaceValues<dim> fe_face_values_neighbor;
+ FESubfaceValues<dim> fe_subface_values_cell;
+
+ std::vector<double> jump_residual;
+ std::vector<double> dual_weights;
+ typename std::vector<Tensor<1,dim> > cell_grads;
+ typename std::vector<Tensor<1,dim> > neighbor_grads;
+ FaceData (const FiniteElement<dim> &fe,
+ const Quadrature<dim-1> &face_quadrature);
+ };
+
+
+
+ // Regarding the evaluation of
+ // the error estimator, we have
+ // two driver functions that do
+ // this: the first is called to
+ // generate the cell-wise
+ // estimates, and splits up the
+ // task in a number of threads
+ // each of which work on a
+ // subset of the cells. The
+ // first function will run the
+ // second for each of these
+ // threads:
+ void estimate_error (Vector<float> &error_indicators) const;
+
+ void estimate_some (const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ const unsigned int n_threads,
+ const unsigned int this_thread,
+ Vector<float> &error_indicators,
+ FaceIntegrals &face_integrals) const;
+
+ // Then we have functions that
+ // do the actual integration of
+ // the error representation
+ // formula. They will treat the
+ // terms on the cell interiors,
+ // on those faces that have no
+ // hanging nodes, and on those
+ // faces with hanging nodes,
+ // respectively:
+ void
+ integrate_over_cell (const active_cell_iterator &cell,
+ const unsigned int cell_index,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ CellData &cell_data,
+ Vector<float> &error_indicators) const;
+
+ void
+ integrate_over_regular_face (const active_cell_iterator &cell,
+ const unsigned int face_no,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ FaceData &face_data,
+ FaceIntegrals &face_integrals) const;
+ void
+ integrate_over_irregular_face (const active_cell_iterator &cell,
+ const unsigned int face_no,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ FaceData &face_data,
+ FaceIntegrals &face_integrals) const;
+ };
+
+
+
+ // In the implementation of this
+ // class, we first have the
+ // constructors of the ``CellData''
+ // and ``FaceData'' member classes,
+ // and the ``WeightedResidual''
+ // constructor. They only
+ // initialize fields to their
+ // correct lengths, so we do not
+ // have to discuss them to length.
+ template <int dim>
+ WeightedResidual<dim>::CellData::
+ CellData (const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &right_hand_side)
+ :
+ mapping (4),
+ fe_values (mapping, fe, quadrature,
+ update_values |
+ update_second_derivatives |
+ update_q_points |
+ update_JxW_values),
+ right_hand_side (&right_hand_side)
+ {
+ const unsigned int n_q_points
+ = quadrature.n_quadrature_points;
+
+ cell_residual.resize(n_q_points);
+ rhs_values.resize(n_q_points);
+ dual_weights.resize(n_q_points);
+ cell_grad_grads.resize(n_q_points);
+ }
+
+
+
+ template <int dim>
+ WeightedResidual<dim>::FaceData::
+ FaceData (const FiniteElement<dim> &fe,
+ const Quadrature<dim-1> &face_quadrature)
+ :
+ mapping (4),
+ fe_face_values_cell (mapping, fe, face_quadrature,
+ update_values |
+ update_gradients |
+ update_JxW_values |
+ update_normal_vectors),
+ fe_face_values_neighbor (mapping, fe, face_quadrature,
+ update_values |
+ update_gradients |
+ update_JxW_values |
+ update_normal_vectors),
+ fe_subface_values_cell (mapping, fe, face_quadrature,
+ update_gradients)
+ {
+ const unsigned int n_face_q_points
+ = face_quadrature.n_quadrature_points;
+
+ jump_residual.resize(n_face_q_points);
+ dual_weights.resize(n_face_q_points);
+ cell_grads.resize(n_face_q_points);
+ neighbor_grads.resize(n_face_q_points);
+ }
+
+
+
+
+ template <int dim>
+ WeightedResidual<dim>::
+ WeightedResidual (Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &primal_fe,
+ const FiniteElement<dim> &dual_fe,
+ const Quadrature<dim> &quadrature,
+ const Quadrature<dim-1> &face_quadrature,
+ const Function<dim> &rhs_function,
+ const Function<dim> &bv,
+ const DualFunctional::DualFunctionalBase<dim> &dual_functional)
+ :
+ Base<dim> (coarse_grid),
+ PrimalSolver<dim> (coarse_grid, primal_fe,
+ quadrature, face_quadrature,
+ rhs_function, bv),
+ DualSolver<dim> (coarse_grid, dual_fe,
+ quadrature, face_quadrature,
+ dual_functional)
+ {}
+
+
+ // The next five functions are
+ // boring, as they simply relay
+ // their work to the base
+ // classes. The first calls the
+ // primal and dual solvers in
+ // parallel, while postprocessing
+ // the solution and retrieving the
+ // number of degrees of freedom is
+ // done by the primal class.
+ template <int dim>
+ void
+ WeightedResidual<dim>::solve_problem ()
+ {
+ Threads::ThreadGroup<> threads;
+ threads += Threads::spawn (*this, &WeightedResidual<dim>::solve_primal_problem)();
+ threads += Threads::spawn (*this, &WeightedResidual<dim>::solve_dual_problem)();
+ threads.join_all ();
+ }
+
+
+ template <int dim>
+ void
+ WeightedResidual<dim>::solve_primal_problem ()
+ {
+ PrimalSolver<dim>::solve_problem ();
+ }
+
+ template <int dim>
+ void
+ WeightedResidual<dim>::solve_dual_problem ()
+ {
+ DualSolver<dim>::solve_problem ();
+ }
+
+
+ template <int dim>
+ void
+ WeightedResidual<dim>::
+ postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
+ {
+ PrimalSolver<dim>::postprocess (postprocessor);
+ }
+
+
+ template <int dim>
+ unsigned int
+ WeightedResidual<dim>::n_dofs () const
+ {
+ return PrimalSolver<dim>::n_dofs();
+ }
+
+
+
+ // Now, it is becoming more
+ // interesting: the ``refine_grid''
+ // function asks the error
+ // estimator to compute the
+ // cell-wise error indicators, then
+ // uses their absolute values for
+ // mesh refinement.
+ template <int dim>
+ void
+ WeightedResidual<dim>::refine_grid ()
+ {
+ // First call the function that
+ // computes the cell-wise and
+ // global error:
+ Vector<float> error_indicators (this->triangulation->n_active_cells());
+ estimate_error (error_indicators);
+
+ // Then note that marking cells
+ // for refinement or coarsening
+ // only works if all indicators
+ // are positive, to allow their
+ // comparison. Thus, drop the
+ // signs on all these indicators:
+ for (Vector<float>::iterator i=error_indicators.begin();
+ i != error_indicators.end(); ++i)
+ *i = std::fabs (*i);
+
+ // Finally, we can select between
+ // different strategies for
+ // refinement. The default here
+ // is to refine those cells with
+ // the largest error indicators
+ // that make up for a total of 80
+ // per cent of the error, while
+ // we coarsen those with the
+ // smallest indicators that make
+ // up for the bottom 2 per cent
+ // of the error.
+ GridRefinement::refine_and_coarsen_fixed_fraction (*this->triangulation,
+ error_indicators,
+ 0.5/std::sqrt(1.+step),
+ 0.2/std::sqrt(1.+step));
+ this->triangulation->execute_coarsening_and_refinement ();
+ }
+
+
+ // Since we want to output both the
+ // primal and the dual solution, we
+ // overload the ``output_solution''
+ // function. The only interesting
+ // feature of this function is that
+ // the primal and dual solutions
+ // are defined on different finite
+ // element spaces, which is not the
+ // format the ``DataOut'' class
+ // expects. Thus, we have to
+ // transfer them to a common finite
+ // element space. Since we want the
+ // solutions only to see them
+ // qualitatively, we contend
+ // ourselves with interpolating the
+ // dual solution to the (smaller)
+ // primal space. For the
+ // interpolation, there is a
+ // library function, that takes a
+ // ``ConstraintMatrix'' object
+ // including the hanging node
+ // constraints. The rest is
+ // standard.
+ //
+ // There is, however, one
+ // work-around worth mentioning: in
+ // this function, as in a couple of
+ // following ones, we have to
+ // access the ``DoFHandler''
+ // objects and solutions of both
+ // the primal as well as of the
+ // dual solver. Since these are
+ // members of the ``Solver'' base
+ // class which exists twice in the
+ // class hierarchy leading to the
+ // present class (once as base
+ // class of the ``PrimalSolver''
+ // class, once as base class of the
+ // ``DualSolver'' class), we have
+ // to disambiguate accesses to them
+ // by telling the compiler a member
+ // of which of these two instances
+ // we want to access. The way to do
+ // this would be identify the
+ // member by pointing a path
+ // through the class hierarchy
+ // which disambiguates the base
+ // class, for example writing
+ // ``PrimalSolver::dof_handler'' to
+ // denote the member variable
+ // ``dof_handler'' from the
+ // ``Solver'' base class of the
+ // ``PrimalSolver''
+ // class. Unfortunately, this
+ // confuses gcc's version 2.96 (a
+ // version that was intended as a
+ // development snapshot, but
+ // delivered as system compiler by
+ // Red Hat in their 7.x releases)
+ // so much that it bails out and
+ // refuses to compile the code.
+ //
+ // Thus, we have to work around
+ // this problem. We do this by
+ // introducing references to the
+ // ``PrimalSolver'' and
+ // ``DualSolver'' components of the
+ // ``WeightedResidual'' object at
+ // the beginning of the
+ // function. Since each of these
+ // has an unambiguous base class
+ // ``Solver'', we can access the
+ // member variables we want through
+ // these references. However, we
+ // are now accessing protected
+ // member variables of these
+ // classes through a pointer other
+ // than the ``this'' pointer (in
+ // fact, this is of course the
+ // ``this'' pointer, but not
+ // explicitly). This finally is the
+ // reason why we had to declare the
+ // present class a friend of the
+ // classes we so access.
+ template <int dim>
+ void
+ WeightedResidual<dim>::output_solution () const
+ {
+ const PrimalSolver<dim> &primal_solver = *this;
+ const DualSolver<dim> &dual_solver = *this;
+
+ ConstraintMatrix primal_hanging_node_constraints;
+ DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
+ primal_hanging_node_constraints);
+ primal_hanging_node_constraints.close();
+ Vector<double> dual_solution (primal_solver.dof_handler.n_dofs());
+ FETools::interpolate (dual_solver.dof_handler,
+ dual_solver.solution,
+ primal_solver.dof_handler,
+ primal_hanging_node_constraints,
+ dual_solution);
+
+ // approximate error, gradient,
+ // and second derivative
+ // information as cell information
+ Vector<float> error_indicators (this->triangulation->n_active_cells());
+ Vector<float> gradient_indicators (this->triangulation->n_active_cells());
+ Vector<float> second_indicators (this->triangulation->n_active_cells());
+ {
+ MappingQ<dim> mapping(4);
+ KellyErrorEstimator<dim>::estimate (mapping, primal_solver.dof_handler,
+ QGauss3<dim-1>(),
+ typename FunctionMap<dim>::type(),
+ primal_solver.solution,
+ error_indicators);
+
+ DerivativeApproximation::
+ approximate_gradient (mapping,
+ primal_solver.dof_handler,
+ primal_solver.solution,
+ gradient_indicators);
+
+ DerivativeApproximation::
+ approximate_second_derivative (mapping,
+ primal_solver.dof_handler,
+ primal_solver.solution,
+ second_indicators);
+
+ }
+ // distribute cell to dof vectors
+ Vector<double> x_error_indicators (primal_solver.dof_handler.n_dofs());
+ Vector<double> x_gradient_indicators (primal_solver.dof_handler.n_dofs());
+ Vector<double> x_second_indicators (primal_solver.dof_handler.n_dofs());
+ DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
+ error_indicators,
+ x_error_indicators);
+ DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
+ gradient_indicators,
+ x_gradient_indicators);
+ DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler,
+ second_indicators,
+ x_second_indicators);
+
+
+
+ // we generate too much output in
+ // 3d. instead of doing it that
+ // way, simply generate a coarser
+ // mesh and output from there
+ Triangulation<dim> coarser_mesh;
+ coarser_mesh.copy_triangulation (*this->triangulation);
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = coarser_mesh.begin_active();
+ cell != coarser_mesh.end(); ++cell)
+ cell->set_coarsen_flag();
+ coarser_mesh.execute_coarsening_and_refinement ();
+
+ // next generate a DoF handler on
+ // that mesh and a map fron one
+ // to the other mesh
+ DoFHandler<dim> coarser_dof_handler (coarser_mesh);
+ coarser_dof_handler.distribute_dofs (primal_solver.dof_handler.get_fe());
+ InterGridMap<DoFHandler<dim> > coarse_to_fine_map;
+ coarse_to_fine_map.make_mapping (coarser_dof_handler,
+ primal_solver.dof_handler);
+
+ // finally we have to transfer
+ // the data vectors
+ Vector<double> coarse_primal_solution (coarser_dof_handler.n_dofs());
+ Vector<double> coarse_dual_solution (coarser_dof_handler.n_dofs());
+ Vector<double> coarse_error_indicators (coarser_dof_handler.n_dofs());
+ Vector<double> coarse_gradient_indicators (coarser_dof_handler.n_dofs());
+ Vector<double> coarse_second_indicators (coarser_dof_handler.n_dofs());
+
+ Vector<double> tmp (coarser_dof_handler.get_fe().dofs_per_cell);
+ for (typename DoFHandler<dim>::active_cell_iterator
+ cell = coarser_dof_handler.begin_active();
+ cell != coarser_dof_handler.end(); ++cell)
+ {
+ coarse_to_fine_map[cell]->get_interpolated_dof_values (primal_solver.solution,tmp);
+ cell->set_dof_values (tmp, coarse_primal_solution);
+
+ coarse_to_fine_map[cell]->get_interpolated_dof_values (dual_solution,tmp);
+ cell->set_dof_values (tmp, coarse_dual_solution);
+
+ coarse_to_fine_map[cell]->get_interpolated_dof_values (x_error_indicators,tmp);
+ cell->set_dof_values (tmp, coarse_error_indicators);
+
+ coarse_to_fine_map[cell]->get_interpolated_dof_values (x_gradient_indicators,tmp);
+ cell->set_dof_values (tmp, coarse_gradient_indicators);
+
+ coarse_to_fine_map[cell]->get_interpolated_dof_values (x_second_indicators,tmp);
+ cell->set_dof_values (tmp, coarse_second_indicators);
+ }
+
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (coarser_dof_handler);
+ data_out.add_data_vector (coarse_primal_solution, "primal_solution");
+ data_out.add_data_vector (coarse_dual_solution, "dual_solution");
+ data_out.add_data_vector (coarse_error_indicators, "errors");
+ data_out.add_data_vector (coarse_gradient_indicators, "gradient");
+ data_out.add_data_vector (coarse_second_indicators, "second_derivatives");
+ data_out.build_patches ();
+
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ostringstream filename;
+#else
+ std::ostrstream filename;
+#endif
+ filename << "spec2006-447.dealII/"
+ << "solution-"
+ << this->refinement_cycle
+ << ".gmv"
+ << std::ends;
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ofstream out (filename.str().c_str());
+#else
+ std::ofstream out (filename.str());
+#endif
+
+ data_out.write_gmv (out);
+ }
+
+ }
+
+
+ // @sect3{Estimating errors}
+
+ // @sect4{Error estimation driver functions}
+ //
+ // As for the actual computation of
+ // error estimates, let's start
+ // with the function that drives
+ // all this, i.e. calls those
+ // functions that actually do the
+ // work, and finally collects the
+ // results.
+
+ template <int dim>
+ void
+ WeightedResidual<dim>::
+ estimate_error (Vector<float> &error_indicators) const
+ {
+ const PrimalSolver<dim> &primal_solver = *this;
+ const DualSolver<dim> &dual_solver = *this;
+
+ // The first task in computing
+ // the error is to set up vectors
+ // that denote the primal
+ // solution, and the weights
+ // (z-z_h)=(z-I_hz), both in the
+ // finite element space for which
+ // we have computed the dual
+ // solution. For this, we have to
+ // interpolate the primal
+ // solution to the dual finite
+ // element space, and to subtract
+ // the interpolation of the
+ // computed dual solution to the
+ // primal finite element
+ // space. Fortunately, the
+ // library provides functions for
+ // the interpolation into larger
+ // or smaller finite element
+ // spaces, so this is mostly
+ // obvious.
+ //
+ // First, let's do that for the
+ // primal solution: it is
+ // cell-wise interpolated into
+ // the finite element space in
+ // which we have solved the dual
+ // problem: But, again as in the
+ // ``WeightedResidual::output_solution''
+ // function we first need to
+ // create a ConstraintMatrix
+ // including the hanging node
+ // constraints, but this time of
+ // the dual finite element space.
+ ConstraintMatrix dual_hanging_node_constraints;
+ DoFTools::make_hanging_node_constraints (dual_solver.dof_handler,
+ dual_hanging_node_constraints);
+ dual_hanging_node_constraints.close();
+ Vector<double> primal_solution (dual_solver.dof_handler.n_dofs());
+ FETools::interpolate (primal_solver.dof_handler,
+ primal_solver.solution,
+ dual_solver.dof_handler,
+ dual_hanging_node_constraints,
+ primal_solution);
+
+ // Then for computing the
+ // interpolation of the
+ // numerically approximated dual
+ // solution z into the finite
+ // element space of the primal
+ // solution and subtracting it
+ // from z: use the
+ // ``interpolate_difference''
+ // function, that gives (z-I_hz)
+ // in the element space of the
+ // dual solution.
+ ConstraintMatrix primal_hanging_node_constraints;
+ DoFTools::make_hanging_node_constraints (primal_solver.dof_handler,
+ primal_hanging_node_constraints);
+ primal_hanging_node_constraints.close();
+ Vector<double> dual_weights (dual_solver.dof_handler.n_dofs());
+ FETools::interpolation_difference (dual_solver.dof_handler,
+ dual_hanging_node_constraints,
+ dual_solver.solution,
+ primal_solver.dof_handler,
+ primal_hanging_node_constraints,
+ dual_weights);
+
+ // Note that this could probably
+ // have been more efficient since
+ // those constraints have been
+ // used previously when
+ // assembling matrix and right
+ // hand side for the primal
+ // problem and writing out the
+ // dual solution. We leave the
+ // optimization of the program in
+ // this respect as an exercise.
+
+ // Having computed the dual
+ // weights we now proceed with
+ // computing the cell and face
+ // residuals of the primal
+ // solution. First we set up a
+ // map between face iterators and
+ // their jump term contributions
+ // of faces to the error
+ // estimator. The reason is that
+ // we compute the jump terms only
+ // once, from one side of the
+ // face, and want to collect them
+ // only afterwards when looping
+ // over all cells a second time.
+ //
+ // We initialize this map already
+ // with a value of -1e20 for all
+ // faces, since this value will
+ // strike in the results if
+ // something should go wrong and
+ // we fail to compute the value
+ // for a face for some
+ // reason. Secondly, we
+ // initialize the map once before
+ // we branch to different threads
+ // since this way the map's
+ // structure is no more modified
+ // by the individual threads,
+ // only existing entries are set
+ // to new values. This relieves
+ // us from the necessity to
+ // synchronise the threads
+ // through a mutex each time they
+ // write to (and modify the
+ // structure of) this map.
+ FaceIntegrals face_integrals;
+ for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
+ cell!=dual_solver.dof_handler.end();
+ ++cell)
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ face_integrals[cell->face(face_no)] = -1e20;
+
+ // Then set up a vector with
+ // error indicators. Reserve one
+ // slot for each cell and set it
+ // to zero.
+ error_indicators.reinit (dual_solver.dof_handler
+ .get_tria().n_active_cells());
+
+ // Now start a number of threads
+ // which compute the error
+ // formula on parts of all the
+ // cells, and once they are all
+ // started wait until they have
+ // all finished:
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+ for (unsigned int i=0; i<n_threads; ++i)
+ threads += Threads::spawn (*this, &WeightedResidual<dim>::estimate_some)
+ (primal_solution,
+ dual_weights,
+ n_threads, i,
+ error_indicators,
+ face_integrals);
+ threads.join_all();
+
+ // Once the error contributions
+ // are computed, sum them up. For
+ // this, note that the cell terms
+ // are already set, and that only
+ // the edge terms need to be
+ // collected. Thus, loop over all
+ // cells and their faces, make
+ // sure that the contributions of
+ // each of the faces are there,
+ // and add them up. Only take
+ // minus one half of the jump
+ // term, since the other half
+ // will be taken by the
+ // neighboring cell.
+ unsigned int present_cell=0;
+ for (active_cell_iterator cell=dual_solver.dof_handler.begin_active();
+ cell!=dual_solver.dof_handler.end();
+ ++cell, ++present_cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ Assert(face_integrals.find(cell->face(face_no)) !=
+ face_integrals.end(),
+ ExcInternalError());
+ error_indicators(present_cell)
+ -= 0.5*face_integrals[cell->face(face_no)];
+ };
+ std::cout << " Estimated error="
+ << std::accumulate (error_indicators.begin(),
+ error_indicators.end(), 0.)
+ << std::endl;
+ }
+
+
+ // @sect4{Estimating on a subset of cells}
+
+ // Next we have the function that
+ // is called to estimate the error
+ // on a subset of cells. The
+ // function may be called multiply
+ // if the library was configured to
+ // use multi-threading. Here it
+ // goes:
+ template <int dim>
+ void
+ WeightedResidual<dim>::
+ estimate_some (const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ const unsigned int n_threads,
+ const unsigned int this_thread,
+ Vector<float> &error_indicators,
+ FaceIntegrals &face_integrals) const
+ {
+ const PrimalSolver<dim> &primal_solver = *this;
+ const DualSolver<dim> &dual_solver = *this;
+
+ // At the beginning, we
+ // initialize two variables for
+ // each thread which may be
+ // running this function. The
+ // reason for these functions was
+ // discussed above, when the
+ // respective classes were
+ // discussed, so we here only
+ // point out that since they are
+ // local to the function that is
+ // spawned when running more than
+ // one thread, the data of these
+ // objects exists actually once
+ // per thread, so we don't have
+ // to take care about
+ // synchronising access to them.
+ CellData cell_data (*dual_solver.fe,
+ *dual_solver.quadrature,
+ *primal_solver.rhs_function);
+ FaceData face_data (*dual_solver.fe,
+ *dual_solver.face_quadrature);
+
+ // Then calculate the start cell
+ // for this thread. We let the
+ // different threads run on
+ // interleaved cells, i.e. for
+ // example if we have 4 threads,
+ // then the first thread treates
+ // cells 0, 4, 8, etc, while the
+ // second threads works on cells 1,
+ // 5, 9, and so on. The reason is
+ // that it takes vastly more time
+ // to work on cells with hanging
+ // nodes than on regular cells, but
+ // such cells are not evenly
+ // distributed across the range of
+ // cell iterators, so in order to
+ // have the different threads do
+ // approximately the same amount of
+ // work, we have to let them work
+ // interleaved to the effect of a
+ // pseudorandom distribution of the
+ // `hard' cells to the different
+ // threads.
+ active_cell_iterator cell=dual_solver.dof_handler.begin_active();
+ for (unsigned int t=0;
+ (t<this_thread) && (cell!=dual_solver.dof_handler.end());
+ ++t, ++cell);
+
+ // If there are no cells for this
+ // thread (for example if there
+ // are a total of less cells than
+ // there are threads), then go
+ // back right now
+ if (cell == dual_solver.dof_handler.end())
+ return;
+
+ // Next loop over all cells. The
+ // check for loop end is done at
+ // the end of the loop, along
+ // with incrementing the loop
+ // index.
+ for (unsigned int cell_index=this_thread; true; )
+ {
+ // First task on each cell is
+ // to compute the cell
+ // residual contributions of
+ // this cell, and put them
+ // into the
+ // ``error_indicators''
+ // variable:
+ integrate_over_cell (cell, cell_index,
+ primal_solution,
+ dual_weights,
+ cell_data,
+ error_indicators);
+
+ // After computing the cell
+ // terms, turn to the face
+ // terms. For this, loop over
+ // all faces of the present
+ // cell, and see whether
+ // something needs to be
+ // computed on it:
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ // First, if this face is
+ // part of the boundary,
+ // then there is nothing
+ // to do. However, to
+ // make things easier
+ // when summing up the
+ // contributions of the
+ // faces of cells, we
+ // enter this face into
+ // the list of faces with
+ // a zero contribution to
+ // the error.
+ if (cell->face(face_no)->at_boundary())
+ {
+ face_integrals[cell->face(face_no)] = 0;
+ continue;
+ };
+
+ // Next, note that since
+ // we want to compute the
+ // jump terms on each
+ // face only once
+ // although we access it
+ // twice (if it is not at
+ // the boundary), we have
+ // to define some rules
+ // who is responsible for
+ // computing on a face:
+ //
+ // First, if the
+ // neighboring cell is on
+ // the same level as this
+ // one, i.e. neither
+ // further refined not
+ // coarser, then the one
+ // with the lower index
+ // within this level does
+ // the work. In other
+ // words: if the other
+ // one has a lower index,
+ // then skip work on this
+ // face:
+ if ((cell->neighbor(face_no)->has_children() == false) &&
+ (cell->neighbor(face_no)->level() == cell->level()) &&
+ (cell->neighbor(face_no)->index() < cell->index()))
+ continue;
+
+ // Likewise, we always
+ // work from the coarser
+ // cell if this and its
+ // neighbor differ in
+ // refinement. Thus, if
+ // the neighboring cell
+ // is less refined than
+ // the present one, then
+ // do nothing since we
+ // integrate over the
+ // subfaces when we visit
+ // the coarse cell.
+ if (cell->at_boundary(face_no) == false)
+ if (cell->neighbor(face_no)->level() < cell->level())
+ continue;
+
+
+ // Now we know that we
+ // are in charge here, so
+ // actually compute the
+ // face jump terms. If
+ // the face is a regular
+ // one, i.e. the other
+ // side's cell is neither
+ // coarser not finer than
+ // this cell, then call
+ // one function, and if
+ // the cell on the other
+ // side is further
+ // refined, then use
+ // another function. Note
+ // that the case that the
+ // cell on the other side
+ // is coarser cannot
+ // happen since we have
+ // decided above that we
+ // handle this case when
+ // we pass over that
+ // other cell.
+ if (cell->face(face_no)->has_children() == false)
+ integrate_over_regular_face (cell, face_no,
+ primal_solution,
+ dual_weights,
+ face_data,
+ face_integrals);
+ else
+ integrate_over_irregular_face (cell, face_no,
+ primal_solution,
+ dual_weights,
+ face_data,
+ face_integrals);
+ };
+
+ // After computing the cell
+ // contributions and looping
+ // over the faces, go to the
+ // next cell for this
+ // thread. Note again that
+ // the cells for each of the
+ // threads are interleaved.
+ // If we are at the end of
+ // our workload, jump out
+ // of the loop.
+ for (unsigned int t=0;
+ ((t<n_threads) && (cell!=dual_solver.dof_handler.end()));
+ ++t, ++cell, ++cell_index);
+ if (cell == dual_solver.dof_handler.end())
+ break;
+ };
+ }
+
+
+ // @sect4{Computing cell term error contributions}
+
+ // As for the actual computation of
+ // the error contributions, first
+ // turn to the cell terms:
+ template <int dim>
+ void WeightedResidual<dim>::
+ integrate_over_cell (const active_cell_iterator &cell,
+ const unsigned int cell_index,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ CellData &cell_data,
+ Vector<float> &error_indicators) const
+ {
+ // The tasks to be done are what
+ // appears natural from looking
+ // at the error estimation
+ // formula: first compute the the
+ // right hand side and the
+ // Laplacian of the numerical
+ // solution at the quadrature
+ // points for the cell residual,
+ cell_data.fe_values.reinit (cell);
+ cell_data.right_hand_side
+ ->value_list (cell_data.fe_values.get_quadrature_points(),
+ cell_data.rhs_values);
+ cell_data.fe_values.get_function_2nd_derivatives (primal_solution,
+ cell_data.cell_grad_grads);
+
+ // ...then get the dual weights...
+ cell_data.fe_values.get_function_values (dual_weights,
+ cell_data.dual_weights);
+
+ // ...and finally build the sum
+ // over all quadrature points and
+ // store it with the present
+ // cell:
+ double sum = 0;
+ for (unsigned int p=0; p<cell_data.fe_values.n_quadrature_points; ++p)
+ sum += ((cell_data.rhs_values[p]+trace(cell_data.cell_grad_grads[p])) *
+ cell_data.dual_weights[p] *
+ cell_data.fe_values.JxW (p));
+ error_indicators(cell_index) += sum;
+ }
+
+
+ // @sect4{Computing edge term error contributions - 1}
+
+ // On the other hand, computation
+ // of the edge terms for the error
+ // estimate is not so
+ // simple. First, we have to
+ // distinguish between faces with
+ // and without hanging
+ // nodes. Because it is the simple
+ // case, we first consider the case
+ // without hanging nodes on a face
+ // (let's call this the `regular'
+ // case):
+ template <int dim>
+ void WeightedResidual<dim>::
+ integrate_over_regular_face (const active_cell_iterator &cell,
+ const unsigned int face_no,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ FaceData &face_data,
+ FaceIntegrals &face_integrals) const
+ {
+ const unsigned int
+ n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+
+ // The first step is to get the
+ // values of the gradients at the
+ // quadrature points of the
+ // finite element field on the
+ // present cell. For this,
+ // initialize the
+ // ``FEFaceValues'' object
+ // corresponding to this side of
+ // the face, and extract the
+ // gradients using that
+ // object.
+ face_data.fe_face_values_cell.reinit (cell, face_no);
+ face_data.fe_face_values_cell.get_function_grads (primal_solution,
+ face_data.cell_grads);
+
+ // The second step is then to
+ // extract the gradients of the
+ // finite element solution at the
+ // quadrature points on the other
+ // side of the face, i.e. from
+ // the neighboring cell.
+ //
+ // For this, do a sanity check
+ // before: make sure that the
+ // neigbor actually exists (yes,
+ // we should not have come here
+ // if the neighbor did not exist,
+ // but in complicated software
+ // there are bugs, so better
+ // check this), and if this is
+ // not the case throw an error.
+ Assert (cell->neighbor(face_no).state() == IteratorState::valid,
+ ExcInternalError());
+ // If we have that, then we need
+ // to find out with which face of
+ // the neighboring cell we have
+ // to work, i.e. the
+ // ``home-many''the neighbor the
+ // present cell is of the cell
+ // behind the present face. For
+ // this, there is a function, and
+ // we put the result into a
+ // variable with the name
+ // ``neighbor_neighbor'':
+ const unsigned int
+ neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
+ // Then define an abbreviation
+ // for the neigbor cell,
+ // initialize the
+ // ``FEFaceValues'' object on
+ // that cell, and extract the
+ // gradients on that cell:
+ const active_cell_iterator neighbor = cell->neighbor(face_no);
+ face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
+ face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+ face_data.neighbor_grads);
+
+ // Now that we have the gradients
+ // on this and the neighboring
+ // cell, compute the jump
+ // residual by multiplying the
+ // jump in the gradient with the
+ // normal vector:
+ for (unsigned int p=0; p<n_q_points; ++p)
+ face_data.jump_residual[p]
+ = ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
+ face_data.fe_face_values_cell.normal_vector(p));
+
+ // Next get the dual weights for
+ // this face:
+ face_data.fe_face_values_cell.get_function_values (dual_weights,
+ face_data.dual_weights);
+
+ // Finally, we have to compute
+ // the sum over jump residuals,
+ // dual weights, and quadrature
+ // weights, to get the result for
+ // this face:
+ double face_integral = 0;
+ for (unsigned int p=0; p<n_q_points; ++p)
+ face_integral += (face_data.jump_residual[p] *
+ face_data.dual_weights[p] *
+ face_data.fe_face_values_cell.JxW(p));
+
+ // Double check that the element
+ // already exists and that it was
+ // not already written to...
+ Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
+ ExcInternalError());
+ Assert (face_integrals[cell->face(face_no)] == -1e20,
+ ExcInternalError());
+
+ // ...then store computed value
+ // at assigned location. Note
+ // that the stored value does not
+ // contain the factor 1/2 that
+ // appears in the error
+ // representation. The reason is
+ // that the term actually does
+ // not have this factor if we
+ // loop over all faces in the
+ // triangulation, but only
+ // appears if we write it as a
+ // sum over all cells and all
+ // faces of each cell; we thus
+ // visit the same face twice. We
+ // take account of this by using
+ // this factor -1/2 later, when we
+ // sum up the contributions for
+ // each cell individually.
+ face_integrals[cell->face(face_no)] = face_integral;
+ }
+
+
+ // @sect4{Computing edge term error contributions - 2}
+
+ // We are still missing the case of
+ // faces with hanging nodes. This
+ // is what is covered in this
+ // function:
+ template <int dim>
+ void WeightedResidual<dim>::
+ integrate_over_irregular_face (const active_cell_iterator &cell,
+ const unsigned int face_no,
+ const Vector<double> &primal_solution,
+ const Vector<double> &dual_weights,
+ FaceData &face_data,
+ FaceIntegrals &face_integrals) const
+ {
+ // First again two abbreviations,
+ // and some consistency checks
+ // whether the function is called
+ // only on faces for which it is
+ // supposed to be called:
+ const unsigned int
+ n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+
+ const typename DoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(face_no);
+ Assert (neighbor.state() == IteratorState::valid,
+ ExcInternalError());
+ Assert (neighbor->has_children(),
+ ExcInternalError());
+
+ // Then find out which neighbor
+ // the present cell is of the
+ // adjacent cell. Note that we
+ // will operator on the children
+ // of this adjacent cell, but
+ // that their orientation is the
+ // same as that of their mother,
+ // i.e. the neigbor direction is
+ // the same.
+ const unsigned int
+ neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
+
+ // Then simply do everything we
+ // did in the previous function
+ // for one face for all the
+ // sub-faces now:
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<dim>::subfaces_per_face;
+ ++subface_no)
+ {
+ const active_cell_iterator neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ // Now start the work by
+ // again getting the gradient
+ // of the solution first at
+ // this side of the
+ // interface,
+ face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no);
+ face_data.fe_subface_values_cell.get_function_grads (primal_solution,
+ face_data.cell_grads);
+ // then at the other side,
+ face_data.fe_face_values_neighbor.reinit (neighbor_child,
+ neighbor_neighbor);
+ face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+ face_data.neighbor_grads);
+
+ // and finally building the
+ // jump residuals. Since we
+ // take the normal vector
+ // from the other cell this
+ // time, revert the sign of
+ // the first term compared to
+ // the other function:
+ for (unsigned int p=0; p<n_q_points; ++p)
+ face_data.jump_residual[p]
+ = ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
+ face_data.fe_face_values_neighbor.normal_vector(p));
+
+ // Then get dual weights:
+ face_data.fe_face_values_neighbor.get_function_values (dual_weights,
+ face_data.dual_weights);
+
+ // At last, sum up the
+ // contribution of this
+ // sub-face, and set it in
+ // the global map:
+ double face_integral = 0;
+ for (unsigned int p=0; p<n_q_points; ++p)
+ face_integral += (face_data.jump_residual[p] *
+ face_data.dual_weights[p] *
+ face_data.fe_face_values_neighbor.JxW(p));
+ face_integrals[neighbor_child->face(neighbor_neighbor)]
+ = face_integral;
+ };
+
+ // Once the contributions of all
+ // sub-faces are computed, loop
+ // over all sub-faces to collect
+ // and store them with the mother
+ // face for simple use when later
+ // collecting the error terms of
+ // cells. Again make safety
+ // checks that the entries for
+ // the sub-faces have been
+ // computed and do not carry an
+ // invalid value.
+ double sum = 0;
+ typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<dim>::subfaces_per_face;
+ ++subface_no)
+ {
+ Assert (face_integrals.find(face->child(subface_no)) !=
+ face_integrals.end(),
+ ExcInternalError());
+ Assert (face_integrals[face->child(subface_no)] != -1e20,
+ ExcInternalError());
+
+ sum += face_integrals[face->child(subface_no)];
+ };
+ // Finally store the value with
+ // the parent face.
+ face_integrals[face] = sum;
+ }
+
+}
+
+
+ // @sect3{A simulation framework}
+
+ // In the previous example program,
+ // we have had two functions that
+ // were used to drive the process of
+ // solving on subsequently finer
+ // grids. We extend this here to
+ // allow for a number of parameters
+ // to be passed to these functions,
+ // and put all of that into framework
+ // class.
+ //
+ // You will have noted that this
+ // program is built up of a number of
+ // small parts (evaluation functions,
+ // solver classes implementing
+ // various refinement methods,
+ // different dual functionals,
+ // different problem and data
+ // descriptions), which makes the
+ // program relatively simple to
+ // extend, but also allows to solve a
+ // large number of different problems
+ // by replacing one part by
+ // another. We reflect this
+ // flexibility by declaring a
+ // structure in the following
+ // framework class that holds a
+ // number of parameters that may be
+ // set to test various combinations
+ // of the parts of this program, and
+ // which can be used to test it at
+ // various problems and
+ // discretizations in a simple way.
+template <int dim>
+struct Framework
+{
+ public:
+ // First, we declare two
+ // abbreviations for simple use
+ // of the respective data types:
+ typedef Evaluation::EvaluationBase<dim> Evaluator;
+ typedef std::list<Evaluator*> EvaluatorList;
+
+
+ // Then we have the structure
+ // which declares all the
+ // parameters that may be set. In
+ // the default constructor of the
+ // structure, these values are
+ // all set to default values, for
+ // simple use.
+ struct ProblemDescription
+ {
+ // First allow for the
+ // degrees of the piecewise
+ // polynomials by which the
+ // primal and dual problems
+ // will be discretized. They
+ // default to (bi-,
+ // tri-)linear ansatz
+ // functions for the primal,
+ // and (bi-, tri-)quadratic
+ // ones for the dual
+ // problem. If a refinement
+ // criterion is chosen that
+ // does not need the solution
+ // of a dual problem, the
+ // value of the dual finite
+ // element degree is of
+ // course ignored.
+ unsigned int primal_fe_degree;
+ unsigned int dual_fe_degree;
+
+ // Then have an object that
+ // describes the problem
+ // type, i.e. right hand
+ // side, domain, boundary
+ // values, etc. The pointer
+ // needed here defaults to
+ // the Null pointer, i.e. you
+ // will have to set it in
+ // actual instances of this
+ // object to make it useful.
+ SmartPointer<const Data::SetUpBase<dim> > data;
+
+ // Since we allow to use
+ // different refinement
+ // criteria (global
+ // refinement, refinement by
+ // the Kelly error indicator,
+ // possibly with a weight,
+ // and using the dual
+ // estimator), define a
+ // number of enumeration
+ // values, and subsequently a
+ // variable of that type. It
+ // will default to
+ // ``dual_weighted_error_estimator''.
+ enum RefinementCriterion {
+ dual_weighted_error_estimator,
+ global_refinement,
+ kelly_indicator,
+ weighted_kelly_indicator
+ };
+
+ RefinementCriterion refinement_criterion;
+
+ // Next, an object that
+ // describes the dual
+ // functional. It is only
+ // needed if the dual
+ // weighted residual
+ // refinement is chosen, and
+ // also defaults to a Null
+ // pointer.
+ SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;
+
+ // Then a list of evaluation
+ // objects. Its default value
+ // is empty, i.e. no
+ // evaluation objects.
+ EvaluatorList evaluator_list;
+
+ // Next to last, a function
+ // that is used as a weight
+ // to the
+ // ``RefinementWeightedKelly''
+ // class. The default value
+ // of this pointer is zero,
+ // but you have to set it to
+ // some other value if you
+ // want to use the
+ // ``weighted_kelly_indicator''
+ // refinement criterion.
+ SmartPointer<const Function<dim> > kelly_weight;
+
+ // Finally, we have a
+ // variable that denotes the
+ // maximum number of degrees
+ // of freedom we allow for
+ // the (primal)
+ // discretization. If it is
+ // exceeded, we stop the
+ // process of solving and
+ // intermittend mesh
+ // refinement. Its default
+ // value is 20,000.
+ unsigned int max_degrees_of_freedom;
+
+ // Finally the default
+ // constructor of this class:
+ ProblemDescription ();
+ };
+
+ // The driver framework class
+ // only has one method which
+ // calls solver and mesh
+ // refinement intermittently, and
+ // does some other small tasks in
+ // between. Since it does not
+ // need data besides the
+ // parameters given to it, we
+ // make it static:
+ static void run (const ProblemDescription &descriptor);
+};
+
+
+ // As for the implementation, first
+ // the constructor of the parameter
+ // object, setting all values to
+ // their defaults:
+template <int dim>
+Framework<dim>::ProblemDescription::ProblemDescription ()
+ :
+ primal_fe_degree (1),
+ dual_fe_degree (2),
+ refinement_criterion (dual_weighted_error_estimator),
+ max_degrees_of_freedom (1000)
+{}
+
+
+
+ // Then the function which drives the
+ // whole process:
+template <int dim>
+void Framework<dim>::run (const ProblemDescription &descriptor)
+{
+ // First create a triangulation
+ // from the given data object,
+ Triangulation<dim>
+ triangulation (Triangulation<dim>::smoothing_on_refinement);
+ descriptor.data->create_coarse_grid (triangulation);
+
+ // then a set of finite elements
+ // and appropriate quadrature
+ // formula:
+ const FE_Q<dim> primal_fe(descriptor.primal_fe_degree);
+ const FE_Q<dim> dual_fe(descriptor.dual_fe_degree);
+ const QGauss<dim> quadrature(descriptor.dual_fe_degree+1);
+ const QGauss<dim-1> face_quadrature(descriptor.dual_fe_degree+1);
+
+ // Next, select one of the classes
+ // implementing different
+ // refinement criteria.
+ LaplaceSolver::Base<dim> * solver = 0;
+ switch (descriptor.refinement_criterion)
+ {
+ case ProblemDescription::dual_weighted_error_estimator:
+ {
+ solver
+ = new LaplaceSolver::WeightedResidual<dim> (triangulation,
+ primal_fe,
+ dual_fe,
+ quadrature,
+ face_quadrature,
+ descriptor.data->get_right_hand_side(),
+ descriptor.data->get_boundary_values(),
+ *descriptor.dual_functional);
+ break;
+ };
+
+ case ProblemDescription::global_refinement:
+ {
+ solver
+ = new LaplaceSolver::RefinementGlobal<dim> (triangulation,
+ primal_fe,
+ quadrature,
+ face_quadrature,
+ descriptor.data->get_right_hand_side(),
+ descriptor.data->get_boundary_values());
+ break;
+ };
+
+ case ProblemDescription::kelly_indicator:
+ {
+ solver
+ = new LaplaceSolver::RefinementKelly<dim> (triangulation,
+ primal_fe,
+ quadrature,
+ face_quadrature,
+ descriptor.data->get_right_hand_side(),
+ descriptor.data->get_boundary_values());
+ break;
+ };
+
+ case ProblemDescription::weighted_kelly_indicator:
+ {
+ solver
+ = new LaplaceSolver::RefinementWeightedKelly<dim> (triangulation,
+ primal_fe,
+ quadrature,
+ face_quadrature,
+ descriptor.data->get_right_hand_side(),
+ descriptor.data->get_boundary_values(),
+ *descriptor.kelly_weight);
+ break;
+ };
+
+ default:
+ AssertThrow (false, ExcInternalError());
+ };
+
+ // Now that all objects are in
+ // place, run the main loop. The
+ // stopping criterion is
+ // implemented at the bottom of the
+ // loop.
+ //
+ // In the loop, first set the new
+ // cycle number, then solve the
+ // problem, output its solution(s),
+ // apply the evaluation objects to
+ // it, then decide whether we want
+ // to refine the mesh further and
+ // solve again on this mesh, or
+ // jump out of the loop.
+ for (step=0; step<=n_steps; ++step)
+ {
+ std::cout << "Refinement cycle: " << step
+ << std::endl;
+
+ solver->set_refinement_cycle (step);
+ solver->solve_problem ();
+ solver->output_solution ();
+
+ std::cout << " Number of degrees of freedom="
+ << solver->n_dofs() << std::endl;
+
+ for (typename EvaluatorList::const_iterator
+ e = descriptor.evaluator_list.begin();
+ e != descriptor.evaluator_list.end(); ++e)
+ {
+ (*e)->set_refinement_cycle (step);
+ solver->postprocess (**e);
+ };
+
+
+ if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
+ solver->refine_grid ();
+ else
+ break;
+ };
+
+ // After the loop has run, clean up
+ // the screen, and delete objects
+ // no more needed:
+ std::cout << std::endl;
+ delete solver;
+ solver = 0;
+}
+
+
+
+
+ // @sect3{The main function}
+
+ // Here finally comes the main
+ // function. It drives the whole
+ // process by specifying a set of
+ // parameters to be used for the
+ // simulation (polynomial degrees,
+ // evaluation and dual functionals,
+ // etc), and passes them packed into
+ // a structure to the frame work
+ // class above.
+int main (int argc, char **argv)
+{
+ // if no argument is given, then do 18
+ // iterations
+ if (argc == 1)
+ n_steps = 18;
+ else
+ if (argc == 2)
+ {
+ n_steps = atoi(argv[1]);
+ if ((n_steps==0) || (n_steps>100))
+ {
+ std::cout << "Please call this program with an argument in the range 1..100"
+ << std::endl;
+ exit (1);
+ }
+ }
+ else
+ {
+ std::cout << "Please call this program with a single argument in the range 1..100"
+ << std::endl;
+ exit (1);
+ }
+
+
+
+
+ deallog.depth_console (0);
+ try
+ {
+ // Describe the problem we want
+ // to solve here by passing a
+ // descriptor object to the
+ // function doing the rest of
+ // the work:
+ const unsigned int dim = 3;
+ Framework<dim>::ProblemDescription descriptor;
+
+ // First set the refinement
+ // criterion we wish to use:
+ descriptor.refinement_criterion
+ = Framework<dim>::ProblemDescription::dual_weighted_error_estimator;
+ // Here, we could as well have
+ // used ``global_refinement''
+ // or
+ // ``weighted_kelly_indicator''. Note
+ // that the information given
+ // about dual finite elements,
+ // dual functional, etc is only
+ // important for the given
+ // choice of refinement
+ // criterion, and is ignored
+ // otherwise.
+
+ // Then set the polynomial
+ // degrees of primal and dual
+ // problem. We choose here
+ // bi-linear and bi-quadratic
+ // ones:
+ descriptor.primal_fe_degree = 1;
+ descriptor.dual_fe_degree = 2;
+
+ // Then set the description of
+ // the test case, i.e. domain,
+ // boundary values, and right
+ // hand side. These are
+ // prepackaged in classes. We
+ // take here the description of
+ // ``Exercise_2_3'', but you
+ // can also use
+ // ``CurvedRidges<dim>'':
+ descriptor.data = new Data::SetUp<Data::Exercise_2_3<dim>,dim> ();
+
+ // Next set first a dual
+ // functional, then a list of
+ // evaluation objects. We
+ // choose as default the
+ // evaluation of the
+ // value at an
+ // evaluation point,
+ // represented by the classes
+ // ``PointValueEvaluation''
+ // in the namespaces of
+ // evaluation and dual
+ // functional classes. You can
+ // also set the
+ // ``PointXDerivativeEvaluation''
+ // classes for the x-derivative
+ // instead of the value
+ // at the evaluation point.
+ //
+ // Note that dual functional
+ // and evaluation objects
+ // should match. However, you
+ // can give as many evaluation
+ // functionals as you want, so
+ // you can have both point
+ // value and derivative
+ // evaluated after each step.
+ // One such additional
+ // evaluation is to output the
+ // grid in each step.
+ const Point<dim> evaluation_point (0., 0., 0.);
+ descriptor.dual_functional
+ = new DualFunctional::PointValueEvaluation<dim> (evaluation_point);
+
+ Evaluation::PointValueEvaluation<dim>
+ postprocessor1 (evaluation_point);
+ Evaluation::GridOutput<dim>
+ postprocessor2 ("grid");
+
+ descriptor.evaluator_list.push_back (&postprocessor1);
+ descriptor.evaluator_list.push_back (&postprocessor2);
+
+ // Set the maximal number of
+ // degrees of freedom after
+ // which we want the program to
+ // stop refining the mesh
+ // further:
+#if defined(SPEC_CPU)
+ // raise from 20000 to 30000. (jfk p6f)
+ descriptor.max_degrees_of_freedom = 30000;
+#else
+ descriptor.max_degrees_of_freedom = 20000;
+#endif
+
+ // Finally pass the descriptor
+ // object to a function that
+ // runs the entire solution
+ // with it:
+ Framework<dim>::run (descriptor);
+ }
+
+ // Catch exceptions to give
+ // information about things that
+ // failed:
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+
+ return 0;
+}