From d4a17c8966573b9e4794e1337fa613168c33473b Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 22 Apr 2002 08:33:43 +0000 Subject: [PATCH] Move slightly forward. git-svn-id: https://svn.dealii.org/trunk@5701 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-14/Makefile | 4 +- deal.II/examples/step-14/step-14.cc | 1161 ++++++++++++++++++++++++--- 2 files changed, 1051 insertions(+), 114 deletions(-) diff --git a/deal.II/examples/step-14/Makefile b/deal.II/examples/step-14/Makefile index c02ef0f148..2780e00497 100644 --- a/deal.II/examples/step-14/Makefile +++ b/deal.II/examples/step-14/Makefile @@ -148,8 +148,8 @@ Makefile.dep: $(target).cc Makefile \ $(include-path-lac)/lac/*.h \ $(include-path-deal2)/*/*.h) @echo ============================ Remaking Makefile - @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ - | perl -pi -e 's!lib/g?o/!!g;' \ + @$(PERL) $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ + | $(PERL) -pi -e 's!lib/g?o/!!g;' \ > Makefile.dep # To make the dependencies known to `make', we finally have to include diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index 0283a294e0..67ee9296b0 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -42,6 +43,8 @@ #include #include #include +#include +#include #ifdef HAVE_STD_STRINGSTREAM # include @@ -146,66 +149,6 @@ namespace Evaluation results_table.add_value ("u(x_0)", point_value); }; - - - - - template - class SolutionOutput : public EvaluationBase - { - public: - SolutionOutput (const std::string &output_name_base, - const typename DataOut::OutputFormat output_format); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - private: - const std::string output_name_base; - const typename DataOut::OutputFormat output_format; - }; - - - template - SolutionOutput:: - SolutionOutput (const std::string &output_name_base, - const typename DataOut::OutputFormat output_format) - : - output_name_base (output_name_base), - output_format (output_format) - {}; - - - template - void - SolutionOutput::operator () (const DoFHandler &dof_handler, - const Vector &solution) const - { - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, "solution"); - data_out.build_patches (); - -#ifdef HAVE_STD_STRINGSTREAM - std::ostringstream filename; -#else - std::ostrstream filename; -#endif - filename << output_name_base << "-" - << refinement_cycle - << data_out.default_suffix (output_format) - << std::ends; -#ifdef HAVE_STD_STRINGSTREAM - std::ofstream out (filename.str().c_str()); -#else - std::ofstream out (filename.str()); -#endif - - data_out.write (out, output_format); - }; - - - - }; @@ -224,9 +167,15 @@ namespace LaplaceSolver 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; protected: const SmartPointer > triangulation; + + unsigned int refinement_cycle; }; @@ -240,6 +189,15 @@ namespace LaplaceSolver template Base::~Base () {}; + + + + template + void + Base::set_refinement_cycle (const unsigned int cycle) + { + refinement_cycle = cycle; + }; @@ -250,6 +208,7 @@ namespace LaplaceSolver Solver (Triangulation &triangulation, const FiniteElement &fe, const Quadrature &quadrature, + const Quadrature &face_quadrature, const Function &boundary_values); virtual ~Solver (); @@ -269,6 +228,7 @@ namespace LaplaceSolver protected: const SmartPointer > fe; const SmartPointer > quadrature; + const SmartPointer > face_quadrature; DoFHandler dof_handler; Vector solution; const SmartPointer > boundary_values; @@ -304,11 +264,13 @@ namespace LaplaceSolver 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) {}; @@ -506,6 +468,7 @@ namespace LaplaceSolver PrimalSolver (Triangulation &triangulation, const FiniteElement &fe, const Quadrature &quadrature, + const Quadrature &face_quadrature, const Function &rhs_function, const Function &boundary_values); @@ -533,12 +496,14 @@ namespace LaplaceSolver 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, boundary_values), + quadrature, face_quadrature, + boundary_values), rhs_function (&rhs_function) {}; @@ -664,6 +629,7 @@ namespace LaplaceSolver RefinementKelly (Triangulation &coarse_grid, const FiniteElement &fe, const Quadrature &quadrature, + const Quadrature &face_quadrature, const Function &rhs_function, const Function &boundary_values); @@ -677,11 +643,13 @@ namespace LaplaceSolver 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) {}; @@ -722,11 +690,15 @@ double Solution::value (const Point &p, const unsigned int /*component*/) const { - double q = p(0); - for (unsigned int i=1; i::value (const Point &p, const unsigned int /*component*/) const { - double q = p(0); - for (unsigned int i=1; i &dof_handler, Vector &rhs) const; + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + protected: const Point evaluation_point; }; @@ -846,6 +827,7 @@ namespace LaplaceSolver DualSolver (Triangulation &triangulation, const FiniteElement &fe, const Quadrature &quadrature, + const Quadrature &face_quadrature, const DualFunctional::DualFunctionalBase &dual_functional); // XXX @@ -868,17 +850,21 @@ namespace LaplaceSolver static const ZeroFunction boundary_values; }; + 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, boundary_values), + quadrature, face_quadrature, + boundary_values), dual_functional (&dual_functional) {}; @@ -925,8 +911,10 @@ namespace LaplaceSolver { public: WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &fe, + 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); @@ -943,25 +931,248 @@ namespace LaplaceSolver unsigned int n_dofs () const; - virtual void refine_grid () {}; + virtual void refine_grid (); + + virtual + void + output_solution () const; + + private: + + /** + * Declare a data type to + * represent the mapping between + * faces and integrated jumps of + * gradients of each of the + * solution vectors. See the + * general documentation of this + * class for more information. + */ + typedef typename std::pair::active_cell_iterator> FaceEntry; + typedef typename std::map::face_iterator,FaceEntry> FaceIntegrals; + + + /** + * Redeclare an active cell iterator. + * This is simply for convenience. + */ + typedef typename DoFHandler::active_cell_iterator active_cell_iterator; + + /** + * All data needed by the several + * functions of the error + * estimator is gathered in this + * struct. It is passed as a + * reference to the separate + * functions in this class. + * + * The reason for invention of + * this object is two-fold: + * first, class member data is + * not possible because no real + * object is created (all + * functions are @p{static}), + * which is a historical + * reason. Second, if we don't + * collect the data the various + * functions need somewhere at a + * central place, that would mean + * that the functions would have + * to allocate them upon + * need. However, then some + * variables would be allocated + * over and over again, which can + * take a significant amount of + * time (10-20 per cent) and most + * importantly, memory allocation + * requires synchronisation in + * multithreaded mode. While that + * is done by the C++ library and + * has not to be handcoded, it + * nevertheless seriously damages + * the ability to efficiently run + * the functions of this class in + * parallel, since they are quite + * often blocked by these + * synchronisation points. + * + * Thus, every thread gets an + * instance of this class to work + * with and needs not allocate + * memory itself, or synchronise + * with other threads. + */ + struct FaceData + { + 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 &dof_handler, + const Quadrature &face_quadrature); + }; + + struct CellData + { + 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 &dof_handler, + const Quadrature &quadrature, + const Function &right_hand_side); + }; + + + + 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; + + 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; + + /** + * Actually do the computation on + * a face which has no hanging + * nodes (it is regular), i.e. + * either on the other side there + * is nirvana (face is at + * boundary), or the other side's + * refinement level is the same + * as that of this side, then + * handle the integration of + * these both cases together. + * + * The meaning of the parameters + * becomes clear when looking at + * the source code. This function + * is only externalized from + * @p{estimate_error} to avoid + * ending up with a function of + * 500 lines of code. + */ + 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; + + + /** + * The same applies as for the + * function above, except that + * integration is over face + * @p{face_no} of @p{cell}, where + * the respective neighbor is + * refined, so that the + * integration is a bit more + * complex. + */ + 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; }; + + + template + WeightedResidual::FaceData:: + FaceData (const FiniteElement &fe, + const Quadrature &face_quadrature) + : + fe_face_values_cell (fe, face_quadrature, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors), + fe_face_values_neighbor (fe, face_quadrature, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors), + fe_subface_values_cell (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::CellData:: + CellData (const FiniteElement &fe, + const Quadrature &quadrature, + const Function &right_hand_side) + : + fe_values (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:: WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &fe, + 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) : Base (coarse_grid), - PrimalSolver (coarse_grid, fe, quadrature, + PrimalSolver (coarse_grid, primal_fe, + quadrature, face_quadrature, rhs_function, boundary_values), - DualSolver (coarse_grid, fe, quadrature, - dual_functional) + DualSolver (coarse_grid, dual_fe, + quadrature, face_quadrature, + dual_functional) {}; @@ -988,10 +1199,688 @@ namespace LaplaceSolver { return PrimalSolver::n_dofs(); }; + + + + template + void + WeightedResidual::refine_grid () + { + Vector error_indicators (triangulation->n_active_cells()); + estimate_error (error_indicators); + std::cout << "Estimated error=" + << std::accumulate (error_indicators.begin(), + error_indicators.end(), 0.) + << std::endl; + DataOut data_out; + ofstream x("x"); + Vector xe (error_indicators.begin(), + error_indicators.end()); + data_out.attach_dof_handler (DualSolver::dof_handler); + data_out.add_data_vector (xe, "e"); + data_out.build_patches (); + data_out.write_gnuplot (x); + + std::transform (error_indicators.begin(), + error_indicators.end(), + error_indicators.begin(), + &fabs); + GridRefinement::refine_and_coarsen_fixed_number (*triangulation, + error_indicators, + 0.3, 0.03); + triangulation->execute_coarsening_and_refinement (); + }; + + + template + void + WeightedResidual::output_solution () const + { + for (unsigned int solution_type=0; solution_type<2; ++solution_type) + { + DataOut data_out; + + switch (solution_type) + { + case 0: + data_out.attach_dof_handler (PrimalSolver::dof_handler); + data_out.add_data_vector (PrimalSolver::solution, + "primal_solution"); + break; + case 1: + data_out.attach_dof_handler (DualSolver::dof_handler); + data_out.add_data_vector (DualSolver::solution, + "dual_solution"); + break; + default: + Assert (false, ExcInternalError()); + }; + data_out.build_patches (); + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream filename; +#else + std::ostrstream filename; +#endif + filename << "solution-" + << (solution_type == 0 ? + "primal-" : "dual-") + << refinement_cycle + << ".gnuplot" + << std::ends; +#ifdef HAVE_STD_STRINGSTREAM + std::ofstream out (filename.str().c_str()); +#else + std::ofstream out (filename.str()); +#endif + + data_out.write (out, DataOut::gnuplot); + }; + }; + + + template + void + WeightedResidual:: + estimate_error (Vector &error_indicators) const + { + // 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 + // these two actions. + Vector primal_solution (DualSolver::dof_handler.n_dofs()); + FETools::interpolate (PrimalSolver::dof_handler, + PrimalSolver::solution, + DualSolver::dof_handler, + primal_solution); + Vector dual_weights (DualSolver::dof_handler.n_dofs()); + FETools::interpolation_difference (DualSolver::dof_handler, + DualSolver::solution, + *PrimalSolver::fe, + dual_weights); + + + // Map of integrals indexed by + // the corresponding face. In this map + // we store the integrated jump of the + // gradient for each face. + // At the end of the function, we again + // loop over the cells and collect the + // contributions of the different faces + // of the cell. + // + // The initial values for all faces + // are set to -1e20. It would cost + // a lot of time to synchronise the + // initialisation (i.e. the + // creation of new keys) of the map + // in multithreaded mode. Negative + // value indicates that the face + // has not yet been processed. + FaceIntegrals face_integrals; + for (active_cell_iterator cell=DualSolver::dof_handler.begin_active(); + cell!=DualSolver::dof_handler.end(); + ++cell) + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + face_integrals[cell->face(face_no)].first = -1e20; + + // reserve one slot for each cell + // and set it to zero + error_indicators.reinit (DualSolver::dof_handler + .get_tria().n_active_cells()); + + + // all the data needed in the error + // estimator by each of the threads + // is gathered in the following + // stuctures + // + // note that if no component mask + // was given, then treat all + // components + const unsigned int n_threads = multithread_info.n_default_threads; + + // split all cells into threads if + // multithreading is used and run + // the whole thing + Threads::ThreadManager thread_manager; + for (unsigned int i=0; i:: + estimate_some) + .collect_args (this, + primal_solution, + dual_weights, + n_threads, i, + error_indicators, + face_integrals)); + thread_manager.wait(); + + // finally add up the + // contributions of the faces for + // each cell + + unsigned int present_cell=0; + for (active_cell_iterator cell=DualSolver::dof_handler.begin_active(); + cell!=DualSolver::dof_handler.end(); + ++cell, ++present_cell) + { + // loop over all faces of this 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()); + if (true || (face_integrals[cell->face(face_no)].second + == + cell)) + error_indicators(present_cell) + += -0.5*face_integrals[cell->face(face_no)].first; + else + error_indicators(present_cell) + -= -0.5*face_integrals[cell->face(face_no)].first; + }; + }; + }; + + + 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 + { + FaceData face_data (*DualSolver::fe, + *DualSolver::face_quadrature); + CellData cell_data (*DualSolver::fe, + *DualSolver::quadrature, + *PrimalSolver::rhs_function); + + // First 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=DualSolver::dof_handler.begin_active(); + for (unsigned int t=0; + (t::dof_handler.end()); + ++t, ++cell); + + + // Then loop over all cells. The + // check for loop end is done at + // the end of the loop, along + // with incrementing the loop + // index. + for (unsigned int cell_index=this_thread; true; ) + { + + integrate_over_cell (cell, cell_index, + primal_solution, + dual_weights, + cell_data, + error_indicators); + + // After computing the cell + // terms, turn to the face + // terms. For this, loop over + // all faces of the present + // cell, and see whether + // something needs to be + // computed on it: + for (unsigned int face_no=0; + face_no::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, and also + // mark the cell on which + // we computed this + // value. + if (cell->face(face_no)->at_boundary()) + { + face_integrals[cell->face(face_no)].first = 0; + face_integrals[cell->face(face_no)].second = cell; + 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::dof_handler.end())); + ++t, ++cell, ++cell_index); + if (cell == DualSolver::dof_handler.end()) + break; + }; + }; + + + template + 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 + { + 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); + cell_data.fe_values.get_function_values (dual_weights, + cell_data.dual_weights); + 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)].first == -1e20, + ExcInternalError()); + + // ...then store computed value + // at assigned location: + face_integrals[cell->face(face_no)].first = face_integral; + face_integrals[cell->face(face_no)].second = cell; + }; + + + + 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); - template class WeightedResidual<2>; + // 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) + { + // Start with some checks + // again: get an iterator + // pointing to the cell + // behind the present subface + // and check whether its face + // is a subface of the one we + // are considering. If that + // were not the case, then + // there would be either a + // bug in the + // ``neighbor_neighbor'' + // function called above, or + // -- worse -- some function + // in the library did not + // keep to some underlying + // assumptions about cells, + // their children, and their + // faces. In any case, even + // though this assertion + // should not be triggered, + // it does not harm to be + // cautious, and in optimized + // mode computations the + // assertion will be removed + // anyway. + const active_cell_iterator neighbor_child + = neighbor->child(GeometryInfo:: + child_cell_on_face(neighbor_neighbor, + subface_no)); + Assert (neighbor_child->face(neighbor_neighbor) == + cell->face(face_no)->child(subface_no), + ExcInternalError()); + + // 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)].first + = face_integral; + face_integrals[neighbor_child->face(neighbor_neighbor)].second + = cell; + }; + // 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)].first != -1e20, + ExcInternalError()); + + sum += face_integrals[face->child(subface_no)].first; + }; + // Finally store the value with + // the parent face. + face_integrals[face].first = sum; + face_integrals[face].second = cell; + }; + }; @@ -1006,9 +1895,13 @@ run_simulation (LaplaceSolver::Base &solver, for (unsigned int step=0; true; ++step) { - std::cout << step << " " << std::flush; + std::cout << step << " Solving " + << solver.n_dofs() + << std::endl; + solver.set_refinement_cycle (step); solver.solve_problem (); + solver.output_solution (); for (typename std::list *>::const_iterator i = postprocessor_list.begin(); @@ -1019,7 +1912,7 @@ run_simulation (LaplaceSolver::Base &solver, }; - if (solver.n_dofs() < 20000) + if (solver.n_dofs() < 5000) solver.refine_grid (); else break; @@ -1030,6 +1923,52 @@ run_simulation (LaplaceSolver::Base &solver, +void +create_triangulation (Triangulation<2> &tria) +{ + const Point<2> + vertices[16] = { Point<2> (-1., -1.), + Point<2> (-1./3, -1.), + Point<2> (+1./3, -1.), + Point<2> (+1, -1.), + Point<2> (-1., -1./3.), + Point<2> (-1./3, -1./3.), + Point<2> (+1./3, -1./3.), + Point<2> (+1, -1./3.), + Point<2> (-1., 1./3.), + Point<2> (-1./3, 1./3.), + Point<2> (+1./3, 1./3.), + Point<2> (+1, 1./3.), + Point<2> (-1., 1.), + Point<2> (-1./3, 1.), + Point<2> (+1./3, 1.), + Point<2> (+1, 1.) }; + + const int cell_vertices[8][4] = {{0, 1, 5, 4}, + {1, 2, 6, 5}, + {2, 3, 7, 6}, + {4, 5, 9, 8}, + {6, 7, 11, 10}, + {8,9,13,12}, + {9,10,14,13}, + {10,11,15,14}}; + + std::vector > cells (8, CellData<2>()); + + for (unsigned int i=0; i<8; ++i) + { + for (unsigned int j=0; j<4; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; + + tria.create_triangulation (std::vector >(&vertices[0], &vertices[16]), + cells, + SubCellData()); // no boundary information +}; + + + template void solve_problem (const std::string &solver_name) { @@ -1039,38 +1978,36 @@ void solve_problem (const std::string &solver_name) << std::string (header.size(), '-') << std::endl; Triangulation triangulation; +// create_triangulation (triangulation); GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (2); - const FE_Q fe(1); + const FE_Q primal_fe(1); + const FE_Q dual_fe(2); const QGauss4 quadrature; + const QGauss4 face_quadrature; const RightHandSide rhs_function; const Solution boundary_values; + const Point evaluation_point(0.5,0.5); + const DualFunctional::PointValueEvaluation + dual_functional (evaluation_point); + LaplaceSolver::Base * solver = 0; - if (solver_name == "global") - solver = new LaplaceSolver::RefinementGlobal (triangulation, fe, - quadrature, - rhs_function, - boundary_values); - else if (solver_name == "kelly") - solver = new LaplaceSolver::RefinementKelly (triangulation, fe, - quadrature, - rhs_function, - boundary_values); - else - AssertThrow (false, ExcNotImplemented()); + solver = new LaplaceSolver::WeightedResidual (triangulation, + primal_fe, + dual_fe, + quadrature, + face_quadrature, + rhs_function, + boundary_values, + dual_functional); TableHandler results_table; Evaluation::PointValueEvaluation postprocessor1 (Point(0.5,0.5), results_table); - Evaluation::SolutionOutput - postprocessor2 (std::string("solution-")+solver_name, - DataOut::gnuplot); - std::list *> postprocessor_list; postprocessor_list.push_back (&postprocessor1); - postprocessor_list.push_back (&postprocessor2); run_simulation (*solver, postprocessor_list); @@ -1089,7 +2026,7 @@ int main () deallog.depth_console (0); solve_problem<2> ("global"); - solve_problem<2> ("kelly"); +// solve_problem<2> ("kelly"); } catch (std::exception &exc) { -- 2.39.5