From: wolf Date: Fri, 29 Oct 1999 09:45:43 +0000 (+0000) Subject: Implement the error estimator for several components at a time, with different weight... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6f7395c1dca0617e7932999aaba850b3852318bd;p=dealii-svn.git Implement the error estimator for several components at a time, with different weights and all other nifty features one needs. git-svn-id: https://svn.dealii.org/trunk@1806 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 85999c4569..fb802549d5 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -16,12 +16,12 @@ /** - * Implementation of the error estimator by Kelly, Gago, Zienkiewicz and - * Babuska. - * This error estimator tries to approximate the error per cell by integration - * of the jump of the gradient of the solution along the faces of each cell. - * It can be understood as a gradient recovery estimator; see the survey - * of Ainsworth for a complete discussion. + * Implementation of the error estimator by Kelly, Gago, Zienkiewicz + * and Babuska. This error estimator tries to approximate the error + * per cell by integration of the jump of the gradient of the + * solution along the faces of each cell. It can be understood as a + * gradient recovery estimator; see the survey of Ainsworth for a + * complete discussion. * * It seem as if this error estimator should only be valid for linear trial * spaces, and there are indications that for higher order trial spaces the @@ -83,6 +83,28 @@ * the diameter of the cell. * * + * \subsection{Vector-valued functions} + * + * If the finite element field for which the error is to be estimated + * is vector-valued, i.e. the finite element has more than one + * component, then all the above can be applied to all or only some + * components at the same time. The main function of this class takes + * a list of flags denoting the components for which components the + * error estimator is to be applied; by default, it is a list of only + * #true#s, meaning that all variables shall be treated. + * + * In case the different components of a field have different + * physical meaning (for example velocity and pressure in the Stokes + * equations), it would be senseless to use the same coefficient for + * all components. In that case, you can pass a function with as many + * components as there are components in the finite element field and + * each component of the error estimator will then be weighted by the + * respective component in this coefficient function. In the other + * case, when all components have the same meaning (for example the + * displacements in Lame's equations of elasticity), you can specify + * a scalar coefficient which will then be used for all components. + * + * * \subsection{Boundary values} * * If the face is at the boundary, i.e. there is no neighboring cell to which @@ -108,7 +130,9 @@ * \item The face belongs to a Neumann boundary. In this case, the * contribution of the face $F\in\partial K$ looks like * $$ \int_F \left|g-a\frac{\partial u_h}{\partial n}\right|^2 ds $$ - * where $g$ is the Neumann boundary function. + * where $g$ is the Neumann boundary function. If the finite element is + * vector-valued, then obviously the function denoting the Neumann boundary + * conditions needs to be vector-valued as well. * * \item No other boundary conditions are considered. * \end{itemize} @@ -160,30 +184,50 @@ class KellyErrorEstimator { typedef map*> FunctionMap; /** - * Implementation of the error estimator - * described above. You may give a - * coefficient, but there is a default - * value which denotes the constant - * coefficient with value one. + * Implementation of the error + * estimator described above. You + * may give a coefficient, but + * there is a default value which + * denotes the constant + * coefficient with value + * one. The coefficient function + * may either be a scalar one, in + * which case it is used for all + * components of the finite + * element, or a vector-valued + * one with as many components as + * there are in the finite + * element; in the latter case, + * each component is weighted by + * the respective component in + * the coefficient. * - * You must give the component if the - * finite element in use by the #dof# - * object has more than one component. - * This number shall be between zero - * and the number of components within - * the finite element. If the finite - * element has only one component, - * then the parameter selecting the - * component shall be zero, which is - * also the default value. + * You might give a list of + * components you want to + * evaluate, in case the finite + * element used by the + * #DoFHandler# object is + * vector-valued. You then have + * to set those entries to true + * in the bit-vector + * #component_mask# for which the + * respective component is to be + * used in the error + * estimator. The default is to + * use all components, which is + * done by either providing a + * bit-vector with all-set + * entries, or an empty + * bit-vector. */ - static void estimate (const DoFHandler &dof, - const Quadrature &quadrature, - const FunctionMap &neumann_bc, - const Vector &solution, - Vector &error, - const Function *coefficient = 0, - const unsigned int selected_component = 0); + static void estimate (const DoFHandler &dof, + const Quadrature &quadrature, + const FunctionMap &neumann_bc, + const Vector &solution, + Vector &error, + const vector &component_mask = vector(), + const Function *coefficients = 0); + /** * Exception */ @@ -191,12 +235,16 @@ class KellyErrorEstimator { /** * Exception */ - DeclException2 (ExcInvalidComponent, - int, int, - << "The component you gave (" << arg1 << ") " - << "is invalid with respect to the number " - << "of components in the finite element " - << "(" << arg2 << ")"); + DeclException0 (ExcInvalidComponentMask); + /** + * Exception + */ + DeclException0 (ExcInvalidCoefficient); + /** + * Exception + */ + DeclException0 (ExcInvalidBoundaryFunction); + private: /** * Declare a data type to represent the @@ -238,8 +286,7 @@ class KellyErrorEstimator { FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, const Vector&solution, - const unsigned int n_components, - const unsigned int selected_component, + const vector &component_mask, const Function *coefficient); /** @@ -257,8 +304,7 @@ class KellyErrorEstimator { FESubfaceValues &fe_subface_values, FaceIntegrals &face_integrals, const Vector &solution, - const unsigned int n_components, - const unsigned int selected_component, + const vector &component_mask, const Function *coefficient); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 835c193af7..36010e7b46 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -12,12 +12,11 @@ #include #include #include -#include #include #include #include - +#include inline static double sqr (const double x) { @@ -34,13 +33,27 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const Function<1> *coefficient, - const unsigned int selected_component) + const vector &component_mask_, + const Function<1> *coefficient) { - Assert (selected_component < dof.get_fe().n_components, - ExcInvalidComponent (selected_component, dof.get_fe().n_components)); - Assert (coefficient->n_components == 1, - ExcInternalError()); + const unsigned int n_components = dof.get_fe().n_components; + + // if no mask given: treat all components + vector component_mask ((component_mask_.size() == 0) ? + vector(n_components, true) : + component_mask_); + Assert (component_mask.size() == n_components, ExcInvalidComponentMask()); + Assert (count(component_mask.begin(), component_mask.end(), true) > 0, + ExcInvalidComponentMask()); + + Assert ((coefficient == 0) || + (coefficient->n_components == n_components) || + (coefficient->n_components == 1), + ExcInvalidCoefficient()); + + for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i) + Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); + const unsigned int dim=1; @@ -48,16 +61,39 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, // it to zero error.reinit (dof.get_tria().n_active_cells()); - // loop over all cells. note that the - // error indicator is only a sum over - // the two contributions from the two + // fields to get the gradients on + // the present and the neighbor cell. + // + // for the neighbor gradient, we + // 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); + + // reserve some space for + // coefficient values at one point. + // if there is no coefficient, then + // we fill it by unity once and for + // all and don't set it any more + Vector coefficient_values (n_components); + if (coefficient == 0) + for (unsigned int c=0; c quadrature; - FEValues fe_values (dof.get_fe(), quadrature, update_gradients); - DoFHandler::active_cell_iterator cell = dof.begin_active(); + FEValues<1> fe_values (dof.get_fe(), quadrature, update_gradients); + DoFHandler<1>::active_cell_iterator cell = dof.begin_active(); for (unsigned int cell_index=0; cell != dof.end(); ++cell, ++cell_index) { - // loop over te two points bounding + error(cell_index) = 0; + // loop over the two points bounding // this line. n==0 is left point, // n==1 is right point for (unsigned int n=0; n<2; ++n) @@ -70,32 +106,70 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, // now get the gradients on the // both sides of the point - vector > > - gradients (2, vector >(dof.get_fe().n_components)); - fe_values.reinit (cell); - fe_values.get_function_grads (solution, gradients); - const double grad_here = gradients[n][selected_component][0]; + fe_values.get_function_grads (solution, gradients_here); - double grad_neighbor; if (neighbor.state() == valid) { fe_values.reinit (neighbor); - fe_values.get_function_grads (solution, gradients); - grad_neighbor = gradients[n==0 ? 1 : 0][selected_component][0]; + fe_values.get_function_grads (solution, gradients_neighbor); + + // extract the + // gradients of all the + // components. [0] + // means: x-derivative, + // which is the only + // one here + for (unsigned int c=0; csecond->value(cell->vertex(0)); + // if Neumann b.c., then fill + // the gradients field which + // will be used later on. + neumann_bc.find(n)->second->vector_value(cell->vertex(0), + grad_neighbor); else - grad_neighbor = 0; + // fill with zeroes. + grad_neighbor.clear (); + + // if there is a + // coefficient, then + // evaluate it at the + // present position. if + // there is none, reuse the + // preset values. + if (coefficient != 0) + { + if (coefficient->n_components == 1) + { + const double c_value = coefficient->value (cell->vertex(n)); + for (unsigned int c=0; cvector_value(cell->vertex(n), + coefficient_values); + }; + + + for (unsigned int component=0; componentvalue(cell->vertex(n)) : - 1); - error(cell_index) += jump*jump * cell->diameter(); + const double jump = ((grad_here - grad_neighbor(component)) * + coefficient_values(component)); + error(cell_index) += jump*jump * cell->diameter(); + }; }; + error(cell_index) = sqrt(error(cell_index)); }; }; @@ -103,19 +177,36 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, #endif + template void KellyErrorEstimator::estimate (const DoFHandler &dof, const Quadrature &quadrature, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const Function *coefficient, - const unsigned int selected_component) + const vector &component_mask_, + const Function *coefficient) { + const unsigned int n_components = dof.get_fe().n_components; + + // if no mask given: treat all components + vector component_mask ((component_mask_.size() == 0) ? + vector(n_components, true) : + component_mask_); + Assert (component_mask.size() == n_components, ExcInvalidComponentMask()); + Assert (count(component_mask.begin(), component_mask.end(), true) > 0, + ExcInvalidComponentMask()); + + Assert ((coefficient == 0) || + (coefficient->n_components == n_components) || + (coefficient->n_components == 1), + ExcInvalidCoefficient()); + Assert (neumann_bc.find(255) == neumann_bc.end(), ExcInvalidBoundaryIndicator()); - Assert (selected_component < dof.get_fe().n_components, - ExcInvalidComponent (selected_component, dof.get_fe().n_components)); + + for (FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i) + Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); // create a map of integrals indexed by // the corresponding face. In this map @@ -206,8 +297,7 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, fe_face_values_neighbor, face_integrals, solution, - dof.get_fe().n_components, - selected_component, + component_mask, coefficient); else // otherwise we need to do some @@ -219,8 +309,7 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof, fe_face_values_cell, fe_subface_values, face_integrals, solution, - dof.get_fe().n_components, - selected_component, + component_mask, coefficient); }; @@ -262,8 +351,7 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iter FEFaceValues<1> &, FaceIntegrals &, const Vector &, - const unsigned int , - const unsigned int , + const vector &, const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -279,8 +367,7 @@ integrate_over_irregular_face (const active_cell_iterator &, FESubfaceValues<1> &, FaceIntegrals &, const Vector &, - const unsigned int , - const unsigned int , + const vector &, const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -299,9 +386,10 @@ integrate_over_regular_face (const active_cell_iterator &cell, FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, const Vector &solution, - const unsigned int n_components, - const unsigned int selected_component, - const Function *coefficient) { + const vector &component_mask, + const Function *coefficient) +{ + const unsigned int n_components = component_mask.size(); const DoFHandler::face_iterator face = cell->face(face_no); // initialize data of the restriction @@ -314,20 +402,13 @@ integrate_over_regular_face (const active_cell_iterator &cell, // points // // let psi be a short name for - // [a grad u_h] - vector > psi(n_q_points); - if (n_components == 1) - fe_face_values_cell.get_function_grads (solution, psi); - else - { - vector > > tmp (n_q_points, - vector >(n_components)); - fe_face_values_cell.get_function_grads (solution, tmp); - for (unsigned int i=0; i > > psi(n_q_points, vector >(n_components)); + fe_face_values_cell.get_function_grads (solution, psi); // now compute over the other side of // the face @@ -338,8 +419,7 @@ integrate_over_regular_face (const active_cell_iterator &cell, Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); unsigned int neighbor_neighbor; - DoFHandler::active_cell_iterator neighbor - = cell->neighbor(face_no); + DoFHandler::active_cell_iterator neighbor = cell->neighbor(face_no); // find which number the current // face has relative to the neighboring @@ -359,29 +439,18 @@ integrate_over_regular_face (const active_cell_iterator &cell, // get a list of the gradients of // the finite element solution // restricted to the neighbor cell - vector > neighbor_psi (n_q_points); - if (n_components == 1) - fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); - else - { - vector > > tmp (n_q_points, - vector >(n_components)); - fe_face_values_neighbor.get_function_grads (solution, tmp); - for (unsigned int i=0; i > > neighbor_psi (n_q_points, + vector >(n_components)); + fe_face_values_neighbor.get_function_grads (solution, neighbor_psi); // compute the jump in the gradients - transform (psi.begin(), psi.end(), - neighbor_psi.begin(), - psi.begin(), - minus >()); + for (unsigned int component=0; component phi(n_q_points,0); + vector > phi(n_q_points, vector(n_components)); const vector > &normal_vectors(fe_face_values_cell. get_normal_vectors()); - - for (unsigned int point=0; point coefficient_values (n_q_points); - coefficient->value_list (fe_face_values_cell.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointn_components == 1) + { + vector coefficient_values (n_q_points); + coefficient->value_list (fe_face_values_cell.get_quadrature_points(), + coefficient_values); + for (unsigned int component=0; component > coefficient_values (n_q_points, + Vector(n_components)); + coefficient->vector_value_list (fe_face_values_cell.get_quadrature_points(), + coefficient_values); + for (unsigned int component=0; component g(n_q_points); + vector > g(n_q_points, Vector(n_components)); neumann_bc.find(boundary_indicator)->second - ->value_list (fe_face_values_cell.get_quadrature_points(), - g); + ->vector_value_list (fe_face_values_cell.get_quadrature_points(), + g); - for (unsigned int point=0; point &JxW_values = fe_face_values_cell.get_JxW_values(); // take the square of the phi[i] - // for integration - transform (phi.begin(), phi.end(), - phi.begin(), ptr_fun(sqr)); + // for integration, and sum up + double face_integral = 0; + for (unsigned int component=0; component &fe_subface_values, FaceIntegrals &face_integrals, const Vector &solution, - const unsigned int n_components, - const unsigned int selected_component, - const Function *coefficient) { + const vector &component_mask, + const Function *coefficient) +{ + const unsigned int n_components = component_mask.size(); + const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); @@ -485,8 +575,12 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // points // // let psi be a short name for - // [a grad u_h] - vector > psi(n_q_points); + // [a grad u_h], where the second + // index be the component of the + // finite element, and the first + // index the number of the + // quadrature point + vector > > psi(n_q_points, vector >(n_components)); // store which number #cell# has in the // list of neighbors of #neighbor# @@ -518,39 +612,29 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // store the gradient of the solution // in psi fe_subface_values.reinit (cell, face_no, subface_no); - if (n_components == 1) - fe_subface_values.get_function_grads (solution, psi); - else - { - vector > > tmp (n_q_points, - vector >(n_components)); - fe_subface_values.get_function_grads (solution, tmp); - for (unsigned int i=0; i > neighbor_psi (n_q_points); + vector > > neighbor_psi (n_q_points, + vector >(n_components)); fe_face_values.reinit (neighbor_child, neighbor_neighbor); - if (n_components == 1) - fe_face_values.get_function_grads (solution, neighbor_psi); - else - { - vector > > tmp (n_q_points, - vector >(n_components)); - fe_face_values.get_function_grads (solution, tmp); - for (unsigned int i=0; i >()); - + for (unsigned int component=0; component phi(n_q_points,0); + vector > phi(n_q_points, vector(n_components)); const vector > &normal_vectors(fe_face_values. get_normal_vectors()); - for (unsigned int point=0; point coefficient_values (n_q_points); - coefficient->value_list (fe_face_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointn_components == 1) + { + vector coefficient_values (n_q_points); + coefficient->value_list (fe_face_values.get_quadrature_points(), + coefficient_values); + for (unsigned int component=0; component > coefficient_values (n_q_points, + Vector(n_components)); + coefficient->vector_value_list (fe_face_values.get_quadrature_points(), + coefficient_values); + for (unsigned int component=0; componentface(neighbor_neighbor)] - = inner_product (phi.begin(), phi.end(), - fe_face_values.get_JxW_values(). - begin(), - 0.0); + const vector &JxW_values = fe_face_values.get_JxW_values(); + + // take the square of the phi[i] + // for integration, and sum up + double face_integral = 0; + for (unsigned int component=0; componentface(neighbor_neighbor)] = face_integral; }; @@ -613,6 +714,6 @@ integrate_over_irregular_face (const active_cell_iterator &cell, -// explicit instantiations +// explicit instantiations template class KellyErrorEstimator;