From 658a1b5517a42b296ce4ca0ca8a40c96acb8cba0 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 20 Sep 2012 03:09:52 +0000 Subject: [PATCH] Do a patch that has happened on mainline but failed on this branch because of local changes. git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@26545 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/benchmarks/old/spec2006-447.dealII.cc | 362 +- tests/benchmarks/spec2006-447.dealII.cc | 4275 ------------------- 2 files changed, 181 insertions(+), 4456 deletions(-) delete mode 100644 tests/benchmarks/spec2006-447.dealII.cc diff --git a/tests/benchmarks/old/spec2006-447.dealII.cc b/tests/benchmarks/old/spec2006-447.dealII.cc index 9fc45eac2c..73ca71506a 100644 --- a/tests/benchmarks/old/spec2006-447.dealII.cc +++ b/tests/benchmarks/old/spec2006-447.dealII.cc @@ -107,7 +107,7 @@ class MassCoefficient : public Function } }; - + // @sect3{Evaluating the solution} @@ -128,16 +128,16 @@ namespace Evaluation { // @sect4{The EvaluationBase class} template - class EvaluationBase + 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; }; @@ -145,9 +145,9 @@ namespace Evaluation template EvaluationBase::~EvaluationBase () {} - - + + template void EvaluationBase::set_refinement_cycle (const unsigned int step) @@ -162,15 +162,15 @@ namespace Evaluation { 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; }; @@ -181,14 +181,14 @@ namespace Evaluation : evaluation_point (evaluation_point) {} - + template void PointValueEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { double point_value = 1e20; @@ -246,15 +246,15 @@ namespace Evaluation { 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; }; @@ -265,7 +265,7 @@ namespace Evaluation : evaluation_point (evaluation_point) {} - + // The more interesting things // happen inside the function doing @@ -274,7 +274,7 @@ namespace Evaluation void PointXDerivativeEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { // This time initialize the // return value with something @@ -294,7 +294,7 @@ namespace Evaluation 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 @@ -429,7 +429,7 @@ namespace Evaluation } - + // @sect4{The GridOutput class} // Since this program has a more @@ -460,10 +460,10 @@ namespace Evaluation { 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; }; @@ -474,7 +474,7 @@ namespace Evaluation : output_name_base (output_name_base) {} - + template void @@ -496,12 +496,12 @@ namespace Evaluation #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 @@ -518,8 +518,8 @@ namespace LaplaceSolver // which requires the class to be // known: template class WeightedResidual; - - + + // @sect4{The Laplace solver base class} // This class is almost unchanged, @@ -552,8 +552,8 @@ namespace LaplaceSolver virtual void set_refinement_cycle (const unsigned int cycle); virtual void output_solution () const = 0; - - + + const SmartPointer > triangulation; unsigned int refinement_cycle; @@ -568,7 +568,7 @@ namespace LaplaceSolver template - Base::~Base () + Base::~Base () {} @@ -579,7 +579,7 @@ namespace LaplaceSolver { refinement_cycle = cycle; } - + // @sect4{The Laplace Solver class} @@ -593,7 +593,7 @@ namespace LaplaceSolver Solver (Triangulation &triangulation, const FiniteElement &fe, const Quadrature &quadrature, - const Quadrature &face_quadrature, + const Quadrature &face_quadrature, const Function &boundary_values); virtual ~Solver (); @@ -609,24 +609,24 @@ namespace LaplaceSolver virtual unsigned int n_dofs () const; - - + + const SmartPointer > fe; const SmartPointer > quadrature; - const SmartPointer > face_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; @@ -655,14 +655,14 @@ namespace LaplaceSolver Base (triangulation), fe (&fe), quadrature (&quadrature), - face_quadrature (&face_quadrature), + face_quadrature (&face_quadrature), dof_handler (triangulation), boundary_values (&boundary_values) {} template - Solver::~Solver () + Solver::~Solver () { dof_handler.clear (); } @@ -696,7 +696,7 @@ namespace LaplaceSolver { return dof_handler.n_dofs(); } - + template void @@ -708,7 +708,7 @@ namespace LaplaceSolver const unsigned int n_threads = multithread_info.n_default_threads; std::vector > - thread_ranges + thread_ranges = Threads::split_range (dof_handler.begin_active (), dof_handler.end (), n_threads); @@ -730,7 +730,7 @@ namespace LaplaceSolver 0, *boundary_values, boundary_value_map); - + threads.join_all (); linear_system.hanging_node_constraints.condense (linear_system.matrix); @@ -749,7 +749,7 @@ namespace LaplaceSolver Threads::ThreadMutex &mutex) const { MappingQ mapping (4); - FEValues fe_values (mapping, *fe, *quadrature, + FEValues fe_values (mapping, *fe, *quadrature, UpdateFlags(update_gradients | update_values | update_q_points | update_JxW_values)); @@ -759,12 +759,12 @@ namespace LaplaceSolver FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (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) { @@ -776,8 +776,8 @@ namespace LaplaceSolver laplace_coefficients); MassCoefficient().value_list (fe_values.get_quadrature_points(), mass_coefficients); - - + + for (unsigned int q_point=0; q_point &, ConstraintMatrix &) = &DoFTools::make_hanging_node_constraints; - + Threads::Thread<> mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints); @@ -928,17 +928,17 @@ namespace LaplaceSolver 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; @@ -1008,16 +1008,16 @@ namespace LaplaceSolver { abort (); } - + template void PrimalSolver:: - assemble_rhs (Vector &rhs) const + assemble_rhs (Vector &rhs) const { MappingQ mapping (4); - FEValues fe_values (mapping, *this->fe, *this->quadrature, + FEValues fe_values (mapping, *this->fe, *this->quadrature, UpdateFlags(update_values | update_q_points | update_JxW_values)); @@ -1027,7 +1027,7 @@ namespace LaplaceSolver Vector cell_rhs (dofs_per_cell); std::vector rhs_values (n_q_points); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); typename DoFHandler::active_cell_iterator cell = this->dof_handler.begin_active(), @@ -1040,7 +1040,7 @@ namespace LaplaceSolver rhs_function->value_list (fe_values.get_quadrature_points(), rhs_values); - + for (unsigned int q_point=0; q_point > weighting_function; }; @@ -1247,7 +1247,7 @@ namespace LaplaceSolver 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); @@ -1392,7 +1392,7 @@ namespace LaplaceSolver namespace Data { // @sect4{The SetUpBase and SetUp classes} - + // Based on the above description, // the ``SetUpBase'' class then // looks as follows. To allow using @@ -1436,12 +1436,12 @@ namespace Data 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; }; @@ -1458,7 +1458,7 @@ namespace Data // functions: template const Function & - SetUp::get_boundary_values () const + SetUp::get_boundary_values () const { return boundary_values; } @@ -1466,7 +1466,7 @@ namespace Data template const Function & - SetUp::get_right_hand_side () const + SetUp::get_right_hand_side () const { return right_hand_side; } @@ -1475,11 +1475,11 @@ namespace Data template void SetUp:: - create_coarse_grid (Triangulation &coarse_grid) const + create_coarse_grid (Triangulation &coarse_grid) const { Traits::create_coarse_grid (coarse_grid); } - + // @sect4{The CurvedRidges class} @@ -1496,7 +1496,7 @@ namespace Data { public: BoundaryValues () : Function () {}; - + virtual double value (const Point &p, const unsigned int component) const; }; @@ -1506,7 +1506,7 @@ namespace Data { public: RightHandSide () : Function () {}; - + virtual double value (const Point &p, const unsigned int component) const; }; @@ -1515,8 +1515,8 @@ namespace Data void create_coarse_grid (Triangulation &coarse_grid); }; - - + + template double CurvedRidges::BoundaryValues:: @@ -1553,7 +1553,7 @@ namespace Data 100*std::sin(10*p(i)+5*p(0)*p(0)); }; t1 = t1*t1; - + return -u*(t1+t2+t3); } @@ -1566,10 +1566,10 @@ namespace Data 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 @@ -1611,7 +1611,7 @@ namespace Data public: RightHandSide () : ConstantFunction (1.) {}; }; - + // Finally a function to // generate the coarse // grid. This is somewhat more @@ -1651,7 +1651,7 @@ namespace Data // enough to do the creation by // hand, rather than using a mesh // generator. -/* +/* template <> void Exercise_2_3<2>:: @@ -1686,27 +1686,27 @@ namespace Data 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> (0., 1.), Point<2> (+1./2, 1.), Point<2> (+1, 1.) }; const unsigned int @@ -1750,7 +1750,7 @@ namespace Data // the material indicator to zero // for all the cells: std::vector > cells (n_cells, CellData()); - for (unsigned int i=0; i::vertices_per_cell; @@ -1772,14 +1772,14 @@ namespace Data 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 <> @@ -1792,7 +1792,7 @@ namespace Data coarse_grid.set_boundary (0, boundary); coarse_grid.refine_global (1); } - + } // @sect4{Discussion} @@ -1908,7 +1908,7 @@ namespace Data namespace DualFunctional { // @sect4{The DualFunctionalBase class} - + // First start with the base class // for dual functionals. Since for // linear problems the @@ -1930,7 +1930,7 @@ namespace DualFunctional // @sect4{The PointValueEvaluation class} - + // As a first application, we // consider the functional // corresponding to the evaluation @@ -1952,13 +1952,13 @@ namespace DualFunctional 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; }; @@ -1969,7 +1969,7 @@ namespace DualFunctional : evaluation_point (evaluation_point) {} - + // As for doing the main purpose of // the class, assembling the right @@ -2042,7 +2042,7 @@ namespace DualFunctional // @sect4{The PointXDerivativeEvaluation class} - + // As second application, we again // consider the evaluation of the // x-derivative of the solution at @@ -2061,13 +2061,13 @@ namespace DualFunctional 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; }; @@ -2078,7 +2078,7 @@ namespace DualFunctional : evaluation_point (evaluation_point) {} - + // What is interesting is the // implementation of this @@ -2130,7 +2130,7 @@ namespace DualFunctional // for the number of quadrature // points and shape functions... QGauss quadrature(4); - MappingQ mapping (4); + MappingQ mapping (4); FEValues fe_values (mapping, dof_handler.get_fe(), quadrature, update_gradients | update_q_points | @@ -2146,7 +2146,7 @@ namespace DualFunctional // shape functions at the // quadrature points: Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); // Finally have a variable in // which we will sum up the @@ -2155,7 +2155,7 @@ namespace DualFunctional // 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 @@ -2178,7 +2178,7 @@ namespace DualFunctional // total area/volume. fe_values.reinit (cell); cell_rhs = 0; - + for (unsigned int q=0; q &postprocessor) const; - + const SmartPointer > dual_functional; virtual void assemble_rhs (Vector &rhs) const; @@ -2340,13 +2340,13 @@ namespace LaplaceSolver { Solver::postprocess(postprocessor); } - + template void DualSolver:: - assemble_rhs (Vector &rhs) const + assemble_rhs (Vector &rhs) const { dual_functional->assemble_rhs (this->dof_handler, rhs); } @@ -2392,7 +2392,7 @@ namespace LaplaceSolver virtual void postprocess (const Evaluation::EvaluationBase &postprocessor) const; - + virtual unsigned int n_dofs () const; @@ -2403,7 +2403,7 @@ namespace LaplaceSolver void output_solution () const; - + // In the private section, we // have two functions that are // used to call the @@ -2420,7 +2420,7 @@ namespace LaplaceSolver // to avoid that we have to // write this lengthy name // over and over again: - + typedef typename DoFHandler::active_cell_iterator active_cell_iterator; @@ -2549,8 +2549,8 @@ namespace LaplaceSolver const SmartPointer > right_hand_side; std::vector cell_residual; - std::vector rhs_values; - std::vector dual_weights; + std::vector rhs_values; + std::vector dual_weights; typename std::vector > cell_grad_grads; CellData (const FiniteElement &fe, const Quadrature &quadrature, @@ -2565,14 +2565,14 @@ namespace LaplaceSolver FESubfaceValues fe_subface_values_cell; std::vector jump_residual; - std::vector dual_weights; + 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 @@ -2652,17 +2652,17 @@ namespace LaplaceSolver 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); + rhs_values.resize(n_q_points); + dual_weights.resize(n_q_points); cell_grad_grads.resize(n_q_points); } - - + + template WeightedResidual::FaceData:: @@ -2682,16 +2682,16 @@ namespace LaplaceSolver 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); + dual_weights.resize(n_face_q_points); cell_grads.resize(n_face_q_points); neighbor_grads.resize(n_face_q_points); } - + @@ -2735,7 +2735,7 @@ namespace LaplaceSolver threads.join_all (); } - + template void WeightedResidual::solve_primal_problem () @@ -2749,7 +2749,7 @@ namespace LaplaceSolver { DualSolver::solve_problem (); } - + template void @@ -2758,8 +2758,8 @@ namespace LaplaceSolver { PrimalSolver::postprocess (postprocessor); } - - + + template unsigned int WeightedResidual::n_dofs () const @@ -2809,11 +2809,11 @@ namespace LaplaceSolver // of the error. GridRefinement::refine_and_coarsen_fixed_fraction (*this->triangulation, error_indicators, - 0.5/std::sqrt(1.+step), + 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 @@ -2905,7 +2905,7 @@ namespace LaplaceSolver { 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); @@ -2915,7 +2915,7 @@ namespace LaplaceSolver dual_solver.solution, primal_solver.dof_handler, primal_hanging_node_constraints, - dual_solution); + dual_solution); // approximate error, gradient, // and second derivative @@ -2936,13 +2936,13 @@ namespace LaplaceSolver 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()); @@ -2957,8 +2957,8 @@ namespace LaplaceSolver 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 @@ -2996,20 +2996,20 @@ namespace LaplaceSolver { 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); @@ -3019,7 +3019,7 @@ namespace LaplaceSolver 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 @@ -3035,10 +3035,10 @@ namespace LaplaceSolver #else std::ofstream out (filename.str()); #endif - + data_out.write_gmv (out); } - + } @@ -3053,7 +3053,7 @@ namespace LaplaceSolver // functions that actually do the // work, and finally collects the // results. - + template void WeightedResidual:: @@ -3105,7 +3105,7 @@ namespace LaplaceSolver dual_solver.dof_handler, dual_hanging_node_constraints, primal_solution); - + // Then for computing the // interpolation of the // numerically approximated dual @@ -3128,7 +3128,7 @@ namespace LaplaceSolver primal_solver.dof_handler, primal_hanging_node_constraints, dual_weights); - + // Note that this could probably // have been more efficient since // those constraints have been @@ -3139,7 +3139,7 @@ namespace LaplaceSolver // 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 @@ -3206,7 +3206,7 @@ namespace LaplaceSolver n_threads, i, error_indicators, face_integrals); - threads.join_all(); + threads.join_all(); // Once the error contributions // are computed, sum them up. For @@ -3222,7 +3222,7 @@ namespace LaplaceSolver // term, since the other half // will be taken by the // neighboring cell. - unsigned int present_cell=0; + 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) @@ -3284,7 +3284,7 @@ namespace LaplaceSolver *dual_solver.quadrature, *primal_solver.rhs_function); FaceData face_data (*dual_solver.fe, - *dual_solver.face_quadrature); + *dual_solver.face_quadrature); // Then calculate the start cell // for this thread. We let the @@ -3320,7 +3320,7 @@ namespace LaplaceSolver // back right now if (cell == dual_solver.dof_handler.end()) return; - + // Next loop over all cells. The // check for loop end is done at // the end of the loop, along @@ -3340,7 +3340,7 @@ namespace LaplaceSolver dual_weights, cell_data, error_indicators); - + // After computing the cell // terms, turn to the face // terms. For this, loop over @@ -3364,12 +3364,12 @@ namespace LaplaceSolver // the list of faces with // a zero contribution to // the error. - if (cell->face(face_no)->at_boundary()) + 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 @@ -3413,7 +3413,7 @@ namespace LaplaceSolver // the coarse cell. if (cell->at_boundary(face_no) == false) if (cell->neighbor(face_no)->level() < cell->level()) - continue; + continue; // Now we know that we @@ -3443,7 +3443,7 @@ namespace LaplaceSolver primal_solution, dual_weights, face_data, - face_integrals); + face_integrals); else integrate_over_irregular_face (cell, face_no, primal_solution, @@ -3518,7 +3518,7 @@ namespace LaplaceSolver // @sect4{Computing edge term error contributions - 1} - + // On the other hand, computation // of the edge terms for the error // estimate is not so @@ -3595,7 +3595,7 @@ namespace LaplaceSolver // 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.reinit (neighbor, neighbor_neighbor); face_data.fe_face_values_neighbor.get_function_grads (primal_solution, face_data.neighbor_grads); @@ -3614,7 +3614,7 @@ namespace LaplaceSolver // this face: face_data.fe_face_values_cell.get_function_values (dual_weights, face_data.dual_weights); - + // Finally, we have to compute // the sum over jump residuals, // dual weights, and quadrature @@ -3657,7 +3657,7 @@ namespace LaplaceSolver // @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 @@ -3680,7 +3680,7 @@ namespace LaplaceSolver n_q_points = face_data.fe_face_values_cell.n_quadrature_points; const typename DoFHandler::cell_iterator - neighbor = cell->neighbor(face_no); + neighbor = cell->neighbor(face_no); Assert (neighbor.state() == IteratorState::valid, ExcInternalError()); Assert (neighbor->has_children(), @@ -3697,7 +3697,7 @@ namespace LaplaceSolver // 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 @@ -3722,7 +3722,7 @@ namespace LaplaceSolver 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 @@ -3738,7 +3738,7 @@ namespace LaplaceSolver // Then get dual weights: face_data.fe_face_values_neighbor.get_function_values (dual_weights, face_data.dual_weights); - + // At last, sum up the // contribution of this // sub-face, and set it in @@ -3767,21 +3767,21 @@ namespace LaplaceSolver typename DoFHandler::face_iterator face = cell->face(face_no); for (unsigned int subface_no=0; subface_no::subfaces_per_face; - ++subface_no) + ++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; } - + } @@ -3837,7 +3837,7 @@ struct Framework // structure, these values are // all set to default values, for // simple use. - struct ProblemDescription + struct ProblemDescription { // First allow for the // degrees of the piecewise @@ -4006,7 +4006,7 @@ void Framework::run (const ProblemDescription &descriptor) *descriptor.dual_functional); break; }; - + case ProblemDescription::global_refinement: { solver @@ -4018,7 +4018,7 @@ void Framework::run (const ProblemDescription &descriptor) descriptor.data->get_boundary_values()); break; }; - + case ProblemDescription::kelly_indicator: { solver @@ -4043,11 +4043,11 @@ void Framework::run (const ProblemDescription &descriptor) *descriptor.kelly_weight); break; }; - + default: AssertThrow (false, ExcInternalError()); }; - + // Now that all objects are in // place, run the main loop. The // stopping criterion is @@ -4066,14 +4066,14 @@ void Framework::run (const ProblemDescription &descriptor) { 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) @@ -4082,7 +4082,7 @@ void Framework::run (const ProblemDescription &descriptor) solver->postprocess (**e); }; - + if (solver->n_dofs() < descriptor.max_degrees_of_freedom) solver->refine_grid (); else @@ -4111,7 +4111,7 @@ void Framework::run (const ProblemDescription &descriptor) // etc), and passes them packed into // a structure to the frame work // class above. -int main (int argc, char **argv) +int main (int argc, char **argv) { // if no argument is given, then do 18 // iterations @@ -4134,10 +4134,10 @@ int main (int argc, char **argv) << std::endl; exit (1); } - - - - + + + + deallog.depth_console (0); try { @@ -4183,7 +4183,7 @@ int main (int argc, char **argv) // can also use // ``CurvedRidges'': descriptor.data = new Data::SetUp,dim> (); - + // Next set first a dual // functional, then a list of // evaluation objects. We @@ -4216,12 +4216,12 @@ int main (int argc, char **argv) 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); @@ -4236,7 +4236,7 @@ int main (int argc, char **argv) #else descriptor.max_degrees_of_freedom = 20000; #endif - + // Finally pass the descriptor // object to a function that // runs the entire solution @@ -4259,7 +4259,7 @@ int main (int argc, char **argv) << std::endl; return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" diff --git a/tests/benchmarks/spec2006-447.dealII.cc b/tests/benchmarks/spec2006-447.dealII.cc deleted file mode 100644 index 7801034c6d..0000000000 --- a/tests/benchmarks/spec2006-447.dealII.cc +++ /dev/null @@ -1,4275 +0,0 @@ -/* $Id$ */ -/* Author: Wolfgang Bangerth, ETH Zurich, 2002 */ - -/* $Id$ */ -/* Version: $Name$ */ -/* */ -/* Copyright (C) 2002, 2003, 2004, 2008, 2010 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, - QGauss(3), - 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... - QGauss quadrature(4); - 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, - QGauss(3), - 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