From: Denis Davydov Date: Wed, 17 Jun 2015 21:31:55 +0000 (+0200) Subject: extend KellyErrorEstimator to calculate the boundary residual term of the reliable... X-Git-Tag: v8.3.0-rc1~103^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F1024%2Fhead;p=dealii.git extend KellyErrorEstimator to calculate the boundary residual term of the reliable hp estimator --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 3bd537bfe5..4a17b396d9 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -43,7 +43,7 @@ inconvenience this causes. but the scalar type of the vector was always double. This has been changed so that the two always need to match.
- (David Davydov, Wolfgang Bangerth, 2015/06/17) + (Denis Davydov, Wolfgang Bangerth, 2015/06/17)
  • Changed: Functions such as FEValuesBase::get_function_values() that @@ -64,7 +64,7 @@ inconvenience this causes. complex-valued data types. This includes compiling PETSc with complex scalars.
    - (David Davydov, 2015/05/08) + (Denis Davydov, 2015/05/08)
  • Removed: The generic, templated vmult, Tvmult, etc. -interfaces of @@ -385,6 +385,12 @@ inconvenience this causes.
      +
    1. New: Implement a modified version of the Kelly error estimator, which + effectively provides the boundary residual term for the hp-FEM error estimators. +
      + (Denis Davydov, 2015/06/17) +
      2. +
    2. New: DerivativeForm() now takes an additional optional template argument specifying the type, similarly to Tensor() classes.
      diff --git a/include/deal.II/numerics/error_estimator.h b/include/deal.II/numerics/error_estimator.h index 8d860b6d77..f0be09f769 100644 --- a/include/deal.II/numerics/error_estimator.h +++ b/include/deal.II/numerics/error_estimator.h @@ -41,17 +41,26 @@ namespace hp /** * Implementation of the error indicator by Kelly, De S. R. Gago, Zienkiewicz and - * Babuska. This error indicator tries to approximate the error per cell by + * Babuska and its modification for the hp-FEM. + * This error indicator 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 and Oden, "A Posteriori Error Estimation in Finite Element * Analysis" (Wiley, 2000) for a complete discussion. * - * @note In spite of the name, this is not truly an a posteriori error + * In the original Kelly error estimator, the contribution of each face to the + * cell error is scaled with the cell diagonal. In the modified version, + * however, we employ a scaling factor which depends on the face diagonal and + * polynomial degrees of the adjacent elements. The choice between the two is + * done by means of the enumerator, defined within the class. + * + * @note In spite of the name, Kelly estimator is not truly an a posteriori error * estimator, even if applied to the Poisson problem only. It gives good hints * for mesh refinement, but the estimate is not to be trusted. For higher * order trial spaces the integrals computed here tend to zero faster than the - * error itself, thus ruling out the values as error estimators. + * error itself, thus ruling out the values as error estimators. However, the + * modified version discussed below can be utilised to obtain the reliable + * error estimator by adding the residual (volume) part. * * The error estimator really only estimates the error for the generalized * Poisson equation $-\nabla\cdot a(x) \nabla u = f$ with either Dirichlet @@ -83,22 +92,24 @@ namespace hp *

      Implementation

      * * In principle, the implementation of the error estimation is simple: let \f[ - * \eta_K^2 = \frac h{24} \int_{\partial K} \left[a \frac{\partial + * \eta_K^2 = \sum_{F\in\partial K} c_F \int_{\partial K_F} \left[a \frac{\partial * u_h}{\partial n}\right]^2 do \f] 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$, but this factor is a bit esoteric, + * Ainsworth $ c_F=\frac h{24} $, but this factor is a bit esoteric, * stemming from interpolation estimates and stability constants which may * hold for the Poisson problem, but may not hold for more general situations. - * In the implementation, this factor is considered, but may lead to wrong - * results. You may scale the vector appropriately afterwards. + * Alternatively, we consider the case when $ c_F=\frac {h_F}_{2p_F} $, + * where $ h_F $ is face diagonal and $ p_F=max(p^+,p^-) $ is the + * maximum polynomial degree of adjacent elements. The choice between the two is + * done by means of the enumerator, provided as the last argument in all functions. * * To perform the integration, use is made of the FEFaceValues and * FESubfaceValues classes. The integration is performed by looping over all * cells and integrating over faces that are not yet treated. This way we * avoid integration on faces twice, once for each time we visit one of the * adjacent cells. In a second loop over all cells, we sum up the - * contributions of the faces (which are the integrated square of the jumps) - * of each cell and take the square root. + * contributions of the faces (which are the integrated square of the jumps times + * some factor) of each cell and take the square root. * * The integration is done using a quadrature formula on the face. For linear * trial functions (FEQ1), the QGauss2 or even the QMidpoint rule will @@ -107,16 +118,17 @@ namespace hp * * We store the contribution of each face in a @p map, as provided by the C++ * standard library, with the iterator pointing to that face being the key - * into the map. In fact, we do not store the indicator per face, but only the - * integral listed above. When looping the second time over all cells, we have - * to sum up the contributions of the faces, multiply them with $\frac h{24}$ - * and take the square root. By doing the multiplication with $h$ in the - * second loop, we avoid problems to decide with which $h$ to multiply, that - * of the cell on the one or that of the cell on the other side of the face. + * into the map. When looping the second time over all cells, + * we have to sum up the contributions of the faces and take the square root. + * For the Kelly estimator, the multiplication with $\frac h{24}$ is done + * in the second loop. By doing so we avoid problems to decide with which $h$ + * to multiply, that of the cell on the one or that of the cell on the other + * side of the face. Whereas for the hp-estimator the @p map stores integrals + * multiplied by $\frac {h_f}_{2p_f}$, which are then summed in the second loop. * - * $h$ is taken to be the greatest length of the diagonals of the cell. For - * more or less uniform cells without deformed angles, this coincides with the - * diameter of the cell. + * $h$ ($h_f$) is taken to be the greatest length of the diagonals of the cell (face). + * For more or less uniform cells (faces) without deformed angles, this coincides + * with the diameter of the cell (face). * * *

      Vector-valued functions

      @@ -164,9 +176,11 @@ namespace hp * different faces to the cells easier. * *
    3. The face belongs to a Neumann boundary. In this case, the - * contribution of the face $F\in\partial K$ looks like \f[ \int_F + * contribution of the face $F\in\partial K$ looks like \f[ n_F\int_F * \left|g-a\frac{\partial u_h}{\partial n}\right|^2 ds \f] where $g$ is the - * Neumann boundary function. If the finite element is vector-valued, then + * Neumann boundary function, $n_F=\frac h_{24}$ and $n_F=\frac {h_F}{p}$ for + * the Kelly and hp-estimator, respectively. + * If the finite element is vector-valued, then * obviously the function denoting the Neumann boundary conditions needs to be * vector-valued as well. * @@ -233,13 +247,25 @@ namespace hp * that accepts several in- and output vectors at the same time. * * @ingroup numerics - * @author Wolfgang Bangerth, 1998, 1999, 2000, 2004, 2006; parallelization by - * Thomas Richter, 2000 + * @author Wolfgang Bangerth, 1998, 1999, 2000, 2004, 2006, Denis Davydov, 2015; + * parallelization by Thomas Richter, 2000 */ template class KellyErrorEstimator { public: + /** + * The enum type given to the class functions to decide on the scaling factors + * of the facial integrals. + */ + enum Strategy + { + //! Kelly error estimator with the factor $\frac h_{24}$. + cell_diameter_over_24 = 0, + //! the boundary residual estimator with the factor $\frac {h_F}_{2 max(p^+,p^-)$. + face_diameter_over_twice_max_degree + }; + /** * Implementation of the error estimator described above. You may give a * coefficient, but there is a default value which denotes the constant @@ -285,6 +311,9 @@ public: * the number of threads determined automatically. The parameter is retained * for compatibility with old versions of the library. * + * The @p strategy parameter is used to choose the scaling factor for + * the integral over cell's faces. + * * @note If the DoFHandler object given as an argument to this function * builds on a parallel::distributed::Triangulation, this function skips * computations on all cells that are not locally owned. In that case, the @@ -308,7 +337,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** * Calls the @p estimate function, see above, with @@ -324,7 +354,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** * Same function as above, but accepts more than one solution vector and @@ -350,7 +381,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** * Calls the @p estimate function, see above, with @@ -366,7 +398,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** @@ -384,7 +417,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** @@ -401,7 +435,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** @@ -419,7 +454,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** @@ -436,8 +472,8 @@ public: const Function *coefficients = 0, const unsigned int n_threads = numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id, - const types::material_id material_id = numbers::invalid_material_id); - + const types::material_id material_id = numbers::invalid_material_id, + const KellyErrorEstimator::Strategy strategy = cell_diameter_over_24); /** * Exception diff --git a/source/numerics/error_estimator.cc b/source/numerics/error_estimator.cc index f73a36ead9..ae5bafe731 100644 --- a/source/numerics/error_estimator.cc +++ b/source/numerics/error_estimator.cc @@ -464,6 +464,90 @@ namespace internal return face_integral; } + /** + * A factor to scale the integral for the face at the boundary. + * Used for Neumann BC. + */ + template + double boundary_face_factor(const typename DH::active_cell_iterator &cell, + const unsigned int face_no, + const dealii::hp::FEFaceValues &fe_face_values_cell, + const typename KellyErrorEstimator::Strategy strategy) + { + switch (strategy) + { + case KellyErrorEstimator::cell_diameter_over_24 : + return 1.0; + case KellyErrorEstimator::face_diameter_over_twice_max_degree : + const double cell_degree = fe_face_values_cell.get_fe_collection()[cell->active_fe_index()].degree; + return cell->face(face_no)->diameter() / cell_degree; + } + } + + + /** + * A factor to scale the integral for the regular face. + */ + template + double regular_face_factor(const typename DH::active_cell_iterator &cell, + const unsigned int face_no, + const dealii::hp::FEFaceValues &fe_face_values_cell, + const dealii::hp::FEFaceValues &fe_face_values_neighbor, + const typename KellyErrorEstimator::Strategy strategy) + { + switch (strategy) + { + case KellyErrorEstimator::cell_diameter_over_24 : + return 1.0; + case KellyErrorEstimator::face_diameter_over_twice_max_degree : + const double cell_degree = fe_face_values_cell.get_fe_collection()[cell->active_fe_index()].degree; + const double neighbor_degree = fe_face_values_neighbor.get_fe_collection()[cell->neighbor(face_no)->active_fe_index()].degree; + return cell->face(face_no)->diameter() / std::max(cell_degree,neighbor_degree) / 2.0; + } + } + + /** + * A factor to scale the integral for the irregular face. + */ + template + double irregular_face_factor(const typename DH::active_cell_iterator &cell, + const typename DH::active_cell_iterator &neighbor_child, + const unsigned int face_no, + const unsigned int subface_no, + const dealii::hp::FEFaceValues &fe_face_values, + dealii::hp::FESubfaceValues &fe_subface_values, + const typename KellyErrorEstimator::Strategy strategy) + { + switch (strategy) + { + case KellyErrorEstimator::cell_diameter_over_24 : + return 1.0; + case KellyErrorEstimator::face_diameter_over_twice_max_degree : + const double cell_degree = fe_face_values.get_fe_collection()[cell->active_fe_index()].degree; + const double neighbor_child_degree = fe_subface_values.get_fe_collection()[neighbor_child->active_fe_index()].degree; + return cell->face(face_no)->child(subface_no)->diameter()/std::max(neighbor_child_degree,cell_degree)/2.0; + } + } + + /** + * A factor used when summing up all the contribution + * from different faces of each cell. + */ + template + double cell_factor(const typename DH::active_cell_iterator &cell, + const unsigned int face_no, + const DH &dof_handler, + const typename KellyErrorEstimator::Strategy strategy) + { + switch (strategy) + { + case KellyErrorEstimator::cell_diameter_over_24 : + return cell->diameter()/24; + case KellyErrorEstimator::face_diameter_over_twice_max_degree : + return 1.0; + } + } + /** @@ -480,7 +564,8 @@ namespace internal const typename DH::active_cell_iterator &cell, const unsigned int face_no, dealii::hp::FEFaceValues &fe_face_values_cell, - dealii::hp::FEFaceValues &fe_face_values_neighbor) + dealii::hp::FEFaceValues &fe_face_values_neighbor, + const typename KellyErrorEstimator::Strategy strategy) { const unsigned int dim = DH::dimension; (void)dim; @@ -500,6 +585,7 @@ namespace internal fe_face_values_cell.get_present_fe_values() .get_function_gradients (*solutions[n], parallel_data.psi[n]); + double factor; // now compute over the other side of the face if (face->at_boundary() == false) // internal face; integrate jump of gradient across this face @@ -522,6 +608,10 @@ namespace internal fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor, cell->active_fe_index()); + factor = regular_face_factor(cell,face_no, + fe_face_values_cell,fe_face_values_neighbor, + strategy); + // get gradients on neighbor cell for (unsigned int n=0; n(cell,face_no, + fe_face_values_cell, + strategy); + } // now go to the generic function that does all the other things local_face_integrals[face] = integrate_over_face (parallel_data, face, fe_face_values_cell); + for (unsigned int i = 0; i < local_face_integrals[face].size(); i++) + local_face_integrals[face][i] *= factor; } @@ -557,7 +655,8 @@ namespace internal const typename DH::active_cell_iterator &cell, const unsigned int face_no, dealii::hp::FEFaceValues &fe_face_values, - dealii::hp::FESubfaceValues &fe_subface_values) + dealii::hp::FESubfaceValues &fe_subface_values, + const typename KellyErrorEstimator::Strategy strategy) { const unsigned int dim = DH::dimension; (void)dim; @@ -603,6 +702,14 @@ namespace internal fe_face_values.reinit (neighbor_child, neighbor_neighbor, cell->active_fe_index()); + const double factor = irregular_face_factor(cell, + neighbor_child, + face_no, + subface_no, + fe_face_values, + fe_subface_values, + strategy); + // store the gradient of the solution in psi for (unsigned int n=0; nface(neighbor_neighbor)] = integrate_over_face (parallel_data, face, fe_face_values); + for (unsigned int i = 0; i < local_face_integrals[neighbor_child->face(neighbor_neighbor)].size(); i++) + local_face_integrals[neighbor_child->face(neighbor_neighbor)][i] *= factor; } // finally loop over all subfaces to collect the contributions of the @@ -648,7 +757,8 @@ namespace internal estimate_one_cell (const typename DH::active_cell_iterator &cell, ParallelData ¶llel_data, std::map > &local_face_integrals, - const std::vector &solutions) + const std::vector &solutions, + const typename KellyErrorEstimator::Strategy strategy) { const unsigned int dim = DH::dimension; const unsigned int n_solution_vectors = solutions.size(); @@ -769,7 +879,8 @@ namespace internal local_face_integrals, cell, face_no, parallel_data.fe_face_values_cell, - parallel_data.fe_face_values_neighbor); + parallel_data.fe_face_values_neighbor, + strategy); else // otherwise we need to do some special computations which do not @@ -779,7 +890,8 @@ namespace internal local_face_integrals, cell, face_no, parallel_data.fe_face_values_cell, - parallel_data.fe_subface_values); + parallel_data.fe_subface_values, + strategy); } } } @@ -805,13 +917,14 @@ estimate (const Mapping &mapping, const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { // just pass on to the other function const std::vector solutions (1, &solution); std::vector*> errors (1, &error); estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id, material_id); + component_mask, coefficients, n_threads, subdomain_id, material_id, strategy); } @@ -828,11 +941,12 @@ estimate (const DH &dof_handler, const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { estimate(StaticMappingQ1::mapping, dof_handler, quadrature, neumann_bc, solution, error, component_mask, coefficients, n_threads, - subdomain_id, material_id); + subdomain_id, material_id, strategy); } @@ -850,13 +964,14 @@ estimate (const Mapping &mapping, const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { // just pass on to the other function const std::vector solutions (1, &solution); std::vector*> errors (1, &error); estimate (mapping, dof_handler, quadrature, neumann_bc, solutions, errors, - component_mask, coefficients, n_threads, subdomain_id, material_id); + component_mask, coefficients, n_threads, subdomain_id, material_id, strategy); } @@ -873,11 +988,12 @@ estimate (const DH &dof_handler, const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { estimate(StaticMappingQ1::mapping, dof_handler, quadrature, neumann_bc, solution, error, component_mask, coefficients, n_threads, - subdomain_id, material_id); + subdomain_id, material_id, strategy); } @@ -897,7 +1013,8 @@ estimate (const Mapping &mapping, const Function *coefficients, const unsigned int , const types::subdomain_id subdomain_id_, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { #ifdef DEAL_II_WITH_P4EST if (dynamic_cast*> @@ -985,7 +1102,7 @@ estimate (const Mapping &mapping, WorkStream::run (dof_handler.begin_active(), static_cast(dof_handler.end()), std_cxx11::bind (&internal::estimate_one_cell, - std_cxx11::_1, std_cxx11::_2, std_cxx11::_3, std_cxx11::ref(solutions)), + std_cxx11::_1, std_cxx11::_2, std_cxx11::_3, std_cxx11::ref(solutions),strategy), std_cxx11::bind (&internal::copy_local_to_global, std_cxx11::_1, std_cxx11::ref(face_integrals)), parallel_data, @@ -1022,7 +1139,10 @@ estimate (const Mapping &mapping, Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(), ExcInternalError()); - const double factor = cell->diameter() / 24; + const double factor = internal::cell_factor(cell, + face_no, + dof_handler, + strategy); for (unsigned int n=0; n &mapping, const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { // forward to the function with the QCollection estimate (mapping, dof_handler, hp::QCollection(quadrature), neumann_bc, solutions, errors, component_mask, coefficients, - n_threads, subdomain_id, material_id); + n_threads, subdomain_id, material_id, strategy); } @@ -1079,11 +1200,12 @@ void KellyErrorEstimator::estimate (const DH const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { estimate(StaticMappingQ1::mapping, dof_handler, quadrature, neumann_bc, solutions, errors, component_mask, coefficients, n_threads, - subdomain_id, material_id); + subdomain_id, material_id, strategy); } @@ -1099,11 +1221,12 @@ void KellyErrorEstimator::estimate (const DH const Function *coefficients, const unsigned int n_threads, const types::subdomain_id subdomain_id, - const types::material_id material_id) + const types::material_id material_id, + const typename KellyErrorEstimator::Strategy strategy) { estimate(StaticMappingQ1::mapping, dof_handler, quadrature, neumann_bc, solutions, errors, component_mask, coefficients, n_threads, - subdomain_id, material_id); + subdomain_id, material_id, strategy); } diff --git a/source/numerics/error_estimator.inst.in b/source/numerics/error_estimator.inst.in index 7aa049106f..6d13601f9a 100644 --- a/source/numerics/error_estimator.inst.in +++ b/source/numerics/error_estimator.inst.in @@ -38,7 +38,8 @@ estimate > (const Mapping *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -53,7 +54,8 @@ estimate > ( const Function *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -68,7 +70,8 @@ estimate > (const Mapping *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -83,7 +86,8 @@ estimate > ( const Function *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -98,7 +102,8 @@ estimate > (const Mapping *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -113,7 +118,8 @@ estimate > ( const Function *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -128,7 +134,8 @@ estimate > (const Mapping *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); template void @@ -143,7 +150,8 @@ estimate > ( const Function *, const unsigned int , const unsigned int , - const types::material_id); + const types::material_id, + const KellyErrorEstimator::Strategy); #endif } diff --git a/tests/deal.II/error_estimator_02.cc b/tests/deal.II/error_estimator_02.cc new file mode 100644 index 0000000000..a1ad203e30 --- /dev/null +++ b/tests/deal.II/error_estimator_02.cc @@ -0,0 +1,769 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +/* Author: Denis Davydov, University of Erlangen-Nuremberg, 2015 */ +// make sure that the modified Kelly error estimator with +// strategy = face_diameter_over_twice_max_degree +// returns correct results. +// To that end consider 3 cases at domain [0,2]^dim and one at domain [0,1]^dim: +// +// 1) +// Neuman BC g(x) = c; // constant +// single element with p =3; u=0; +// \eta = h * A * c^2 / p; +// +// +// 2) +// +// --------------- +// | | | +// | | | +// | p1 | p2 | +// | | | +// | | | +// --------------- +// +// u (left) = 0; +// u (right) = kx; +// +// \eta = h * A * k^2 / 2 max(p1,p2) +// +// +// 3) +// ---------------------- +// | | | | +// | p3 | p3 f_1 | +// |------------| p1 | +// | p2 | p2 f_2 | +// | | | | +// ---------------------- +// solution is the same as above; +// +// \eta_1 = h * A * k^2 / 2 max(p3,p1) +// \eta_2 = h * A * k^2 / 2 max(p2,p1) +// \eta_3 = \eta_1 + \eta_2 +// +// 4) just arbitrary mesh +// ------------------- +// | | | +// | 1 | 1 | +// | | | +// |------------------ +// | 3 | 2 | | +// |--------| 1 | +// | 3 | 2 | | +// ------------------ +// +// and evaluate error of the interpolated to it function. + + + +#include "../tests.h" + +// dealii +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// c++ +#include +#include + +using namespace dealii; + +template +class MyFunction : public dealii::Function +{ +public: + MyFunction(const double k); + + virtual double value(const dealii::Point & point, + const unsigned int component = 0 ) const; + + double get_k() const; + +private: + const double k; +}; + +template +MyFunction::MyFunction(const double k) +: +dealii::Function(1), +k(k) +{ + +} + +template +double MyFunction::value(const dealii::Point & point, + const unsigned int ) const +{ + const double x = point[0]-1.0; + + if (x < 0) + return 0.0; + else + return k * x; +} + +template +double MyFunction::get_k() const +{ + return k; +} + +// neuman bc +template +class NeumanBC : public dealii::Function +{ +public: + NeumanBC(const double c); + + virtual double value(const dealii::Point & point, + const unsigned int component = 0 ) const; + + double get_c() const; + +private: + const double c; +}; + +template +NeumanBC::NeumanBC(const double c) +: +dealii::Function(1), +c(c) +{ +} + +template +double NeumanBC::value(const dealii::Point & point, + const unsigned int ) const +{ + return c; +} + +template +double NeumanBC::get_c() const +{ + return c; +} + +// helper function to get diagonal and +// area of the squared element with lenght h +template +void get_h_area(double &h, double &a, const double L); + +template<> +void get_h_area<2>(double &h, double &a, const double L) +{ + h = L; + a = L; +} + +template<> +void get_h_area<3>(double &h, double &a, const double L) +{ + h = std::sqrt(2.0)*L; + a = L*L; +} + +// helper function to get diagonal and area of the +// h-refined face. +template +void get_h_area_sub(double &h, double &a, const double L); + +template<> +void get_h_area_sub<2>(double &h, double &a, const double L) +{ + h = L/2; + a = L/2; +} + +template<> +void get_h_area_sub<3>(double &h, double &a, const double L) +{ + h = std::sqrt(2.0)*L/2; + a = L*L/4.0; +} + +// output for inspection +template +void output(const std::string name, + const Triangulation &triangulation, + const hp::DoFHandler &dof_handler, + const Vector &values, + const Vector &error) +{ + dealii::Vector fe_degrees(triangulation.n_active_cells()); + { + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + fe_degrees(index) = dof_handler.get_fe()[cell->active_fe_index()].degree; + } + + dealii::DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector(values, + std::string("function_interpolation")); + data_out.add_data_vector(fe_degrees, + std::string("fe_degree")); + data_out.add_data_vector(error, + std::string("error")); + data_out.build_patches (); + + std::ofstream output (name); + data_out.write_vtu (output); +} + +// case 1) +template +void test_neumann(const NeumanBC &func) +{ + deallog << "NeumanBC case:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p = 3; + + fe_collection.push_back(dealii::FE_Q(p)); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + + const double L = 2.0; + + // set-up domain + { + GridGenerator::hyper_cube (triangulation,.0,L,/*colorize*/true); + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector values(dof_handler.n_dofs()); + values = 0.0; + + dealii::deallog << "dof values:"<::type function_map; + function_map[0] = &func; + + dealii::Vector error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + function_map, + values, + error, + dealii::ComponentMask(), + 0, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + const double expected_value_squared = h*A*std::pow(func.get_c(),2)/p; + dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl; + + AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[0] ) < 1e-5, dealii::ExcInternalError()); + + dof_handler.clear(); +} + +// case 2) +template +void test_regular(const MyFunction &func) +{ + deallog << std::endl; + deallog << "Regular face:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p1 = 1; + const unsigned int p2 = 2; + std::vector p_degree; + p_degree.push_back(p1); + p_degree.push_back(p2); + + for (unsigned int i=0;i(p)); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + } + + const double L = 2.0; + + // set-up domain + { + std::vector repetitions(dim,1); + repetitions[0] = 2; + dealii::Point p1; + dealii::Point p2; + for (unsigned int d = 0; d < dim; d++) + { + p1[d] = 0.0; + p2[d] = L; + } + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + p1, + p2, + /*colorize*/ false); + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (;cell != endc; cell++) + if (cell->center()[0] > 1.0) + { + cell->set_active_fe_index(1); + break; + } + + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector values(dof_handler.n_dofs()); + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap::type(), + values, + error, + dealii::ComponentMask(), + 0, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + const double expected_value_squared = h*A*std::pow(func.get_k(),2)/2.0/std::max(p1,p2); + dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl; + for (unsigned int i = 0; i < error.size(); i++) + AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[i] ) < 1e-6, dealii::ExcInternalError()); + + dof_handler.clear(); +} + +// case 3) +template +void test_irregular(const MyFunction &func) +{ + deallog << std::endl; + deallog << "Irregular face:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p1 = 1; + const unsigned int p2 = 2; + const unsigned int p3 = 3; + std::vector p_degree; + p_degree.push_back(p1); + p_degree.push_back(p2); + p_degree.push_back(p3); + + for (unsigned int i=0;i(p)); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + } + + const double L = 2.0; + + // set-up domain + { + std::vector repetitions(dim,1); + repetitions[0] = 2; + dealii::Point p1; + dealii::Point p2; + for (unsigned int d = 0; d < dim; d++) + { + p1[d] = 0.0; + p2[d] = L; + } + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + p1, + p2, + /*colorize*/ false); + // refine left side + { + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + } + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (;cell != endc; cell++) + if (cell->center()[0] > 1.0) // right + { + cell->set_active_fe_index(0); + } + else if (cell->center()[1] > 1.0) // top + { + cell->set_active_fe_index(2); + } + else + { + cell->set_active_fe_index(1); + } + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector values(dof_handler.n_dofs()); + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap::type(), + values, + error, + dealii::ComponentMask(), + 0, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + // + const double expected_squared_1 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p3,p1); + const double expected_squared_2 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p2,p1); + const double expected_squared_3 = (dim==2) ? + expected_squared_1 + expected_squared_2: + 2*expected_squared_1 + 2*expected_squared_2; + + std::vector expected_error(error.size(),0.0); + + expected_error[0] = std::sqrt(expected_squared_3); + expected_error[2] = std::sqrt(expected_squared_2); + expected_error[4] = std::sqrt(expected_squared_1); + + if (dim ==3) + { + expected_error[6] = expected_error[2]; + expected_error[8] = expected_error[4]; + } + + dealii::deallog << "expected:"<< std::endl; + for (unsigned int i = 0; i < expected_error.size(); i++) + deallog<<" " << expected_error[i]; + deallog < +class MySecondFunction : public dealii::Function +{ +public: + MySecondFunction(); + + virtual double value(const dealii::Point & point, + const unsigned int component = 0 ) const; +}; + +template +MySecondFunction::MySecondFunction() +: +dealii::Function(1) +{ + +} + +template +double MySecondFunction::value(const dealii::Point & point, + const unsigned int ) const +{ + double f = 0.0; + const double &x = point[0]; + Assert (dim>1, dealii::ExcNotImplemented()); + const double &y = point[1]; + + return (1.-x)*(1.-y)*(1.-y)+std::pow(1.0-y,4)*std::exp(-x); +} + +template +void test(const MySecondFunction &func) +{ + deallog << std::endl; + deallog << "More complicated mesh:"< triangulation; + dealii::hp::DoFHandler dof_handler(triangulation); + dealii::hp::FECollection fe_collection; + dealii::hp::QCollection quadrature_formula; + dealii::hp::QCollection face_quadrature_formula; + dealii::ConstraintMatrix constraints; + for (unsigned int p = 1; p <=3; p++) + { + fe_collection.push_back(dealii::FE_Q(p)); + quadrature_formula.push_back(dealii::QGauss(p+1)); + face_quadrature_formula.push_back(dealii::QGauss(p+1)); + } + dealii::GridGenerator::hyper_cube (triangulation,0.0,1.0); // reference cell + + // refine + { + triangulation.refine_global(1); + + // otherwise the next set_active_fe_index + // will not carry to the child cells. + dof_handler.distribute_dofs (fe_collection); + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (;cell != endc; cell++) + { + bool in_top_left = true; + for (unsigned int d=0; d< dim; d++) + in_top_left = in_top_left && (cell->center()[d] < 0.5); + + if (in_top_left) + { + cell->set_active_fe_index(1); + cell->set_refine_flag(); + break; + } + } + + triangulation.prepare_coarsening_and_refinement(); + + triangulation.execute_coarsening_and_refinement(); + + cell = dof_handler.begin_active(); + + for (;cell != endc; cell++) + { + if (cell->center()[0] < 0.25) + { + cell->set_active_fe_index(2); + } + } + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector values(dof_handler.n_dofs()); + + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap::type (), + values, + error, + dealii::ComponentMask(), + 0, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"< func(1.25); + test_neumann(func); + } + + { + MyFunction<2> func(0.25); + test_regular(func); + } + + { + MyFunction<2> func(0.75); + test_irregular(func); + } + + deallog << "===3d==="< func(1.25); + test_neumann(func); + } + + { + MyFunction<3> func(0.25); + test_regular(func); + } + + { + MyFunction<3> func(0.75); + test_irregular(func); + } + + { + MySecondFunction<2> function; + test(function); + } + + + + dealii::deallog << "Ok"<