From 0bf751a9a312a765a542723dd763c89539104cda Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 31 Jul 1998 15:00:12 +0000 Subject: [PATCH] Extend error estimator for variable coefficients. git-svn-id: https://svn.dealii.org/trunk@465 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/error_estimator.h | 29 +++++-- .../source/numerics/error_estimator.cc | 79 +++++++++++++------ 2 files changed, 77 insertions(+), 31 deletions(-) diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 6783955f48..ee1a03a3a0 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -36,7 +36,12 @@ class dVector; * integrals computed here show superconvergence properties, i.e. they tend * to zero faster than the error itself, thus ruling out the values as error * indicators. - * + * + * The error estimator really only estimates the error for the generalized + * Poisson equation $-\nabla\cdot a(x) \nabla u = f$ with either Dirichlet + * boundary conditions or generalized Neumann boundary conditions involving + * the conormal derivative $a\frac{du}{dn} = g$. + * * The error estimator returns a vector of estimated errors per cell which * can be used to feed the #Triangulation::refine_*# functions. * @@ -45,7 +50,7 @@ class dVector; * * In principle, the implementation of the error estimation is simple: let * $$ \eta_K^2 = - * \frac h{24} \int_{\partial K} \left[\frac{\partial u_h}{\partial n}\right]^2 do + * \frac h{24} \int_{\partial K} \left[a \frac{\partial u_h}{\partial n}\right]^2 do * $$ * be the error estimator for cell $K$. $[\cdot]$ denotes the jump of the * argument at the face. In the paper of Ainsworth, $h$ is divided by $24$, @@ -107,7 +112,7 @@ class dVector; * * \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-\frac{\partial u_h}{\partial n}\right| ds $$ + * $$ \int_F \left|g-a\frac{\partial u_h}{\partial n}\right|^2 ds $$ * where $g$ is the Neumann boundary function. * * \item No other boundary conditions are considered. @@ -156,14 +161,22 @@ class KellyErrorEstimator { * guaranteed by the #map# data type. */ 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. + */ static void estimate_error (const DoFHandler &dof, const Quadrature &quadrature, const FiniteElement &fe, const Boundary &boundary, const FunctionMap &neumann_bc, const dVector &solution, - dVector &error); + dVector &error, + const Function *coefficient = 0); /** * Exception @@ -220,7 +233,8 @@ class KellyErrorEstimator { FEFaceValues &fe_face_values_cell, FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, - const dVector &solution); + const dVector &solution, + const Function *coefficient); /** * The same applies as for the function @@ -238,7 +252,8 @@ class KellyErrorEstimator { FEFaceValues &fe_face_values, FESubfaceValues &fe_subface_values, FaceIntegrals &face_integrals, - const dVector &solution); + const dVector &solution, + 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 bbc196654f..a03e614147 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -32,7 +32,8 @@ void KellyErrorEstimator<1>::estimate_error (const DoFHandler<1> &, const Boundary<1> &, const FunctionMap &, const dVector &, - dVector &) { + dVector &, + const Function<1> *) { Assert(false, ExcNotImplemented()); }; @@ -46,7 +47,8 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, const Boundary &boundary, const FunctionMap &neumann_bc, const dVector &solution, - dVector &error) { + dVector &error, + const Function *coefficient) { Assert (neumann_bc.find(255) == neumann_bc.end(), ExcInvalidBoundaryIndicator()); @@ -138,7 +140,8 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, fe_face_values_cell, fe_face_values_neighbor, face_integrals, - solution); + solution, + coefficient); else // otherwise we need to do some // special computations which do @@ -148,7 +151,8 @@ void KellyErrorEstimator::estimate_error (const DoFHandler &dof, n_q_points, fe_face_values_cell, fe_subface_values, - face_integrals, solution); + face_integrals, solution, + coefficient); }; @@ -190,7 +194,8 @@ void KellyErrorEstimator<1>::integrate_over_regular_face (const active_cell_iter FEFaceValues<1> &, FEFaceValues<1> &, FaceIntegrals &, - const dVector &) { + const dVector &, + const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -206,7 +211,8 @@ integrate_over_irregular_face (const active_cell_iterator &, FEFaceValues<1> &, FESubfaceValues<1> &, FaceIntegrals &, - const dVector &) { + const dVector &, + const Function<1> *) { Assert (false, ExcInternalError()); }; @@ -225,7 +231,8 @@ integrate_over_regular_face (const active_cell_iterator &cell, FEFaceValues &fe_face_values_cell, FEFaceValues &fe_face_values_neighbor, FaceIntegrals &face_integrals, - const dVector &solution) { + const dVector &solution, + const Function *coefficient) { const DoFHandler::face_iterator face = cell->face(face_no); // initialize data of the restriction @@ -238,7 +245,7 @@ integrate_over_regular_face (const active_cell_iterator &cell, // points // // let psi be a short name for - // [grad u_h] + // [a grad u_h] vector > psi(n_q_points); fe_face_values_cell.get_function_grads (solution, psi); @@ -311,7 +318,7 @@ integrate_over_regular_face (const active_cell_iterator &cell, for (unsigned int point=0; pointat_boundary() == true) // neumann boundary face. compute @@ -334,11 +341,22 @@ integrate_over_regular_face (const active_cell_iterator &cell, phi[point] -= g[point]; }; + // if a coefficient was given: use that + // to scale the jump in the gradient + if (coefficient != 0) + { + vector coefficient_values (n_q_points); + coefficient->value_list (fe_face_values_cell.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; point &fe_face_values, FESubfaceValues &fe_subface_values, FaceIntegrals &face_integrals, - const dVector &solution) { + const dVector &solution, + const Function *coefficient) { const DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); Assert (neighbor.state() == valid, ExcInternalError()); Assert (neighbor->has_children(), ExcInternalError()); @@ -378,7 +397,7 @@ integrate_over_irregular_face (const active_cell_iterator &cell, // points // // let psi be a short name for - // [grad u_h] + // [a grad u_h] vector > psi(n_q_points); // store which number #cell# has in the @@ -428,23 +447,35 @@ integrate_over_irregular_face (const active_cell_iterator &cell, psi.begin(), minus >()); - // next we have to multiply this with - // the normal vector. Since we have - // taken the difference of gradients - // for internal faces, we may chose - // the normal vector of one cell, - // taking that of the neighbor - // would only change the sign. We take - // the outward normal. - // - // let phi be the name of the integrand + // next we have to multiply this with + // the normal vector. Since we have + // taken the difference of gradients + // for internal faces, we may chose + // the normal vector of one cell, + // taking that of the neighbor + // would only change the sign. We take + // the outward normal. + // + // let phi be the name of the integrand vector phi(n_q_points,0); 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; point