From 628442c30d4e53fe4c5237d73df4893637d1a4ba Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 11 Sep 2000 14:00:35 +0000 Subject: [PATCH] Implement estimating the error for several solution vectors at once. git-svn-id: https://svn.dealii.org/trunk@3318 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/error_estimator.h | 146 +++++-- .../source/numerics/error_estimator.cc | 397 +++++++++++------- 2 files changed, 367 insertions(+), 176 deletions(-) diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index e17a9e4b8d..3f0e5a6dd9 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -177,8 +177,30 @@ template class FESubfaceValues; * cells at hand and need not think about explicitely determining whether * a face was refined or not. The same applies for boundary faces, see * above. - * - * @author Wolfgang Bangerth, 1998, 1999; parallelization by Thomas Richter, 2000 + * + * + * @sect3{Multiple solution vectors} + * + * In some cases, for example in time-dependent problems, one would + * like to compute the error estimates for several solution vectors + * on the same grid at once, with the same coefficients, etc, + * e.g. for the solutions on several successive time steps. One could + * then call the functions of this class several times for each + * solution. However, it is observed that the largest factor in the + * computation of the error estimates (in terms of computing time) is + * initialization of @ref{FEFaceValues} and @ref{FESubFaceValues} + * objects, and iterating through all faces and subfaces. If the + * solution vectors live on the same grid, this effort can be reduced + * significantly by treating all solution vectors at the same time, + * initializing the @ref{FEFaceValues} objects only once per cell and + * for all solution vectors at once, and also only looping through + * the triangulation only once. For this reason, besides the + * @p{estimate} function in this class that takes a single input + * vector and returns a single output vector, there is also a + * function that accepts several in- and output vectors at the same + * time. + * + * @author Wolfgang Bangerth, 1998, 1999, 2000; parallelization by Thomas Richter, 2000 */ template class KellyErrorEstimator @@ -254,7 +276,44 @@ class KellyErrorEstimator Vector &error, const vector &component_mask = vector(), const Function *coefficients = 0, - unsigned int n_threads = multithread_info.n_default_threads); + unsigned int n_threads = multithread_info.n_default_threads); + + /** + * Same function as above, but + * accepts more than one solution + * vectors and returns one error + * vector for each solution + * vector. For the reason of + * existence of this function, + * see the general documentation + * of this class. + * + * Since we do not want to force + * the user of this function to + * copy around their solution + * vectors, the vector of + * solution vectors takes + * pointers to the solutions, + * rather than being a vector of + * vectors. This makes it simpler + * to have the solution vectors + * somewhere in memory, rather + * than to have them collected + * somewhere special. (Note that + * it is not possible to + * construct of vector of + * references, so we had to use a + * vector of pointers.) + */ + static void estimate (const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const vector*> &solutions, + vector*> &errors, + const vector &component_mask = vector(), + const Function *coefficients = 0, + unsigned int n_threads = multithread_info.n_default_threads); + /** * Exception @@ -272,7 +331,21 @@ class KellyErrorEstimator * Exception */ DeclException0 (ExcInvalidBoundaryFunction); - + /** + * Exception + */ + DeclException2 (ExcIncompatibleNumberOfElements, + int, int, + << "The number of elements " << arg1 << " and " << arg2 + << " of the vectors do not match!"); + /** + * Exception + */ + DeclException0 (ExcInvalidSolutionVector); + /** + * Exception + */ + DeclException0 (ExcNoSolutions); private: @@ -281,11 +354,12 @@ class KellyErrorEstimator * Declare a data type to * represent the mapping between * faces and integrated jumps of - * gradients. See the general - * documentation of this class - * for more information. + * gradients of each of the + * solution vectors. See the + * general documentation of this + * class for more information. */ - typedef map::face_iterator,double> FaceIntegrals; + typedef map::face_iterator,vector > FaceIntegrals; /** @@ -341,13 +415,14 @@ class KellyErrorEstimator */ struct Data { - const DoFHandler &dof_handler; - const Quadrature &quadrature; - const FunctionMap &neumann_bc; - const Vector &solution; - const vector component_mask; - const Function *coefficients; - const unsigned int n_threads; + const DoFHandler &dof_handler; + const Quadrature &quadrature; + const FunctionMap &neumann_bc; + const vector*> &solutions; + const vector component_mask; + const Function *coefficients; + const unsigned int n_threads; + const unsigned int n_solution_vectors; /** * Reference to the global @@ -359,8 +434,9 @@ class KellyErrorEstimator /** * A vector to store the jump * of the normal vectors in - * the quadrature points - * (i.e. a temporary + * the quadrature points for + * each of the solution + * vectors (i.e. a temporary * value). This vector is not * allocated inside the * functions that use it, but @@ -371,26 +447,30 @@ class KellyErrorEstimator * synchronisation makes * things even slower. */ - vector > phi; + vector > > phi; /** * A vector for the gradients of * the finite element function * on one cell * - * Let psi be a short name for - * @p{a grad u_h}, where the second - * index be the component of the - * finite element, and the first + * Let psi be a short name + * for @p{a grad u_h}, where + * the third index be the + * component of the finite + * element, and the second * index the number of the - * quadrature point. + * quadrature point. The + * first index denotes the + * index of the solution + * vector. */ - vector > > psi; + vector > > > psi; /** * The same vector for a neighbor cell */ - vector > > neighbor_psi; + vector > > > neighbor_psi; /** * The normal vectors of the finite @@ -420,14 +500,14 @@ class KellyErrorEstimator * class Data. All variables are * passed as references. */ - Data(const DoFHandler &dof, - const Quadrature &quadrature, - const FunctionMap &neumann_bc, - const Vector &solution, - const vector &component_mask, - const Function *coefficients, - const unsigned int n_threads, - FaceIntegrals &face_integrals); + Data(const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const vector*> &solutions, + const vector &component_mask, + const Function *coefficients, + const unsigned int n_threads, + FaceIntegrals &face_integrals); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index a117fc16d1..0eb8b29e8f 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -46,18 +46,18 @@ double sqr (const double x) #if deal_II_dimension == 1 template <> -KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, - const Quadrature<0> &, - const FunctionMap &, - const Vector &, - const vector &, - const Function<1> *, - const unsigned int , - FaceIntegrals &): +KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, + const Quadrature<0> &, + const FunctionMap &, + const vector*> &, + const vector &, + const Function<1> *, + const unsigned int , + FaceIntegrals &): dof_handler(* static_cast *> (0)), quadrature(* static_cast *> (0)), neumann_bc(* static_cast (0)), - solution(* static_cast *> (0)), + solutions(* static_cast *> *> (0)), face_integrals (* static_cast (0)) { Assert (false, ExcInternalError()); @@ -66,21 +66,22 @@ KellyErrorEstimator<1>::Data::Data(const DoFHandler<1> &, #else template -KellyErrorEstimator::Data::Data(const DoFHandler &dof_handler, - const Quadrature &quadrature, - const FunctionMap &neumann_bc, - const Vector &solution, - const vector &component_mask, - const Function *coefficients, - const unsigned int n_threads, - FaceIntegrals &face_integrals): +KellyErrorEstimator::Data::Data(const DoFHandler &dof_handler, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const vector*> &solutions, + const vector &component_mask, + const Function *coefficients, + const unsigned int n_threads, + FaceIntegrals &face_integrals): dof_handler (dof_handler), quadrature (quadrature), neumann_bc (neumann_bc), - solution (solution), + solutions (solutions), component_mask (component_mask), coefficients (coefficients), n_threads (n_threads), + n_solution_vectors (solutions.size()), face_integrals (face_integrals) { const unsigned int n_components = dof_handler.get_fe().n_components(); @@ -105,27 +106,60 @@ KellyErrorEstimator::Data::Data(const DoFHandler &dof_handler, // needed in the calculations once // per thread. const unsigned int n_q_points = quadrature.n_quadrature_points; - phi.resize(n_q_points); - psi.resize(n_q_points); - neighbor_psi.resize(n_q_points); + normal_vectors.resize(n_q_points); coefficient_values1.resize(n_q_points); coefficient_values.resize(n_q_points); - JxW_values.resize(n_q_points); - - for (unsigned int qp=0;qp +void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const Vector &solution, + Vector &error, + const vector &component_mask, + const Function *coefficients, + unsigned int n_threads) +{ + // just pass on to the other function + const vector*> solutions (1, &solution); + vector*> errors (1, &error); + estimate (dof_handler, quadrature, neumann_bc, solutions, errors, + component_mask, coefficients, n_threads); +}; + + + + #if deal_II_dimension == 1 template <> @@ -139,17 +173,27 @@ void KellyErrorEstimator<1>::estimate_some (Data &, const unsigned int) template <> -void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler, - const Quadrature<0> &, - const FunctionMap &neumann_bc, - const Vector &solution, - Vector &error, - const vector &component_mask_, - const Function<1> *coefficient, +void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler, + const Quadrature<0> &, + const FunctionMap &neumann_bc, + const vector*> &solutions, + vector*> &errors, + const vector &component_mask_, + const Function<1> *coefficient, const unsigned int) { - const unsigned int n_components = dof_handler.get_fe().n_components(); - + const unsigned int n_components = dof_handler.get_fe().n_components(); + const unsigned int n_solution_vectors = solutions.size(); + + // sanity checks + Assert (solutions.size() > 0, + ExcNoSolutions()); + Assert (solutions.size() == errors.size(), + ExcIncompatibleNumberOfElements(solutions.size(), errors.size())); + for (unsigned int n=0; nsize() == dof_handler.n_dofs(), + ExcInvalidSolutionVector()); + // if no mask given: treat all components vector component_mask ((component_mask_.size() == 0) ? vector(n_components, true) : @@ -171,7 +215,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler, // reserve one slot for each cell and set // it to zero - error.reinit (dof_handler.get_tria().n_active_cells()); + for (unsigned int n=0; n::estimate (const DoFHandler<1> &dof_handler, // need several auxiliary fields, // depending on the way we get it // (see below) - vector > > gradients_here (2, vector >(n_components)); - vector > > gradients_neighbor (gradients_here); - Vector grad_neighbor (n_components); + vector > > > + gradients_here (n_solution_vectors, + vector > >(2, vector >(n_components))); + vector > > > + gradients_neighbor (gradients_here); + vector > + grad_neighbor (n_solution_vectors, Vector(n_components)); // reserve some space for // coefficient values at one point. @@ -204,7 +253,9 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_handler, active_cell_iterator cell = dof_handler.begin_active(); for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index) { - error(cell_index) = 0; + for (unsigned int n=0; n::estimate (const DoFHandler<1> &dof_handler, // now get the gradients on the // both sides of the point fe_values.reinit (cell); - fe_values.get_function_grads (solution, gradients_here); + + for (unsigned int s=0; s::estimate (const DoFHandler<1> &dof_handler, // means: x-derivative, // which is the only // one here - for (unsigned int c=0; csecond->vector_value(cell->vertex(0), - grad_neighbor); + for (unsigned int s=0; ssecond->vector_value(cell->vertex(0), + grad_neighbor[s]); else // fill with zeroes. - grad_neighbor.clear (); + for (unsigned int s=0; s::estimate (const DoFHandler<1> &dof_handler, }; - for (unsigned int component=0; componentdiameter(); - }; + for (unsigned int s=0; sdiameter(); + }; }; - error(cell_index) = sqrt(error(cell_index)); + for (unsigned int s=0; s void KellyErrorEstimator::estimate_some (Data &data, const unsigned int this_thread) { + const unsigned int n_solution_vectors = data.n_solution_vectors; // make up a fe face values object for the // restriction of the finite element function @@ -353,9 +416,14 @@ void KellyErrorEstimator::estimate_some (Data &data, // loop over all faces of this cell for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - // if we already visited this - // face: do nothing - if (data.face_integrals[cell->face(face_no)] >=0) + // if we already visited + // this face: do + // nothing. only check + // component for first + // solution vector, as we + // treat them all at the + // same time + if (data.face_integrals[cell->face(face_no)][0] >=0) continue; @@ -381,7 +449,8 @@ void KellyErrorEstimator::estimate_some (Data &data, && data.neumann_bc.find(boundary_indicator)==data.neumann_bc.end()) { - data.face_integrals[cell->face(face_no)] = 0; + for (unsigned int n=0; nface(face_no)][n] = 0; continue; }; @@ -423,21 +492,33 @@ void KellyErrorEstimator::estimate_some (Data &data, template -void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, - const Quadrature &quadrature, - const FunctionMap &neumann_bc, - const Vector &solution, - Vector &error, - const vector &component_mask, - const Function *coefficients, - unsigned int n_threads) +void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const vector*> &solutions, + vector*> &errors, + const vector &component_mask, + const Function *coefficients, + unsigned int n_threads) { + // sanity checks + Assert (solutions.size() > 0, + ExcNoSolutions()); + Assert (solutions.size() == errors.size(), + ExcIncompatibleNumberOfElements(solutions.size(), errors.size())); + for (unsigned int n=0; nsize() == dof_handler.n_dofs(), + ExcInvalidSolutionVector()); + + // if NOT multithreaded, set n_threads to one #ifndef DEAL_II_USE_MT n_threads = 1; #endif Assert (n_threads > 0, ExcInternalError()); + const unsigned int n_solution_vectors = solutions.size(); + // Map of integrals indexed by // the corresponding face. In this map // we store the integrated jump of the @@ -455,10 +536,11 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, // in multithreaded mode. Negative // value indicates that the face // has not yet been processed. + vector default_face_integrals (n_solution_vectors, -10e20); FaceIntegrals face_integrals; for (active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - face_integrals[cell->face(face_no)]=-10e20; + face_integrals[cell->face(face_no)] = default_face_integrals; // all the data needed in the error @@ -474,7 +556,7 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, data_structures[i] = new Data (dof_handler, quadrature, neumann_bc, - solution, + solutions, ((component_mask.size() == 0) ? vector(dof_handler.get_fe().n_components(), true) : @@ -509,9 +591,12 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, // reserve one slot for each cell and set // it to zero - error.reinit (dof_handler.get_tria().n_active_cells()); - for (unsigned int i=0;i::estimate (const DoFHandler &dof_handler, { Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(), ExcInternalError()); - error(present_cell) += (face_integrals[cell->face(face_no)] * - cell->diameter() / 24); + const double factor = cell->diameter() / 24; + + for (unsigned int n=0; nface(face_no)][n] * + factor); }; - error(present_cell) = sqrt(error(present_cell)); + + for (unsigned int n=0; n &fe_face_values_neighbor) { const DoFHandler::face_iterator face = cell->face(face_no); - const unsigned int n_q_points = data.quadrature.n_quadrature_points, - n_components = data.dof_handler.get_fe().n_components(); + const unsigned int n_q_points = data.quadrature.n_quadrature_points, + n_components = data.dof_handler.get_fe().n_components(), + n_solution_vectors = data.n_solution_vectors; // initialize data of the restriction @@ -581,7 +672,8 @@ integrate_over_regular_face (Data &data, // get gradients of the finite element // function on this cell - fe_face_values_cell.get_function_grads (data.solution, data.psi); + for (unsigned int n=0; nvalue_list (fe_face_values_cell.get_quadrature_points(), data.coefficient_values1); + for (unsigned int n=0; nvector_value_list (fe_face_values_cell.get_quadrature_points(), - data.coefficient_values); - for (unsigned int component=0; componentvector_value_list (fe_face_values_cell.get_quadrature_points(), g); - for (unsigned int component=0; component face_integral (n_solution_vectors, 0); + for (unsigned int n=0; n &fe_subface_values) { const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); - const unsigned int n_q_points = data.quadrature.n_quadrature_points, - n_components = data.dof_handler.get_fe().n_components(); + const unsigned int n_q_points = data.quadrature.n_quadrature_points, + n_components = data.dof_handler.get_fe().n_components(), + n_solution_vectors = data.n_solution_vectors; Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); @@ -768,20 +869,25 @@ integrate_over_irregular_face (Data &data, // store the gradient of the solution // in psi fe_subface_values.reinit (cell, face_no, subface_no); - fe_subface_values.get_function_grads (data.solution, data.psi); + + for (unsigned int n=0; nvalue_list (fe_face_values.get_quadrature_points(), data.coefficient_values1); - for (unsigned int component=0; componentvector_value_list (fe_face_values.get_quadrature_points(), data.coefficient_values); - for (unsigned int component=0; component face_integral (n_solution_vectors, 0); + for (unsigned int n=0; nface(neighbor_neighbor)] = face_integral; }; @@ -856,7 +965,7 @@ integrate_over_irregular_face (Data &data, // collect the contributions of the // subfaces and store them with the // mother face - double sum=0; + vector sum (n_solution_vectors, 0); DoFHandler::face_iterator face = cell->face(face_no); for (unsigned int subface_no=0; subface_no::subfaces_per_face; ++subface_no) @@ -864,9 +973,11 @@ integrate_over_irregular_face (Data &data, Assert (data.face_integrals.find(face->child(subface_no)) != data.face_integrals.end(), ExcInternalError()); - Assert (data.face_integrals[face->child(subface_no)]>=0, + Assert (data.face_integrals[face->child(subface_no)][0] >= 0, ExcInternalError()); - sum += data.face_integrals[face->child(subface_no)]; + + for (unsigned int n=0; nchild(subface_no)][n]; }; data.face_integrals[face] = sum; -- 2.39.5