From: Marc Fehling Date: Fri, 20 Jun 2025 12:50:27 +0000 (+0200) Subject: Use MappingCollection for hp variant of KellyErrorEstimator. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F18604%2Fhead;p=dealii.git Use MappingCollection for hp variant of KellyErrorEstimator. --- diff --git a/doc/news/changes/incompatibilities/20250620Fehling b/doc/news/changes/incompatibilities/20250620Fehling new file mode 100644 index 0000000000..a30c5bd3b8 --- /dev/null +++ b/doc/news/changes/incompatibilities/20250620Fehling @@ -0,0 +1,5 @@ +Changed: The hp-variants of KellyErrorEstimator::estimate now take +hp::MappingCollection objects. The versions that require Mapping objects +are now deprecated. +
+(Marc Fehling, 2025/06/20) diff --git a/include/deal.II/numerics/error_estimator.h b/include/deal.II/numerics/error_estimator.h index 7738becbff..fcc9b34c00 100644 --- a/include/deal.II/numerics/error_estimator.h +++ b/include/deal.II/numerics/error_estimator.h @@ -43,7 +43,10 @@ namespace hp { template class QCollection; -} + + template + class MappingCollection; +} // namespace hp #endif @@ -427,11 +430,37 @@ public: /** * Equivalent to the set of functions above, except that this one takes a - * quadrature collection for hp-finite element dof handlers. + * mapping and quadrature collection for hp-finite element dof handlers. */ template static void estimate( + const hp::MappingCollection &mapping, + const DoFHandler &dof, + const hp::QCollection &quadrature, + const std::map *> + &neumann_bc, + const ReadVector &solution, + Vector &error, + const ComponentMask &component_mask = {}, + 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 Strategy strategy = cell_diameter_over_24); + + + /** + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp-finite element dof handlers. + * + * @deprecated Use the version of this function which takes a + * hp::MappingCollection instead. + */ + template + DEAL_II_DEPRECATED_EARLY_WITH_COMMENT( + "Use the version of this function which takes a hp::MappingCollection instead.") + static void estimate( const Mapping &mapping, const DoFHandler &dof, const hp::QCollection &quadrature, @@ -470,11 +499,37 @@ public: /** * Equivalent to the set of functions above, except that this one takes a - * quadrature collection for hp-finite element dof handlers. + * mapping and quadrature collection for hp-finite element dof handlers. */ template static void estimate( + const hp::MappingCollection &mapping, + const DoFHandler &dof, + const hp::QCollection &quadrature, + const std::map *> + &neumann_bc, + const ArrayView *> &solutions, + ArrayView *> &errors, + const ComponentMask &component_mask = {}, + 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 Strategy strategy = cell_diameter_over_24); + + + /** + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp-finite element dof handlers. + * + * @deprecated Use the version of this function which takes a + * hp::MappingCollection instead. + */ + template + DEAL_II_DEPRECATED_EARLY_WITH_COMMENT( + "Use the version of this function which takes a hp::MappingCollection instead.") + static void estimate( const Mapping &mapping, const DoFHandler &dof, const hp::QCollection &quadrature, @@ -714,11 +769,38 @@ public: /** * Equivalent to the set of functions above, except that this one takes a - * quadrature collection for hp-finite element dof handlers. + * mapping and quadrature collection for hp-finite element dof handlers. */ template static void estimate( + const hp::MappingCollection<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof, + const hp::QCollection<0> &quadrature, + const std::map *> + &neumann_bc, + const ReadVector &solution, + Vector &error, + const ComponentMask &component_mask = {}, + 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 Strategy strategy = cell_diameter_over_24); + + + + /** + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp-finite element dof handlers. + * + * @deprecated Use the version of this function which takes a + * hp::MappingCollection instead. + */ + template + DEAL_II_DEPRECATED_EARLY_WITH_COMMENT( + "Use the version of this function which takes a hp::MappingCollection instead.") + static void estimate( const Mapping<1, spacedim> &mapping, const DoFHandler<1, spacedim> &dof, const hp::QCollection<0> &quadrature, @@ -759,11 +841,38 @@ public: /** * Equivalent to the set of functions above, except that this one takes a - * quadrature collection for hp-finite element dof handlers. + * mapping and quadrature collection for hp-finite element dof handlers. */ template static void estimate( + const hp::MappingCollection<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof, + const hp::QCollection<0> &quadrature, + const std::map *> + &neumann_bc, + const ArrayView *> &solutions, + ArrayView *> &errors, + const ComponentMask &component_mask = {}, + 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 Strategy strategy = cell_diameter_over_24); + + + + /** + * Equivalent to the set of functions above, except that this one takes a + * quadrature collection for hp-finite element dof handlers. + * + * @deprecated Use the version of this function which takes a + * hp::MappingCollection instead. + */ + template + DEAL_II_DEPRECATED_EARLY_WITH_COMMENT( + "Use the version of this function which takes a hp::MappingCollection instead.") + static void estimate( const Mapping<1, spacedim> &mapping, const DoFHandler<1, spacedim> &dof, const hp::QCollection<0> &quadrature, diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h index 723c2d62c9..c7a3d9b7dd 100644 --- a/include/deal.II/numerics/error_estimator.templates.h +++ b/include/deal.II/numerics/error_estimator.templates.h @@ -1063,6 +1063,43 @@ KellyErrorEstimator::estimate( } +template +template +void +KellyErrorEstimator::estimate( + const hp::MappingCollection &mapping, + const DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const std::map *> + &neumann_bc, + const ReadVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id, + const Strategy strategy) +{ + // just pass on to the other function + std::vector *> solutions(1, &solution); + std::vector *> errors(1, &error); + ArrayView *> error_view = make_array_view(errors); + estimate(mapping, + dof_handler, + quadrature, + neumann_bc, + make_array_view(solutions), + error_view, + component_mask, + coefficients, + n_threads, + subdomain_id, + material_id, + strategy); +} + + template template void @@ -1081,11 +1118,13 @@ KellyErrorEstimator::estimate( const types::material_id material_id, const Strategy strategy) { + // DEPRECATED // just pass on to the other function std::vector *> solutions(1, &solution); std::vector *> errors(1, &error); ArrayView *> error_view = make_array_view(errors); - estimate(mapping, + const hp::MappingCollection mapping_collection(mapping); + estimate(mapping_collection, dof_handler, quadrature, neumann_bc, @@ -1117,7 +1156,9 @@ KellyErrorEstimator::estimate( const types::material_id material_id, const Strategy strategy) { - estimate(get_default_linear_mapping(dof_handler.get_triangulation()), + const hp::MappingCollection mapping( + get_default_linear_mapping(dof_handler.get_triangulation())); + estimate(mapping, dof_handler, quadrature, neumann_bc, @@ -1137,9 +1178,9 @@ template template void KellyErrorEstimator::estimate( - const Mapping &mapping, - const DoFHandler &dof_handler, - const hp::QCollection &face_quadratures, + const hp::MappingCollection &mapping, + const DoFHandler &dof_handler, + const hp::QCollection &face_quadratures, const std::map *> &neumann_bc, const ArrayView *> &solutions, @@ -1211,11 +1252,10 @@ KellyErrorEstimator::estimate( // all the data needed in the error estimator by each of the threads is // gathered in the following structures - const hp::MappingCollection mapping_collection(mapping); const internal::ParallelData parallel_data( dof_handler.get_fe_collection(), face_quadratures, - mapping_collection, + mapping, (!neumann_bc.empty() || (coefficients != nullptr)), solutions.size(), subdomain_id, @@ -1299,6 +1339,43 @@ KellyErrorEstimator::estimate( +template +template +void +KellyErrorEstimator::estimate( + const Mapping &mapping, + const DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const std::map *> + &neumann_bc, + const ArrayView *> &solutions, + ArrayView *> &errors, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id, + const Strategy strategy) +{ + // DEPRECATED + // just pass on to the other function + const hp::MappingCollection mapping_collection(mapping); + estimate(mapping_collection, + dof_handler, + quadrature, + neumann_bc, + solutions, + errors, + component_mask, + coefficients, + n_threads, + subdomain_id, + material_id, + strategy); +} + + + template template void @@ -1317,10 +1394,11 @@ KellyErrorEstimator::estimate( const types::material_id material_id, const Strategy strategy) { - // forward to the function with the QCollection - estimate(mapping, + const hp::MappingCollection mapping_collection(mapping); + const hp::QCollection quadrature_collection(quadrature); + estimate(mapping_collection, dof_handler, - hp::QCollection(quadrature), + quadrature_collection, neumann_bc, solutions, errors, @@ -1383,7 +1461,9 @@ KellyErrorEstimator::estimate( const types::material_id material_id, const Strategy strategy) { - estimate(get_default_linear_mapping(dof_handler.get_triangulation()), + const hp::MappingCollection mapping( + get_default_linear_mapping(dof_handler.get_triangulation())); + estimate(mapping, dof_handler, quadrature, neumann_bc, diff --git a/source/numerics/error_estimator.inst.in b/source/numerics/error_estimator.inst.in index 8741c991af..7a8fd99a99 100644 --- a/source/numerics/error_estimator.inst.in +++ b/source/numerics/error_estimator.inst.in @@ -92,6 +92,25 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; const KellyErrorEstimator::Strategy); + template void + KellyErrorEstimator::estimate< + S>( + const hp::MappingCollection &, + const DoFHandler &, + const hp::QCollection &, + const std::map *> &, + const ReadVector &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int, + const types::subdomain_id, + const types::material_id, + const KellyErrorEstimator::Strategy); + + // DEPRECATED template void KellyErrorEstimator::estimate< S>(const Mapping &, @@ -125,6 +144,25 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; const KellyErrorEstimator::Strategy); + template void + KellyErrorEstimator::estimate< + S>( + const hp::MappingCollection &, + const DoFHandler &, + const hp::QCollection &, + const std::map *> &, + const ArrayView *> &, + ArrayView *> &, + const ComponentMask &, + const Function *, + const unsigned int, + const types::subdomain_id, + const types::material_id, + const KellyErrorEstimator::Strategy); + + // DEPRECATED template void KellyErrorEstimator::estimate< S>(const Mapping &, diff --git a/source/numerics/error_estimator_1d.cc b/source/numerics/error_estimator_1d.cc index ebc93cc7f8..e3de1b2071 100644 --- a/source/numerics/error_estimator_1d.cc +++ b/source/numerics/error_estimator_1d.cc @@ -163,9 +163,9 @@ template template void KellyErrorEstimator<1, spacedim>::estimate( - const Mapping<1, spacedim> &mapping, - const DoFHandler<1, spacedim> &dof_handler, - const hp::QCollection<0> &quadrature, + const hp::MappingCollection<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof_handler, + const hp::QCollection<0> &quadrature, const std::map *> &neumann_bc, const ReadVector &solution, @@ -197,6 +197,46 @@ KellyErrorEstimator<1, spacedim>::estimate( +template +template +void +KellyErrorEstimator<1, spacedim>::estimate( + const Mapping<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof_handler, + const hp::QCollection<0> &quadrature, + const std::map *> + &neumann_bc, + const ReadVector &solution, + Vector &error, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id, + const Strategy strategy) +{ + // DEPRECATED + // just pass on to the other function + std::vector *> solutions(1, &solution); + std::vector *> errors(1, &error); + ArrayView *> error_view = make_array_view(errors); + const hp::MappingCollection<1, spacedim> mapping_collection(mapping); + estimate(mapping_collection, + dof_handler, + quadrature, + neumann_bc, + make_array_view(solutions), + error_view, + component_mask, + coefficients, + n_threads, + subdomain_id, + material_id, + strategy); +} + + + template template void @@ -215,7 +255,9 @@ KellyErrorEstimator<1, spacedim>::estimate( const Strategy strategy) { const auto reference_cell = ReferenceCells::Line; - estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(), + const hp::MappingCollection<1, spacedim> mapping( + reference_cell.template get_default_linear_mapping<1, spacedim>()); + estimate(mapping, dof_handler, quadrature, neumann_bc, @@ -249,7 +291,9 @@ KellyErrorEstimator<1, spacedim>::estimate( const Strategy strategy) { const auto reference_cell = ReferenceCells::Line; - estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(), + const hp::MappingCollection<1, spacedim> mapping( + reference_cell.template get_default_linear_mapping<1, spacedim>()); + estimate(mapping, dof_handler, quadrature, neumann_bc, @@ -269,8 +313,8 @@ template template void KellyErrorEstimator<1, spacedim>::estimate( - const Mapping<1, spacedim> &mapping, - const DoFHandler<1, spacedim> &dof_handler, + const hp::MappingCollection<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof_handler, const hp::QCollection<0> &, const std::map *> &neumann_bc, @@ -388,15 +432,13 @@ KellyErrorEstimator<1, spacedim>::estimate( const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection(); - hp::MappingCollection<1, spacedim> mapping_collection; - mapping_collection.push_back(mapping); - hp::FEValues<1, spacedim> fe_values(mapping_collection, + hp::FEValues<1, spacedim> fe_values(mapping, fe, q_collection, update_gradients); hp::FEFaceValues<1, spacedim> fe_face_values( - /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors); + /*mapping,*/ fe, q_face_collection, update_normal_vectors); // loop over all cells and do something on the cells which we're told to // work on. note that the error indicator is only a sum over the two @@ -522,6 +564,43 @@ KellyErrorEstimator<1, spacedim>::estimate( +template +template +void +KellyErrorEstimator<1, spacedim>::estimate( + const Mapping<1, spacedim> &mapping, + const DoFHandler<1, spacedim> &dof_handler, + const hp::QCollection<0> &quadrature, + const std::map *> + &neumann_bc, + const ArrayView *> &solutions, + ArrayView *> &errors, + const ComponentMask &component_mask, + const Function *coefficients, + const unsigned int n_threads, + const types::subdomain_id subdomain_id, + const types::material_id material_id, + const Strategy strategy) +{ + // DEPRECATED + // just pass on to the other function + const hp::MappingCollection<1, spacedim> mapping_collection(mapping); + estimate(mapping_collection, + dof_handler, + quadrature, + neumann_bc, + solutions, + errors, + component_mask, + coefficients, + n_threads, + subdomain_id, + material_id, + strategy); +} + + + template template void @@ -540,8 +619,9 @@ KellyErrorEstimator<1, spacedim>::estimate( const types::material_id material_id, const Strategy strategy) { - const hp::QCollection<0> quadrature_collection(quadrature); - estimate(mapping, + const hp::MappingCollection<1, spacedim> mapping_collection(mapping); + const hp::QCollection<0> quadrature_collection(quadrature); + estimate(mapping_collection, dof_handler, quadrature_collection, neumann_bc, diff --git a/source/numerics/error_estimator_1d.inst.in b/source/numerics/error_estimator_1d.inst.in index a43b76fb77..e1f967ee2b 100644 --- a/source/numerics/error_estimator_1d.inst.in +++ b/source/numerics/error_estimator_1d.inst.in @@ -80,6 +80,24 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; const types::material_id, const Strategy); + template void + KellyErrorEstimator::estimate< + S>( + const hp::MappingCollection &, + const DoFHandler &, + const hp::QCollection &, + const std::map *> &, + const ReadVector &, + Vector &, + const ComponentMask &, + const Function *, + const unsigned int, + const types::subdomain_id, + const types::material_id, + const Strategy); + + // DEPRECATED template void KellyErrorEstimator::estimate< S>(const Mapping &, @@ -111,6 +129,24 @@ for (S : REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; const types::material_id, const Strategy); + template void + KellyErrorEstimator::estimate< + S>( + const hp::MappingCollection &, + const DoFHandler &, + const hp::QCollection &, + const std::map *> &, + const ArrayView *> &, + ArrayView *> &, + const ComponentMask &, + const Function *, + const unsigned int, + const types::subdomain_id, + const types::material_id, + const Strategy); + + // DEPRECATED template void KellyErrorEstimator::estimate< S>(const Mapping &, diff --git a/tests/hp/error_estimator_01.cc b/tests/hp/error_estimator_01.cc new file mode 100644 index 0000000000..2f54e19790 --- /dev/null +++ b/tests/hp/error_estimator_01.cc @@ -0,0 +1,110 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2025 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// Test that we pick the correct mapping in KellyErrorEstimator +// for hp applications. + + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +void +test() +{ + // We investigate the effect of different mappings in an hp-context. + // Thus we only feed the MappingCollection with different mappings, + // and keep the FEs and quadratures the same. + hp::FECollection fes(FE_Q(1), FE_Q(1)); + hp::QCollection quads_face(QGauss(2), QGauss(2)); + hp::MappingCollection mappings(MappingQ(1), MappingQ(2)); + + // Set up hyper_ball and assign different active FE indices + // to inner and outer cells. + Triangulation tria; + DoFHandler dofh(tria); + + GridGenerator::hyper_ball(tria); + for (const auto &cell : dofh.active_cell_iterators()) + if (cell->center() != Point()) + cell->set_active_fe_index(1); + tria.refine_global(1); + + dofh.distribute_dofs(fes); + + // Interpolate a transcendental function and estimate the error via Kelly + // with different mappings. + Functions::CosineFunction function; + Vector interpolation(dofh.n_dofs()); + + Vector error(tria.n_active_cells()); + + // MappingQ1 + hp::MappingCollection mapping_q1(mappings[0]); + VectorTools::interpolate(mapping_q1, dofh, function, interpolation); + KellyErrorEstimator::estimate( + mapping_q1, dofh, quads_face, {}, interpolation, error); + deallog << "MappingQ1 : global error estimate=" << error.l2_norm() + << std::endl; + + // MappingQ2 + hp::MappingCollection mapping_q2(mappings[1]); + VectorTools::interpolate(mapping_q2, dofh, function, interpolation); + KellyErrorEstimator::estimate( + mapping_q2, dofh, quads_face, {}, interpolation, error); + deallog << "MappingQ2 : global error estimate=" << error.l2_norm() + << std::endl; + + // MappingQ1&2 + VectorTools::interpolate(mappings, dofh, function, interpolation); + KellyErrorEstimator::estimate( + mappings, dofh, quads_face, {}, interpolation, error); + deallog << "MappingQ1&2: global error estimate=" << error.l2_norm() + << std::endl; +} + + +int +main() +{ + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/hp/error_estimator_01.output b/tests/hp/error_estimator_01.output new file mode 100644 index 0000000000..30d9bae57e --- /dev/null +++ b/tests/hp/error_estimator_01.output @@ -0,0 +1,7 @@ + +DEAL:2d::MappingQ1 : global error estimate=0.585534 +DEAL:2d::MappingQ2 : global error estimate=0.465513 +DEAL:2d::MappingQ1&2: global error estimate=0.465513 +DEAL:3d::MappingQ1 : global error estimate=0.922872 +DEAL:3d::MappingQ2 : global error estimate=0.603987 +DEAL:3d::MappingQ1&2: global error estimate=0.603987