From 4cee6eb1f5ffb7e88eab7e9c36cb2bfb91548039 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 28 Feb 2018 16:00:32 +0100 Subject: [PATCH] make KellyErrorEstimator<1> consistent with KellyErrorEstimator<2|3> --- include/deal.II/numerics/error_estimator.h | 40 +++++++++++++++++----- source/numerics/error_estimator_1d.cc | 37 ++++++++++++-------- source/numerics/error_estimator_1d.inst.in | 24 ++++++++----- 3 files changed, 70 insertions(+), 31 deletions(-) diff --git a/include/deal.II/numerics/error_estimator.h b/include/deal.II/numerics/error_estimator.h index 0084909464..e25438cd47 100644 --- a/include/deal.II/numerics/error_estimator.h +++ b/include/deal.II/numerics/error_estimator.h @@ -253,7 +253,7 @@ class KellyErrorEstimator public: /** * The enum type given to the class functions to decide on the scaling - * factors of the facial integrals. + * factors of the face integrals. */ enum Strategy { @@ -549,6 +549,20 @@ template class KellyErrorEstimator<1,spacedim> { public: + /** + * The enum type given to the class functions to decide on the scaling + * factors of the face integrals. + */ + enum Strategy + { + //! Kelly error estimator with the factor $\frac {h_K}{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, + //! Kelly error estimator with the factor $h_K$. + cell_diameter + }; + /** * Implementation of the error estimator described above. You may give a * coefficient, but there is a default value which denotes the constant @@ -583,7 +597,8 @@ public: const Function *coefficient = nullptr, 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 Strategy strategy = cell_diameter_over_24); /** * Call the @p estimate function, see above, with @@ -600,7 +615,8 @@ public: const Function *coefficients = nullptr, 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 Strategy strategy = cell_diameter_over_24); /** * Same function as above, but accepts more than one solution vectors and @@ -627,7 +643,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 Strategy strategy = cell_diameter_over_24); /** * Call the @p estimate function, see above, with @@ -644,7 +661,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 Strategy strategy = cell_diameter_over_24); /** @@ -663,7 +681,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 Strategy strategy = cell_diameter_over_24); /** @@ -681,7 +700,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 Strategy strategy = cell_diameter_over_24); /** @@ -700,7 +720,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 Strategy strategy = cell_diameter_over_24); /** @@ -718,7 +739,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 Strategy strategy = cell_diameter_over_24); /** * Exception diff --git a/source/numerics/error_estimator_1d.cc b/source/numerics/error_estimator_1d.cc index 71137c2fe4..6a62fdaec7 100644 --- a/source/numerics/error_estimator_1d.cc +++ b/source/numerics/error_estimator_1d.cc @@ -65,13 +65,14 @@ estimate (const Mapping<1,spacedim> &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 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); } @@ -89,10 +90,11 @@ estimate (const DoFHandlerType &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 Strategy strategy) { estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id, material_id); + error, component_mask, coefficients, n_threads, subdomain_id, material_id,strategy); } @@ -110,10 +112,11 @@ estimate (const DoFHandlerType &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 Strategy strategy) { estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id, material_id); + errors, component_mask, coefficients, n_threads, subdomain_id, material_id,strategy); } @@ -132,13 +135,14 @@ estimate (const Mapping<1,spacedim> &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 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); } @@ -155,10 +159,11 @@ estimate (const DoFHandlerType &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 Strategy strategy) { estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solution, - error, component_mask, coefficients, n_threads, subdomain_id, material_id); + error, component_mask, coefficients, n_threads, subdomain_id, material_id,strategy); } @@ -176,10 +181,11 @@ estimate (const DoFHandlerType &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 Strategy strategy) { estimate(StaticMappingQ1<1,spacedim>::mapping, dof_handler, quadrature, neumann_bc, solutions, - errors, component_mask, coefficients, n_threads, subdomain_id, material_id); + errors, component_mask, coefficients, n_threads, subdomain_id, material_id,strategy); } @@ -198,7 +204,8 @@ estimate (const Mapping<1,spacedim> & /*mapping*/, const Function * /*coefficient*/, const unsigned int, const types::subdomain_id /*subdomain_id*/, - const types::material_id /*material_id*/) + const types::material_id /*material_id*/, + const Strategy /*strategy*/) { Assert (false, ExcInternalError()); } @@ -218,8 +225,10 @@ estimate (const Mapping<1,spacedim> &mapping, const Function *coefficient, const unsigned int, const types::subdomain_id subdomain_id_, - const types::material_id material_id) + const types::material_id material_id, + const Strategy strategy) { + AssertThrow(strategy==cell_diameter_over_24, ExcNotImplemented()); typedef typename InputVector::value_type number; #ifdef DEAL_II_WITH_P4EST if (dynamic_cast*> diff --git a/source/numerics/error_estimator_1d.inst.in b/source/numerics/error_estimator_1d.inst.in index 47f69aed1b..5b1ba24d36 100644 --- a/source/numerics/error_estimator_1d.inst.in +++ b/source/numerics/error_estimator_1d.inst.in @@ -31,7 +31,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -46,7 +47,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -61,7 +63,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -76,7 +79,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -91,7 +95,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -106,7 +111,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -121,7 +127,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); template void @@ -136,7 +143,8 @@ for (VEC : VECTOR_TYPES ; deal_II_dimension : DIMENSIONS; deal_II_space_dimensio const Function *, const unsigned int, const types::subdomain_id, - const types::material_id); + const types::material_id, + const Strategy); #endif } -- 2.39.5