From aae8646ee614a69e9a928dd6adb8277489710de2 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 27 Apr 2024 09:02:38 +0200 Subject: [PATCH] NonMatching::FEValues: accept CellAccessor --- include/deal.II/non_matching/fe_values.h | 59 +++++++++++- source/non_matching/fe_values.cc | 89 ++++++++++++++++--- source/non_matching/fe_values.inst.in | 4 +- tests/non_matching/fe_values.cc | 20 +++-- .../fe_values.with_lapack=on.output | 54 +++++++++++ 5 files changed, 206 insertions(+), 20 deletions(-) diff --git a/include/deal.II/non_matching/fe_values.h b/include/deal.II/non_matching/fe_values.h index 88588a7f70..6a014fba8a 100644 --- a/include/deal.II/non_matching/fe_values.h +++ b/include/deal.II/non_matching/fe_values.h @@ -216,14 +216,57 @@ namespace NonMatching * Reinitialize the various FEValues-like objects for the 3 different * regions of the cell. After calling this function an FEValues-like object * can be retrieved for each region using the functions - * get_inside_fe_values(), - * get_outside_fe_values(), or + * get_inside_fe_values(), get_outside_fe_values(), or * get_surface_fe_values(). + * + * If the @p q_index argument is left at its default value, then we use + * that quadrature formula within the hp::QCollection passed to the + * constructor of this class with index given by + * cell-@>active_fe_index(), i.e. the same index as that of + * the finite element. In this case, there should be a corresponding + * quadrature formula for each finite element in the hp::FECollection. As + * a special case, if the quadrature collection contains only a single + * element (a frequent case if one wants to use the same quadrature object + * for all finite elements in an hp-discretization, even if that may not + * be the most efficient), then this single quadrature is used unless a + * different value for this argument is specified. On the other hand, if a + * value is given for this argument, it overrides the choice of + * cell-@>active_fe_index() or the choice for the single + * quadrature. + * + * If the @p mapping_index argument is left at its default value, then we + * use that mapping object within the hp::MappingCollection passed to the + * constructor of this class with index given by + * cell-@>active_fe_index(), i.e. the same index as that of + * the finite element. As above, if the mapping collection contains only a + * single element (a frequent case if one wants to use a $Q_1$ mapping for + * all finite elements in an hp-discretization), then this single mapping + * is used unless a different value for this argument is specified. */ template void reinit( - const TriaIterator> &cell); + const TriaIterator> &cell, + const unsigned int q_index = numbers::invalid_unsigned_int, + const unsigned int mapping_index = numbers::invalid_unsigned_int); + + /** + * Like the previous function, but for non-DoFHandler iterators. The reason + * this function exists is so that one can use NonMatching::FEValues for + * Triangulation objects too. + * + * Since cell-@>active_fe_index() doesn't make sense for + * triangulation iterators, this function chooses the zero-th finite + * element, mapping, and quadrature object from the relevant constructions + * passed to the constructor of this object. The only exception is if you + * specify a value different from the default value for any of these last + * three arguments. + */ + void + reinit(const TriaIterator> &cell, + const unsigned int q_index = numbers::invalid_unsigned_int, + const unsigned int mapping_index = numbers::invalid_unsigned_int, + const unsigned int fe_index = numbers::invalid_unsigned_int); /** * Return an dealii::FEValues object reinitialized with a quadrature for the @@ -258,6 +301,16 @@ namespace NonMatching get_surface_fe_values() const; private: + /** + * Internal function called by the reinit() functions. + */ + template + void + reinit_internal(const CellIteratorType &cell, + const unsigned int q_index, + const unsigned int mapping_index, + const unsigned int fe_index); + /** * Do work common to the constructors. The incoming QCollection should be * quadratures integrating over $[0, 1]^{dim}$. These will be used on the diff --git a/source/non_matching/fe_values.cc b/source/non_matching/fe_values.cc index df9c45d2fe..56675f15ce 100644 --- a/source/non_matching/fe_values.cc +++ b/source/non_matching/fe_values.cc @@ -148,10 +148,72 @@ namespace NonMatching template void FEValues::reinit( - const TriaIterator> &cell) + const TriaIterator> &cell, + const unsigned int q_index, + const unsigned int mapping_index) + { + this->reinit_internal(cell, + q_index, + mapping_index, + cell->active_fe_index()); + } + + + + template + void + FEValues::reinit(const TriaIterator> &cell, + const unsigned int q_index, + const unsigned int mapping_index, + const unsigned int fe_index) + { + this->reinit_internal(cell, q_index, mapping_index, fe_index); + } + + + + template + template + void + FEValues::reinit_internal(const CellIteratorType &cell, + const unsigned int q_index_in, + const unsigned int mapping_index_in, + const unsigned int fe_index_in) { current_cell_location = mesh_classifier->location_to_level_set(cell); - active_fe_index = cell->active_fe_index(); + + if (fe_index_in == numbers::invalid_unsigned_int) + this->active_fe_index = 0; + else + this->active_fe_index = fe_index_in; + + unsigned int mapping_index = mapping_index_in; + unsigned int q_index = q_index_in; + unsigned int q_index_1D = q_index_in; + + if (mapping_index == numbers::invalid_unsigned_int) + { + if (mapping_collection->size() > 1) + mapping_index = active_fe_index; + else + mapping_index = 0; + } + + if (q_index == numbers::invalid_unsigned_int) + { + if (fe_values_inside_full_quadrature.size() > 1) + q_index = active_fe_index; + else + q_index = 0; + } + + if (q_index_1D == numbers::invalid_unsigned_int) + { + if (q_collection_1D.size() > 1) + q_index_1D = active_fe_index; + else + q_index_1D = 0; + } // These objects were created with a quadrature based on the previous cell // and are thus no longer valid. @@ -163,22 +225,29 @@ namespace NonMatching { case LocationToLevelSet::inside: { - fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell); + Assert((active_fe_index == mapping_index) || + ((mapping_collection->size() == 1) && + (mapping_index == 0)), + ExcNotImplemented()); + Assert(active_fe_index == q_index, ExcNotImplemented()); + + fe_values_inside_full_quadrature.at(q_index)->reinit(cell); break; } case LocationToLevelSet::outside: { - fe_values_outside_full_quadrature.at(active_fe_index)->reinit(cell); + Assert((active_fe_index == mapping_index) || + ((mapping_collection->size() == 1) && + (mapping_index == 0)), + ExcNotImplemented()); + Assert(active_fe_index == q_index, ExcNotImplemented()); + + fe_values_outside_full_quadrature.at(q_index)->reinit(cell); break; } case LocationToLevelSet::intersected: { - const unsigned int mapping_index = - mapping_collection->size() > 1 ? active_fe_index : 0; - - const unsigned int q1D_index = - q_collection_1D.size() > 1 ? active_fe_index : 0; - quadrature_generator.set_1D_quadrature(q1D_index); + quadrature_generator.set_1D_quadrature(q_index_1D); quadrature_generator.generate(cell); const Quadrature &inside_quadrature = diff --git a/source/non_matching/fe_values.inst.in b/source/non_matching/fe_values.inst.in index 91c0172018..31a3b1a3f5 100644 --- a/source/non_matching/fe_values.inst.in +++ b/source/non_matching/fe_values.inst.in @@ -69,7 +69,9 @@ for (deal_II_dimension : DIMENSIONS; lda : BOOL) { template void FEValues::reinit( const TriaIterator< - DoFCellAccessor> &); + DoFCellAccessor> &, + const unsigned int, + const unsigned int); template void FEInterfaceValues::do_reinit( const TriaIterator< diff --git a/tests/non_matching/fe_values.cc b/tests/non_matching/fe_values.cc index b64d2289b8..825a59667b 100644 --- a/tests/non_matching/fe_values.cc +++ b/tests/non_matching/fe_values.cc @@ -88,9 +88,10 @@ private: void setup_discrete_level_set(); + template void - test_fe_values_reinitializes_correctly( - NonMatching::FEValues &fe_values) const; + test_fe_values_reinitializes_correctly(NonMatching::FEValues &fe_values, + IteratorType cell) const; Triangulation triangulation; hp::FECollection fe_collection; @@ -139,7 +140,10 @@ Test::run() mesh_classifier, dof_handler, level_set); - test_fe_values_reinitializes_correctly(fe_values); + test_fe_values_reinitializes_correctly(fe_values, + triangulation.begin_active()); + test_fe_values_reinitializes_correctly(fe_values, + dof_handler.begin_active()); } { // Test with the "more advanced" constructor. @@ -151,7 +155,10 @@ Test::run() mesh_classifier, dof_handler, level_set); - test_fe_values_reinitializes_correctly(fe_values); + test_fe_values_reinitializes_correctly(fe_values, + triangulation.begin_active()); + test_fe_values_reinitializes_correctly(fe_values, + dof_handler.begin_active()); } } @@ -197,13 +204,14 @@ Test::setup_discrete_level_set() template +template void Test::test_fe_values_reinitializes_correctly( - NonMatching::FEValues &fe_values) const + NonMatching::FEValues &fe_values, + IteratorType cell) const { // The first is inside so only the inside FEValues object should be // initialized. - auto cell = dof_handler.begin_active(); fe_values.reinit(cell); assert_cells_are_the_same(cell, diff --git a/tests/non_matching/fe_values.with_lapack=on.output b/tests/non_matching/fe_values.with_lapack=on.output index bd38077635..ae5896955b 100644 --- a/tests/non_matching/fe_values.with_lapack=on.output +++ b/tests/non_matching/fe_values.with_lapack=on.output @@ -18,6 +18,24 @@ DEAL::OK DEAL::OK DEAL::OK DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK DEAL:: DEAL::dim = 2 DEAL::OK @@ -38,6 +56,24 @@ DEAL::OK DEAL::OK DEAL::OK DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK DEAL:: DEAL::dim = 3 DEAL::OK @@ -58,4 +94,22 @@ DEAL::OK DEAL::OK DEAL::OK DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK DEAL:: -- 2.39.5