From 5e6bd84812f82655f1e956e0b73104259cf44fd5 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 25 Jan 2008 00:00:18 +0000 Subject: [PATCH] Add the SPEC CPU2006 benchmark 447.dealII with 2 or 3 minor modifications to make sure it still compiles. git-svn-id: https://svn.dealii.org/trunk@15681 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/benchmarks/spec2006-447.dealII.cc | 4275 +++++++++++++++++++++++ 1 file changed, 4275 insertions(+) create mode 100644 tests/benchmarks/spec2006-447.dealII.cc diff --git a/tests/benchmarks/spec2006-447.dealII.cc b/tests/benchmarks/spec2006-447.dealII.cc new file mode 100644 index 0000000000..617b60d228 --- /dev/null +++ b/tests/benchmarks/spec2006-447.dealII.cc @@ -0,0 +1,4275 @@ +/* $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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#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 +class LaplaceCoefficient : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int) const + { + return 1.+p*p; + } +}; + + + /** + * Same for the coefficient in front + * of the mass term. + */ +template +class MassCoefficient : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int) const + { + double factor = 1; + for (unsigned int d=0; d + class EvaluationBase + { + public: + virtual ~EvaluationBase (); + + void set_refinement_cycle (const unsigned int refinement_cycle); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const = 0; + + unsigned int refinement_cycle; + }; + + + template + EvaluationBase::~EvaluationBase () + {} + + + + template + void + EvaluationBase::set_refinement_cycle (const unsigned int step) + { + refinement_cycle = step; + } + + + // @sect4{The PointValueEvaluation class} + template + class PointValueEvaluation : public EvaluationBase + { + public: + PointValueEvaluation (const Point &evaluation_point); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + const Point evaluation_point; + }; + + + template + PointValueEvaluation:: + PointValueEvaluation (const Point &evaluation_point) + : + evaluation_point (evaluation_point) + {} + + + + template + void + PointValueEvaluation:: + operator () (const DoFHandler &dof_handler, + const Vector &solution) const + { + double point_value = 1e20; + + typename DoFHandler::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::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 + class PointXDerivativeEvaluation : public EvaluationBase + { + public: + PointXDerivativeEvaluation (const Point &evaluation_point); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + const Point evaluation_point; + }; + + + template + PointXDerivativeEvaluation:: + PointXDerivativeEvaluation (const Point &evaluation_point) + : + evaluation_point (evaluation_point) + {} + + + // The more interesting things + // happen inside the function doing + // the actual evaluation: + template + void + PointXDerivativeEvaluation:: + operator () (const DoFHandler &dof_handler, + const Vector &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 vertex_quadrature; + MappingQ mapping (4); + FEValues fe_values (mapping, dof_handler.get_fe(), + vertex_quadrature, + update_gradients | update_q_points); + std::vector > + 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::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::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 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 + class GridOutput : public EvaluationBase + { + public: + GridOutput (const std::string &output_name_base); + + virtual void operator () (const DoFHandler &dof_handler, + const Vector &solution) const; + + const std::string output_name_base; + }; + + + template + GridOutput:: + GridOutput (const std::string &output_name_base) + : + output_name_base (output_name_base) + {} + + + template + void + GridOutput::operator () (const DoFHandler &dof_handler, + const Vector &/*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 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 + class Base + { + public: + Base (Triangulation &coarse_grid); + virtual ~Base (); + + virtual void solve_problem () = 0; + virtual void postprocess (const Evaluation::EvaluationBase &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; + + unsigned int refinement_cycle; + }; + + + template + Base::Base (Triangulation &coarse_grid) + : + triangulation (&coarse_grid) + {} + + + template + Base::~Base () + {} + + + + template + void + Base::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 + class Solver : public virtual Base + { + public: + Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &boundary_values); + virtual + ~Solver (); + + virtual + void + solve_problem (); + + virtual + void + postprocess (const Evaluation::EvaluationBase &postprocessor) const; + + virtual + unsigned int + n_dofs () const; + + + const SmartPointer > fe; + const SmartPointer > quadrature; + const SmartPointer > face_quadrature; + DoFHandler dof_handler; + Vector solution; + const SmartPointer > boundary_values; + + virtual void assemble_rhs (Vector &rhs) const = 0; + + + struct LinearSystem + { + LinearSystem (const DoFHandler &dof_handler); + + void solve (Vector &solution) const; + + ConstraintMatrix hanging_node_constraints; + SparsityPattern sparsity_pattern; + SparseMatrix matrix; + Vector rhs; + }; + + void + assemble_linear_system (LinearSystem &linear_system); + + void + assemble_matrix (LinearSystem &linear_system, + const typename DoFHandler::active_cell_iterator &begin_cell, + const typename DoFHandler::active_cell_iterator &end_cell, + Threads::ThreadMutex &mutex) const; + }; + + + + template + Solver::Solver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &boundary_values) + : + Base (triangulation), + fe (&fe), + quadrature (&quadrature), + face_quadrature (&face_quadrature), + dof_handler (triangulation), + boundary_values (&boundary_values) + {} + + + template + Solver::~Solver () + { + dof_handler.clear (); + } + + + template + void + Solver::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 + void + Solver:: + postprocess (const Evaluation::EvaluationBase &postprocessor) const + { + postprocessor (dof_handler, solution); + } + + + template + unsigned int + Solver::n_dofs () const + { + return dof_handler.n_dofs(); + } + + + template + void + Solver::assemble_linear_system (LinearSystem &linear_system) + { + typedef + typename DoFHandler::active_cell_iterator + active_cell_iterator; + + const unsigned int n_threads = multithread_info.n_default_threads; + std::vector > + thread_ranges + = Threads::split_range (dof_handler.begin_active (), + dof_handler.end (), + n_threads); + + Threads::ThreadMutex mutex; + Threads::ThreadGroup<> threads; + for (unsigned int thread=0; thread::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 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 + void + Solver::assemble_matrix (LinearSystem &linear_system, + const typename DoFHandler::active_cell_iterator &begin_cell, + const typename DoFHandler::active_cell_iterator &end_cell, + Threads::ThreadMutex &mutex) const + { + MappingQ mapping (4); + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + std::vector laplace_coefficients (fe_values.n_quadrature_points); + std::vector mass_coefficients (fe_values.n_quadrature_points); + + + for (typename DoFHandler::active_cell_iterator cell=begin_cell; + cell!=end_cell; ++cell) + { + cell_matrix = 0; + + fe_values.reinit (cell); + + LaplaceCoefficient().value_list (fe_values.get_quadrature_points(), + laplace_coefficients); + MassCoefficient().value_list (fe_values.get_quadrature_points(), + mass_coefficients); + + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + Threads::ThreadMutex::ScopedLock lock (mutex); + for (unsigned int i=0; i + Solver::LinearSystem:: + LinearSystem (const DoFHandler &dof_handler) + { + hanging_node_constraints.clear (); + + void (*mhnc_p) (const DoFHandler &, + 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 + void + Solver::LinearSystem::solve (Vector &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::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 + class PrimalSolver : public Solver + { + public: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values); + + virtual + void solve_problem (); + + virtual + unsigned int n_dofs () const; + + virtual + void postprocess (const Evaluation::EvaluationBase &postprocessor) const; + + virtual + void output_solution () const; + + + const SmartPointer > rhs_function; + virtual void assemble_rhs (Vector &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; + }; + + + template + PrimalSolver:: + PrimalSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (triangulation), + Solver (triangulation, fe, + quadrature, face_quadrature, + boundary_values), + rhs_function (&rhs_function) + {} + + + template + void + PrimalSolver::solve_problem () + { + Solver::solve_problem (); + } + + + + template + unsigned int + PrimalSolver::n_dofs() const + { + return Solver::n_dofs(); + } + + + template + void + PrimalSolver:: + postprocess (const Evaluation::EvaluationBase &postprocessor) const + { + Solver::postprocess(postprocessor); + } + + + template + void + PrimalSolver::output_solution () const + { + abort (); + } + + + + template + void + PrimalSolver:: + assemble_rhs (Vector &rhs) const + { + MappingQ mapping (4); + FEValues 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 cell_rhs (dofs_per_cell); + std::vector rhs_values (n_q_points); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i + class RefinementGlobal : public PrimalSolver + { + public: + RefinementGlobal (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values); + + virtual void refine_grid (); + }; + + + + template + RefinementGlobal:: + RefinementGlobal (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + face_quadrature, rhs_function, + boundary_values) + {} + + + + template + void + RefinementGlobal::refine_grid () + { + this->triangulation->refine_global (1); + } + + + + template + class RefinementKelly : public PrimalSolver + { + public: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values); + + virtual void refine_grid (); + }; + + + + template + RefinementKelly:: + RefinementKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + face_quadrature, + rhs_function, boundary_values) + {} + + + + template + void + RefinementKelly::refine_grid () + { + Vector estimated_error_per_cell (this->triangulation->n_active_cells()); + KellyErrorEstimator::estimate (this->dof_handler, + QGauss3(), + typename FunctionMap::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 + class RefinementWeightedKelly : public PrimalSolver + { + public: + RefinementWeightedKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values, + const Function &weighting_function); + + virtual void refine_grid (); + + + const SmartPointer > weighting_function; + }; + + + + template + RefinementWeightedKelly:: + RefinementWeightedKelly (Triangulation &coarse_grid, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values, + const Function &weighting_function) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, fe, quadrature, + face_quadrature, + rhs_function, boundary_values), + weighting_function (&weighting_function) + {} + + + + // Now, here comes the main + // function, including the + // weighting: + template + void + RefinementWeightedKelly::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 estimated_error (this->triangulation->n_active_cells()); + KellyErrorEstimator::estimate (this->dof_handler, + *this->face_quadrature, + typename FunctionMap::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::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'', + // 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 + struct SetUpBase : public Subscriptor + { + virtual + const Function & get_boundary_values () const = 0; + + virtual + const Function & get_right_hand_side () const = 0; + + virtual + void create_coarse_grid (Triangulation &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 + struct SetUp : public SetUpBase + { + SetUp () {}; + + virtual + const Function & get_boundary_values () const; + + virtual + const Function & get_right_hand_side () const; + + + virtual + void create_coarse_grid (Triangulation &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 + const typename Traits::BoundaryValues SetUp::boundary_values; + template + const typename Traits::RightHandSide SetUp::right_hand_side; + + // And definitions of the member + // functions: + template + const Function & + SetUp::get_boundary_values () const + { + return boundary_values; + } + + + template + const Function & + SetUp::get_right_hand_side () const + { + return right_hand_side; + } + + + template + void + SetUp:: + create_coarse_grid (Triangulation &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 + struct CurvedRidges + { + class BoundaryValues : public Function + { + public: + BoundaryValues () : Function () {}; + + virtual double value (const Point &p, + const unsigned int component) const; + }; + + + class RightHandSide : public Function + { + public: + RightHandSide () : Function () {}; + + virtual double value (const Point &p, + const unsigned int component) const; + }; + + static + void + create_coarse_grid (Triangulation &coarse_grid); + }; + + + template + double + CurvedRidges::BoundaryValues:: + value (const Point &p, + const unsigned int /*component*/) const + { + double q = p(0); + for (unsigned int i=1; i + double + CurvedRidges::RightHandSide::value (const Point &p, + const unsigned int /*component*/) const + { + double q = p(0); + for (unsigned int i=1; i + void + CurvedRidges:: + create_coarse_grid (Triangulation &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 + 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 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 + { + public: + RightHandSide () : ConstantFunction (1.) {}; + }; + + // Finally a function to + // generate the coarse + // grid. This is somewhat more + // complicated here, see + // immediately below. + static + void + create_coarse_grid (Triangulation &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 > 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::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 > cells (n_cells, CellData()); + for (unsigned int i=0; i::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'' where + // ``SomeSetUp'' denotes a class, or + // even to explicitly specialize + // ``Data::SetUp''. 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 + class DualFunctionalBase : public Subscriptor + { + public: + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + Vector &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 + class PointValueEvaluation : public DualFunctionalBase + { + public: + PointValueEvaluation (const Point &evaluation_point); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + Vector &rhs) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + + const Point evaluation_point; + }; + + + template + PointValueEvaluation:: + PointValueEvaluation (const Point &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 + void + PointValueEvaluation:: + assemble_rhs (const DoFHandler &dof_handler, + Vector &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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; + vertex::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 + class PointXDerivativeEvaluation : public DualFunctionalBase + { + public: + PointXDerivativeEvaluation (const Point &evaluation_point); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + Vector &rhs) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + + const Point evaluation_point; + }; + + + template + PointXDerivativeEvaluation:: + PointXDerivativeEvaluation (const Point &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 + void + PointXDerivativeEvaluation:: + assemble_rhs (const DoFHandler &dof_handler, + Vector &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 quadrature; + MappingQ mapping (4); + FEValues 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 cell_rhs (dofs_per_cell); + std::vector 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::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; qget_dof_indices (local_dof_indices); + for (unsigned int i=0; i 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 + class DualSolver : public Solver + { + public: + DualSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const DualFunctional::DualFunctionalBase &dual_functional); + + virtual + void + solve_problem (); + + virtual + unsigned int + n_dofs () const; + + virtual + void + postprocess (const Evaluation::EvaluationBase &postprocessor) const; + + + const SmartPointer > dual_functional; + virtual void assemble_rhs (Vector &rhs) const; + + static const ZeroFunction boundary_values; + + // Same as above -- make a + // derived class a friend of + // this one: + friend class WeightedResidual; + }; + + template + const ZeroFunction DualSolver::boundary_values; + + template + DualSolver:: + DualSolver (Triangulation &triangulation, + const FiniteElement &fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const DualFunctional::DualFunctionalBase &dual_functional) + : + Base (triangulation), + Solver (triangulation, fe, + quadrature, face_quadrature, + boundary_values), + dual_functional (&dual_functional) + {} + + + template + void + DualSolver::solve_problem () + { + Solver::solve_problem (); + } + + + + template + unsigned int + DualSolver::n_dofs() const + { + return Solver::n_dofs(); + } + + + template + void + DualSolver:: + postprocess (const Evaluation::EvaluationBase &postprocessor) const + { + Solver::postprocess(postprocessor); + } + + + + template + void + DualSolver:: + assemble_rhs (Vector &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 + class WeightedResidual : public PrimalSolver, + public DualSolver + { + public: + WeightedResidual (Triangulation &coarse_grid, + const FiniteElement &primal_fe, + const FiniteElement &dual_fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &boundary_values, + const DualFunctional::DualFunctionalBase &dual_functional); + + virtual + void + solve_problem (); + + virtual + void + postprocess (const Evaluation::EvaluationBase &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::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::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 mapping; + FEValues fe_values; + const SmartPointer > right_hand_side; + + std::vector cell_residual; + std::vector rhs_values; + std::vector dual_weights; + typename std::vector > cell_grad_grads; + CellData (const FiniteElement &fe, + const Quadrature &quadrature, + const Function &right_hand_side); + }; + + struct FaceData + { + MappingQ mapping; + FEFaceValues fe_face_values_cell; + FEFaceValues fe_face_values_neighbor; + FESubfaceValues fe_subface_values_cell; + + std::vector jump_residual; + std::vector dual_weights; + typename std::vector > cell_grads; + typename std::vector > neighbor_grads; + FaceData (const FiniteElement &fe, + const Quadrature &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 &error_indicators) const; + + void estimate_some (const Vector &primal_solution, + const Vector &dual_weights, + const unsigned int n_threads, + const unsigned int this_thread, + Vector &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 &primal_solution, + const Vector &dual_weights, + CellData &cell_data, + Vector &error_indicators) const; + + void + integrate_over_regular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const Vector &primal_solution, + const Vector &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 &primal_solution, + const Vector &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 + WeightedResidual::CellData:: + CellData (const FiniteElement &fe, + const Quadrature &quadrature, + const Function &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 + WeightedResidual::FaceData:: + FaceData (const FiniteElement &fe, + const Quadrature &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 + WeightedResidual:: + WeightedResidual (Triangulation &coarse_grid, + const FiniteElement &primal_fe, + const FiniteElement &dual_fe, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const Function &rhs_function, + const Function &bv, + const DualFunctional::DualFunctionalBase &dual_functional) + : + Base (coarse_grid), + PrimalSolver (coarse_grid, primal_fe, + quadrature, face_quadrature, + rhs_function, bv), + DualSolver (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 + void + WeightedResidual::solve_problem () + { + Threads::ThreadGroup<> threads; + threads += Threads::spawn (*this, &WeightedResidual::solve_primal_problem)(); + threads += Threads::spawn (*this, &WeightedResidual::solve_dual_problem)(); + threads.join_all (); + } + + + template + void + WeightedResidual::solve_primal_problem () + { + PrimalSolver::solve_problem (); + } + + template + void + WeightedResidual::solve_dual_problem () + { + DualSolver::solve_problem (); + } + + + template + void + WeightedResidual:: + postprocess (const Evaluation::EvaluationBase &postprocessor) const + { + PrimalSolver::postprocess (postprocessor); + } + + + template + unsigned int + WeightedResidual::n_dofs () const + { + return PrimalSolver::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 + void + WeightedResidual::refine_grid () + { + // First call the function that + // computes the cell-wise and + // global error: + Vector 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::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 + void + WeightedResidual::output_solution () const + { + const PrimalSolver &primal_solver = *this; + const DualSolver &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 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 error_indicators (this->triangulation->n_active_cells()); + Vector gradient_indicators (this->triangulation->n_active_cells()); + Vector second_indicators (this->triangulation->n_active_cells()); + { + MappingQ mapping(4); + KellyErrorEstimator::estimate (mapping, primal_solver.dof_handler, + QGauss3(), + typename FunctionMap::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 x_error_indicators (primal_solver.dof_handler.n_dofs()); + Vector x_gradient_indicators (primal_solver.dof_handler.n_dofs()); + Vector 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 coarser_mesh; + coarser_mesh.copy_triangulation (*this->triangulation); + for (typename Triangulation::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 coarser_dof_handler (coarser_mesh); + coarser_dof_handler.distribute_dofs (primal_solver.dof_handler.get_fe()); + InterGridMap > 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 coarse_primal_solution (coarser_dof_handler.n_dofs()); + Vector coarse_dual_solution (coarser_dof_handler.n_dofs()); + Vector coarse_error_indicators (coarser_dof_handler.n_dofs()); + Vector coarse_gradient_indicators (coarser_dof_handler.n_dofs()); + Vector coarse_second_indicators (coarser_dof_handler.n_dofs()); + + Vector tmp (coarser_dof_handler.get_fe().dofs_per_cell); + for (typename DoFHandler::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 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 + void + WeightedResidual:: + estimate_error (Vector &error_indicators) const + { + const PrimalSolver &primal_solver = *this; + const DualSolver &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 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 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::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::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::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 + void + WeightedResidual:: + estimate_some (const Vector &primal_solution, + const Vector &dual_weights, + const unsigned int n_threads, + const unsigned int this_thread, + Vector &error_indicators, + FaceIntegrals &face_integrals) const + { + const PrimalSolver &primal_solver = *this; + const DualSolver &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::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 + void WeightedResidual:: + integrate_over_cell (const active_cell_iterator &cell, + const unsigned int cell_index, + const Vector &primal_solution, + const Vector &dual_weights, + CellData &cell_data, + Vector &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 + void WeightedResidual:: + integrate_over_regular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const Vector &primal_solution, + const Vector &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; pface(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 + void WeightedResidual:: + integrate_over_irregular_face (const active_cell_iterator &cell, + const unsigned int face_no, + const Vector &primal_solution, + const Vector &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::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::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; pface(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::face_iterator face = cell->face(face_no); + for (unsigned int subface_no=0; + subface_no::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 +struct Framework +{ + public: + // First, we declare two + // abbreviations for simple use + // of the respective data types: + typedef Evaluation::EvaluationBase Evaluator; + typedef std::list 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 > 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 > 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 > 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 +Framework::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 +void Framework::run (const ProblemDescription &descriptor) +{ + // First create a triangulation + // from the given data object, + Triangulation + triangulation (Triangulation::smoothing_on_refinement); + descriptor.data->create_coarse_grid (triangulation); + + // then a set of finite elements + // and appropriate quadrature + // formula: + const FE_Q primal_fe(descriptor.primal_fe_degree); + const FE_Q dual_fe(descriptor.dual_fe_degree); + const QGauss quadrature(descriptor.dual_fe_degree+1); + const QGauss face_quadrature(descriptor.dual_fe_degree+1); + + // Next, select one of the classes + // implementing different + // refinement criteria. + LaplaceSolver::Base * solver = 0; + switch (descriptor.refinement_criterion) + { + case ProblemDescription::dual_weighted_error_estimator: + { + solver + = new LaplaceSolver::WeightedResidual (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 (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 (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 (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::ProblemDescription descriptor; + + // First set the refinement + // criterion we wish to use: + descriptor.refinement_criterion + = Framework::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'': + descriptor.data = new Data::SetUp,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 evaluation_point (0., 0., 0.); + descriptor.dual_functional + = new DualFunctional::PointValueEvaluation (evaluation_point); + + Evaluation::PointValueEvaluation + postprocessor1 (evaluation_point); + Evaluation::GridOutput + 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::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; +} -- 2.39.5