From aa4a28f4ecd3d9487cea76d8f8911a42959f0d37 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Fri, 29 Apr 2022 10:29:56 +0200 Subject: [PATCH] Add class NonMatching::FEInterfaceValues This class works the same way as NonMatching::FEValues but for assembling interface terms on faces. This is done using the DisceteFaceQuadratureGenerator class to create quadrature rules on the intersected faces and then creating dealii::FEInterfaceValues objects with these quadratures. --- .../non_matching_fe_interface_values.svg | 591 ++++++++++++++++++ include/deal.II/non_matching/fe_values.h | 374 +++++++++++ source/non_matching/fe_values.cc | 211 +++++++ source/non_matching/fe_values.inst.in | 31 +- tests/non_matching/fe_interface_values.cc | 227 +++++++ tests/non_matching/fe_interface_values.output | 73 +++ ...fe_interface_values_non_standard_meshes.cc | 206 ++++++ ...nterface_values_non_standard_meshes.output | 193 ++++++ 8 files changed, 1905 insertions(+), 1 deletion(-) create mode 100644 doc/doxygen/images/non_matching_fe_interface_values.svg create mode 100644 tests/non_matching/fe_interface_values.cc create mode 100644 tests/non_matching/fe_interface_values.output create mode 100644 tests/non_matching/fe_interface_values_non_standard_meshes.cc create mode 100644 tests/non_matching/fe_interface_values_non_standard_meshes.output diff --git a/doc/doxygen/images/non_matching_fe_interface_values.svg b/doc/doxygen/images/non_matching_fe_interface_values.svg new file mode 100644 index 0000000000..3c7de307fa --- /dev/null +++ b/doc/doxygen/images/non_matching_fe_interface_values.svg @@ -0,0 +1,591 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/deal.II/non_matching/fe_values.h b/include/deal.II/non_matching/fe_values.h index b31ba20717..8fda2693ce 100644 --- a/include/deal.II/non_matching/fe_values.h +++ b/include/deal.II/non_matching/fe_values.h @@ -22,6 +22,7 @@ #include +#include #include #include @@ -368,6 +369,379 @@ namespace NonMatching DiscreteQuadratureGenerator quadrature_generator; }; + + + /** + * This class is intended to facilitate assembling interface terms on faces in + * immersed (in the sense of cut) finite element methods. These types of terms + * occur mainly in cut discontinuous Galerkin methods. This class works + * analogously to NonMatching::FEValues. The domain is assumed to be described + * by a level set function, + * $\psi : \mathbb{R}^{dim} \to \mathbb{R}$, + * and this class assumes that we want to integrate over two different regions + * of each face, $F$: + * @f[ + * N = \{x \in F : \psi(x) < 0 \}, \\ + * P = \{x \in F : \psi(x) > 0 \}, + * @f] + * which we as before refer to as the "inside" and "outside" regions of the + * face. + * @image html non_matching_fe_interface_values.svg + * When calling one of the reinit functions of this class, this class + * will create quadrature rules over these regions in the background and set + * up dealii::FEInterfaceValues objects with these quadratures. These objects + * can then be accessed through one of the functions: get_inside_fe_values() + * or get_outside_fe_values(). For the same reason as described in + * NonMatching::FEValues, the returned the FEInterfaceValues objects that + * get_inside/outside_fe_values() returns is wrapped in a + * std_cxx17::optional. Assembling using this class would then typically be + * done in the following way: + * + * @code + * NonMatching::FEInterfaceValues fe_interface_values(...) + * + * for (const auto &cell : dof_handler.active_cell_iterators()) + * { + * for (unsigned int face_index : cell->face_indices()) + * { + * if (cell->at_boundary(face_index)) + * fe_interface_values.reinit(cell, face_index); + * else + * fe_interface_values.reinit(cell, + * face_index, + * subface_index, + * cell->neighbor(face_index), + * cell->neighbor_of_neighbor(face_index), + * neighbor_subface_index); + * + * const std_cxx17::optional> + * &inside_fe_values = fe_interface_values.get_inside_fe_values(); + * + * // Check if the returned optional contains a value + * if (inside_fe_values) + * { + * // Assemble locally + * } + * } + * } + * @endcode + * + * To reduce the amount of work, the reinit() function of this class uses the + * MeshClassifier passed to the constructor to check how the incoming cell + * relates to the level set function. The immersed quadrature rules are only + * generated if the cell is intersected. If the cell is completely inside or + * outside, it returns a cached FEInterfaceValues object created with a + * quadrature over the reference cell: $[0, 1]^{dim-1}$. + */ + template + class FEInterfaceValues + { + public: + using AdditionalData = + typename FaceQuadratureGenerator::AdditionalData; + + /** + * Constructor. + * + * @param fe_collection Collection of FiniteElements to be used. + * @param quadrature 1-dimensional quadrature rule used to generate the + * immersed quadrature rules. See the QuadratureGenerator class. On the non + * intersected faces a tensor product of this quadrature will be used. + * @param mesh_classifier Object used to determine when the immersed + * quadrature rules need to be generated. + * @param region_update_flags Struct storing UpdateFlags for the + * inside/outside region of the cell. + * @param dof_handler The DoFHandler associated with the discrete level set + * function. + * @param level_set The degrees of freedom of the discrete level set function. + * @param additional_data Additional data passed to the QuadratureGenerator + * class. + * + * @note Pointers to @p fe_collection, @p mesh_classifier, @p dof_handler, + * and @p level_set are stored internally, so these need to have a longer life + * span than the instance of this class. + */ + template + FEInterfaceValues(const hp::FECollection &fe_collection, + const Quadrature<1> & quadrature, + const RegionUpdateFlags region_update_flags, + const MeshClassifier & mesh_classifier, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData &additional_data = AdditionalData()); + + /** + * Constructor. + * + * @param mapping_collection Collection of Mappings to be used. + * @param fe_collection Collection of FiniteElements to be used. + * @param q_collection Collection of Quadrature rules over $[0, 1]^{dim-1}$ + * that should be used when a face is not intersected and we do not need to + * generate immersed quadrature rules. + * @param q_collection_1D Collection of 1-dimensional quadrature rules used + * to generate the immersed quadrature rules. See the QuadratureGenerator + * class. + * @param mesh_classifier Object used to determine when the immersed + * quadrature rules need to be generated. + * @param region_update_flags Struct storing UpdateFlags for the + * inside/outside region of the cell. + * @param dof_handler The DoFHandler associated with the discrete level set + * function. + * @param level_set The degrees of freedom of the discrete level set function. + * @param additional_data Additional data passed to the QuadratureGenerator + * class. + * + * @note Pointers to @p mapping_collection, @p fe_collection, + * @p mesh_classifier, @p dof_handler, and @p level_set are stored + * internally, so these need to have a longer life span than the instance of + * this class. + */ + template + FEInterfaceValues(const hp::MappingCollection &mapping_collection, + const hp::FECollection & fe_collection, + const hp::QCollection & q_collection, + const hp::QCollection<1> & q_collection_1D, + const RegionUpdateFlags region_update_flags, + const MeshClassifier & mesh_classifier, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData &additional_data = AdditionalData()); + + /** + * Reinitialize on the shared face between two neighboring cells. + * After calling this function, an FEInterfaceValues object can be retrieved + * for each region using the functions get_inside_fe_values() or + * get_outside_fe_values(). + * + * @note Currently @p cell and @p cell_neighbor need to have the same + * finite element associated with them. + * + * @note Specifying sub_face_no/sub_face_no_neighbor is currently not + * implemented, the passed values of these must equal + * numbers::invalid_unsigned_int. + */ + template + void + reinit(const CellIteratorType & cell, + const unsigned int face_no, + const unsigned int sub_face_no, + const CellNeighborIteratorType &cell_neighbor, + const unsigned int face_no_neighbor, + const unsigned int sub_face_no_neighbor); + + + /** + * Reinitialize on a face on the boundary of the domain, + * when there is no neighbor and a single face @p face_no of + * the @p cell. After calling this function, an + * FEInterfaceValues object can be retrieved for each region using the + * functions get_inside_fe_values() or get_outside_fe_values(). + */ + template + void + reinit(const CellIteratorType &cell, const unsigned int face_no); + + /** + * Return an dealii::FEInterfaceValues object reinitialized with a + * quadrature for the inside region of the cell: + * $\{x \in K : \psi(x) < 0 \}$. + * + * @note If the quadrature rule over the region is empty, e.g. because the + * cell is completely located in the outside domain, the returned + * optional will not contain a value. + */ + const std_cxx17::optional> & + get_inside_fe_values() const; + + /** + * Return an dealii::FEInterfaceValues object reinitialized with a + * quadrature for the outside region of the cell: + * $\{x \in K : \psi(x) > 0 \}$. + * + * @note If the quadrature rule over the region is empty, e.g. because the + * cell is completely located in the inside domain, the returned + * optional will not contain a value. + */ + const std_cxx17::optional> & + get_outside_fe_values() const; + + private: + /** + * Do work common to the constructors. The incoming QCollection should be + * quadratures integrating over $[0, 1]^{dim-1}$. These will be used on the + * non-intersected cells. + */ + void + initialize(const hp::QCollection &q_collection); + + /** + * Do work that is common two the two reinit functions of the class. + * @p call_reinit is a std::function that describes how we should call + * reinit on a single dealii::FEInterfaceValues object, which is what + * differs between the two reinit functions. + */ + template + void + do_reinit( + const TriaIterator> &cell, + const unsigned int face_no, + const std::function &)> &call_reinit); + + /** + * A pointer to the collection of mappings to be used. + */ + const SmartPointer> mapping_collection; + + /** + * A pointer to the collection of finite elements to be used. + */ + const SmartPointer> fe_collection; + + /** + * Collection of 1-dimensional quadrature rules that are used by + * QuadratureGenerator as base for generating the immersed quadrature rules. + */ + const hp::QCollection<1> q_collection_1D; + + /** + * Location of the last face that reinit was called with. + */ + LocationToLevelSet current_face_location; + + /** + * Active fe index of the last cell that reinit was called with. + */ + unsigned int active_fe_index; + + /** + * The update flags passed to the constructor. + */ + const RegionUpdateFlags region_update_flags; + + /** + * Pointer to the MeshClassifier passed to the constructor. + */ + const SmartPointer> mesh_classifier; + + /** + * For each element in the FECollection passed to the constructor, + * this object contains an dealii::FEInterfaceValues object created with a + * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and + * UpdateFlags for the inside region. Thus, these optionals should always + * contain a value. + * + * When LocationToLevelSet of the cell is INSIDE (and we do not need + * to generate an immersed quadrature), we return the + * dealii::FEInterfaceValues object in this container corresponding to the + * cell's active_fe_index. + * + * This container is a std::deque, which is compatible with the + * `FEInterfaceValues` class that does not have a copy-constructor. + */ + std::deque>> + fe_values_inside_full_quadrature; + + /** + * For each element in the FECollection passed to the constructor, + * this object contains an dealii::FEInterfaceValues object created with a + * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and + * UpdateFlags for the outside region. Thus, these optionals should always + * contain a value. + * + * When LocationToLevelSet of the cell is OUTSIDE (and we do not need + * to generate an immersed quadrature), we return the dealii::FEValues + * object in this container corresponding to the cell's active_fe_index. + * + * This container is a std::deque, which is compatible with the + * `FEInterfaceValues` class that does not have a copy-constructor. + */ + std::deque>> + fe_values_outside_full_quadrature; + + /** + * FEInterfaceValues object created with a quadrature rule integrating over + * the inside region, $\{x \in B : \psi(x) < 0 \}$, that was generated in + * the last call to reinit(..). If the cell in the last call was not + * intersected or if 0 quadrature points were generated, this optional will + * not contain a value. + */ + std_cxx17::optional> fe_values_inside; + + /** + * FEInterfaceValues object created with a quadrature rule integrating over + * the outside region, $\{x \in B : \psi(x) > 0 \}$, that was generated in + * the last call to reinit(..). If the cell in the last call was not + * intersected or if 0 quadrature points were generated, this optional will + * not contain a value. + */ + std_cxx17::optional> fe_values_outside; + + /** + * Object that generates the immersed quadrature rules. + */ + DiscreteFaceQuadratureGenerator face_quadrature_generator; + }; + + +#ifndef DOXYGEN + + /*---------------------- Inline functions ---------------------*/ + + template + template + inline void + FEInterfaceValues::reinit(const CellIteratorType &cell, + const unsigned int face_no) + { + // Lambda describing how we should call reinit on a single + // dealii::FEInterfaceValues object. + const auto reinit_operation = + [&cell, face_no](dealii::FEInterfaceValues &fe_interface_values) { + fe_interface_values.reinit(cell, face_no); + }; + + do_reinit(cell, face_no, reinit_operation); + } + + + + template + template + inline void + FEInterfaceValues::reinit(const CellIteratorType & cell, + const unsigned int face_no, + const unsigned int sub_face_no, + const CellNeighborIteratorType &cell_neighbor, + const unsigned int face_no_neighbor, + const unsigned int sub_face_no_neighbor) + { + Assert(sub_face_no == numbers::invalid_unsigned_int, ExcNotImplemented()); + Assert(sub_face_no_neighbor == numbers::invalid_unsigned_int, + ExcNotImplemented()); + + // Lambda describing how we should call reinit on a single + // dealii::FEInterfaceValues object. + const auto reinit_operation = + [&cell, + face_no, + sub_face_no, + &cell_neighbor, + face_no_neighbor, + sub_face_no_neighbor]( + dealii::FEInterfaceValues &fe_interface_values) { + fe_interface_values.reinit(cell, + face_no, + sub_face_no, + cell_neighbor, + face_no_neighbor, + sub_face_no_neighbor); + }; + + do_reinit(cell, face_no, reinit_operation); + } + +#endif + } // namespace NonMatching DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/fe_values.cc b/source/non_matching/fe_values.cc index 75d82d2526..fcbb399120 100644 --- a/source/non_matching/fe_values.cc +++ b/source/non_matching/fe_values.cc @@ -268,6 +268,217 @@ namespace NonMatching } + + template + template + FEInterfaceValues::FEInterfaceValues( + const hp::FECollection &fe_collection, + const Quadrature<1> & quadrature, + const RegionUpdateFlags region_update_flags, + const MeshClassifier & mesh_classifier, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData & additional_data) + : mapping_collection(&dealii::hp::StaticMappingQ1::mapping_collection) + , fe_collection(&fe_collection) + , q_collection_1D(quadrature) + , region_update_flags(region_update_flags) + , mesh_classifier(&mesh_classifier) + , face_quadrature_generator(q_collection_1D, + dof_handler, + level_set, + additional_data) + { + // Tensor products of each quadrature in q_collection_1D. Used on the + // non-intersected cells. + hp::QCollection q_collection; + for (unsigned int i = 0; i < q_collection_1D.size(); ++i) + q_collection.push_back(Quadrature(q_collection_1D[i])); + + initialize(q_collection); + } + + + + template + template + FEInterfaceValues::FEInterfaceValues( + const hp::MappingCollection &mapping_collection, + const hp::FECollection & fe_collection, + const hp::QCollection & q_collection, + const hp::QCollection<1> & q_collection_1D, + const RegionUpdateFlags region_update_flags, + const MeshClassifier & mesh_classifier, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData & additional_data) + : mapping_collection(&mapping_collection) + , fe_collection(&fe_collection) + , q_collection_1D(q_collection_1D) + , region_update_flags(region_update_flags) + , mesh_classifier(&mesh_classifier) + , face_quadrature_generator(q_collection_1D, + dof_handler, + level_set, + additional_data) + { + initialize(q_collection); + } + + + + template + void + FEInterfaceValues::initialize( + const hp::QCollection &q_collection) + { + current_face_location = LocationToLevelSet::unassigned; + active_fe_index = numbers::invalid_unsigned_int; + + Assert(fe_collection->size() > 0, + ExcMessage("Incoming hp::FECollection can not be empty.")); + Assert( + mapping_collection->size() == fe_collection->size() || + mapping_collection->size() == 1, + ExcMessage( + "Size of hp::MappingCollection must be the same as hp::FECollection or 1.")); + Assert( + q_collection.size() == fe_collection->size() || q_collection.size() == 1, + ExcMessage( + "Size of hp::QCollection must be the same as hp::FECollection or 1.")); + Assert( + q_collection_1D.size() == fe_collection->size() || + q_collection_1D.size() == 1, + ExcMessage( + "Size of hp::QCollection<1> must be the same as hp::FECollection or 1.")); + + // For each element in fe_collection, create dealii::FEInterfaceValues + // objects to use on the non-intersected cells. + fe_values_inside_full_quadrature.resize(fe_collection->size()); + fe_values_outside_full_quadrature.resize(fe_collection->size()); + for (unsigned int fe_index = 0; fe_index < fe_collection->size(); + ++fe_index) + { + const unsigned int mapping_index = + mapping_collection->size() > 1 ? fe_index : 0; + const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0; + + fe_values_inside_full_quadrature[fe_index].emplace( + (*mapping_collection)[mapping_index], + (*fe_collection)[fe_index], + q_collection[q_index], + region_update_flags.inside); + fe_values_outside_full_quadrature[fe_index].emplace( + (*mapping_collection)[mapping_index], + (*fe_collection)[fe_index], + q_collection[q_index], + region_update_flags.outside); + } + } + + + + template + template + void + FEInterfaceValues::do_reinit( + const TriaIterator> &cell, + const unsigned int face_no, + const std::function &)> &call_reinit) + { + current_face_location = + mesh_classifier->location_to_level_set(cell, face_no); + active_fe_index = cell->active_fe_index(); + + // These objects were created with a quadrature based on the previous cell + // and are thus no longer valid. + fe_values_inside.reset(); + fe_values_outside.reset(); + + switch (current_face_location) + { + case LocationToLevelSet::inside: + { + call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index)); + break; + } + case LocationToLevelSet::outside: + { + call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index)); + 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; + + face_quadrature_generator.set_1D_quadrature(q1D_index); + face_quadrature_generator.generate(cell, face_no); + + const Quadrature &inside_quadrature = + face_quadrature_generator.get_inside_quadrature(); + const Quadrature &outside_quadrature = + face_quadrature_generator.get_outside_quadrature(); + + // Even if a cell is formally intersected the number of created + // quadrature points can be 0. Avoid creating an FEInterfaceValues + // object if that is the case. + if (inside_quadrature.size() > 0) + { + fe_values_inside.emplace((*mapping_collection)[mapping_index], + (*fe_collection)[active_fe_index], + inside_quadrature, + region_update_flags.inside); + + call_reinit(*fe_values_inside); + } + + if (outside_quadrature.size() > 0) + { + fe_values_outside.emplace((*mapping_collection)[mapping_index], + (*fe_collection)[active_fe_index], + outside_quadrature, + region_update_flags.outside); + + call_reinit(*fe_values_outside); + } + break; + } + default: + { + Assert(false, ExcInternalError()); + break; + } + } + } + + + + template + const std_cxx17::optional> & + FEInterfaceValues::get_inside_fe_values() const + { + if (current_face_location == LocationToLevelSet::inside) + return fe_values_inside_full_quadrature.at(active_fe_index); + else + return fe_values_inside; + } + + + + template + const std_cxx17::optional> & + FEInterfaceValues::get_outside_fe_values() const + { + if (current_face_location == LocationToLevelSet::outside) + return fe_values_outside_full_quadrature.at(active_fe_index); + else + return fe_values_outside; + } + + #include "fe_values.inst" } // namespace NonMatching diff --git a/source/non_matching/fe_values.inst.in b/source/non_matching/fe_values.inst.in index e9cdf571c4..e8bccf8e59 100644 --- a/source/non_matching/fe_values.inst.in +++ b/source/non_matching/fe_values.inst.in @@ -17,6 +17,8 @@ for (deal_II_dimension : DIMENSIONS) { template class FEValues; + + template class FEInterfaceValues; } for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) @@ -41,12 +43,39 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) const DoFHandler &, const VEC &, const typename FEValues::AdditionalData &); + + template FEInterfaceValues::FEInterfaceValues( + const hp::FECollection &, + const Quadrature<1> &, + const RegionUpdateFlags, + const MeshClassifier &, + const DoFHandler &, + const VEC &, + const FEInterfaceValues::AdditionalData &); + + template FEInterfaceValues::FEInterfaceValues( + const hp::MappingCollection &, + const hp::FECollection &, + const hp::QCollection &, + const hp::QCollection<1> &, + const RegionUpdateFlags, + const MeshClassifier &, + const DoFHandler &, + const VEC &, + const AdditionalData &); } -// Template reinit function +// Template reinit functions for (deal_II_dimension : DIMENSIONS; lda : BOOL) { template void FEValues::reinit( const TriaIterator< DoFCellAccessor> &); + + template void FEInterfaceValues::do_reinit( + const TriaIterator< + DoFCellAccessor> &, + const unsigned int, + const std::function &)> + &); } diff --git a/tests/non_matching/fe_interface_values.cc b/tests/non_matching/fe_interface_values.cc new file mode 100644 index 0000000000..6677f119bd --- /dev/null +++ b/tests/non_matching/fe_interface_values.cc @@ -0,0 +1,227 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +// Set up a triangulation with 2 elements in a row: |-0-|-1-|, over +// [0,2]x[0,1]^{dim-1}, and a level set function with a cut in the plane +// $x_[dim-1] = 0.5$, such that it cuts through the middle of both cells in 2D +// and 3D and thought the middle of cell 0 in 1D. Test that we can create an +// object of type NonMatching::FEInterfaceValues. Call reinit on all the faces +// of cell 0 and check that the dealii::FEInterfaceValues objects we get from +// get_inside/outside_fe_values are set up correctly. That is, the +// std::optionals that should contain a value do and the ones that should not +// contain a value do not. +template +class Test +{ +public: + Test(); + + void + run(); + +private: + void + setup_mesh(); + + void + setup_discrete_level_set(); + + void + print_which_optionals_have_values_on_cell_0( + NonMatching::FEInterfaceValues &fe_values); + + Triangulation triangulation; + hp::FECollection fe_collection; + DoFHandler dof_handler; + + hp::MappingCollection mapping_collection; + hp::QCollection q_collection; + hp::QCollection<1> q_collection1D; + const NonMatching::RegionUpdateFlags region_update_flags; + + Vector level_set; + NonMatching::MeshClassifier mesh_classifier; +}; + + + +template +Test::Test() + : dof_handler(triangulation) + , mesh_classifier(dof_handler, level_set) +{ + fe_collection.push_back(FE_Q(1)); + mapping_collection.push_back(MappingCartesian()); + + const unsigned int n_quadrature_points = 1; + q_collection.push_back(QGauss(n_quadrature_points)); + q_collection1D.push_back(QGauss<1>(n_quadrature_points)); +} + + + +template +void +Test::run() +{ + setup_mesh(); + setup_discrete_level_set(); + mesh_classifier.reclassify(); + + // Test to create an object with each constructor. + { + deallog << "Simple constructor:" << std::endl; + NonMatching::FEInterfaceValues fe_values(fe_collection, + q_collection1D[0], + region_update_flags, + mesh_classifier, + dof_handler, + level_set); + print_which_optionals_have_values_on_cell_0(fe_values); + } + { + deallog << "Advanced constructor:" << std::endl; + NonMatching::FEInterfaceValues fe_values(mapping_collection, + fe_collection, + q_collection, + q_collection1D, + region_update_flags, + mesh_classifier, + dof_handler, + level_set); + print_which_optionals_have_values_on_cell_0(fe_values); + } +} + + + +template +void +Test::setup_mesh() +{ + const Point lower_left; + Point upper_right; + + std::vector repetitions; + upper_right[0] = 2; + repetitions.push_back(2); + for (unsigned int d = 1; d < dim; ++d) + { + upper_right[d] = 1; + repetitions.push_back(1); + } + + GridGenerator::subdivided_hyper_rectangle(triangulation, + repetitions, + lower_left, + upper_right); +} + + + +template +void +Test::setup_discrete_level_set() +{ + dof_handler.distribute_dofs(fe_collection); + level_set.reinit(dof_handler.n_dofs()); + + Point point_on_zero_contour; + point_on_zero_contour[dim - 1] = 0.5; + const Functions::SignedDistance::Plane analytical_levelset( + point_on_zero_contour, Point::unit_vector(dim - 1)); + + VectorTools::interpolate(dof_handler, analytical_levelset, level_set); +} + + + +template +void +Test::print_which_optionals_have_values_on_cell_0( + NonMatching::FEInterfaceValues &fe_values) +{ + const auto cell = dof_handler.begin_active(); + + for (unsigned int face_index : cell->face_indices()) + { + deallog << "face " << face_index << std::endl; + + if (cell->at_boundary(face_index)) + { + fe_values.reinit(cell, face_index); + } + else + { + const unsigned int invalid_subface = + dealii::numbers::invalid_unsigned_int; + + fe_values.reinit(cell, + face_index, + invalid_subface, + cell->neighbor(face_index), + cell->neighbor_of_neighbor(face_index), + invalid_subface); + } + + if (fe_values.get_inside_fe_values()) + deallog << "inside has value" << std::endl; + + if (fe_values.get_outside_fe_values()) + deallog << "outside has value" << std::endl; + } +} + + + +template +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + Test test; + test.run(); + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/non_matching/fe_interface_values.output b/tests/non_matching/fe_interface_values.output new file mode 100644 index 0000000000..52bca2f643 --- /dev/null +++ b/tests/non_matching/fe_interface_values.output @@ -0,0 +1,73 @@ + +DEAL::dim = 1 +DEAL::Simple constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::face 1 +DEAL::outside has value +DEAL::Advanced constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::face 1 +DEAL::outside has value +DEAL:: +DEAL::dim = 2 +DEAL::Simple constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::outside has value +DEAL::face 1 +DEAL::inside has value +DEAL::outside has value +DEAL::face 2 +DEAL::inside has value +DEAL::face 3 +DEAL::outside has value +DEAL::Advanced constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::outside has value +DEAL::face 1 +DEAL::inside has value +DEAL::outside has value +DEAL::face 2 +DEAL::inside has value +DEAL::face 3 +DEAL::outside has value +DEAL:: +DEAL::dim = 3 +DEAL::Simple constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::outside has value +DEAL::face 1 +DEAL::inside has value +DEAL::outside has value +DEAL::face 2 +DEAL::inside has value +DEAL::outside has value +DEAL::face 3 +DEAL::inside has value +DEAL::outside has value +DEAL::face 4 +DEAL::inside has value +DEAL::face 5 +DEAL::outside has value +DEAL::Advanced constructor: +DEAL::face 0 +DEAL::inside has value +DEAL::outside has value +DEAL::face 1 +DEAL::inside has value +DEAL::outside has value +DEAL::face 2 +DEAL::inside has value +DEAL::outside has value +DEAL::face 3 +DEAL::inside has value +DEAL::outside has value +DEAL::face 4 +DEAL::inside has value +DEAL::face 5 +DEAL::outside has value +DEAL:: diff --git a/tests/non_matching/fe_interface_values_non_standard_meshes.cc b/tests/non_matching/fe_interface_values_non_standard_meshes.cc new file mode 100644 index 0000000000..fcfa25cde8 --- /dev/null +++ b/tests/non_matching/fe_interface_values_non_standard_meshes.cc @@ -0,0 +1,206 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + + + +// For all combinations of 2-cell triangulations that can be created by +// GridGenerator::non_standard_orientation_mesh, set up a discrete level set +// function corresponding to +// (y - 0.25) + (z - 0.25) = 0 +// and use NonMatching::FEInterfaceValues to generate the face quadratures on +// the shared face at x = 0. Check that the quadrature points on both sides of +// the face are the same. + + +#include +#include +#include +#include + +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +class Test +{ +public: + Test(const Triangulation &triangulation); + + void + run(); + +private: + void + setup_discrete_level_set(); + + /* + * For the incoming FEInterfaceValues object, print the distance between the + * real space quadrature points of the underlying FEFaceValues objects, for + * all quadrature points. + */ + void + print_quadrature_point_difference_on_both_sides( + const dealii::FEInterfaceValues &fe_interface_values); + + hp::FECollection fe_collection; + DoFHandler dof_handler; + Vector level_set; + + NonMatching::MeshClassifier mesh_classifier; + + QGauss<1> quadrature_1D; +}; + + + +template +Test::Test(const Triangulation &triangulation) + : dof_handler(triangulation) + , mesh_classifier(dof_handler, level_set) + , quadrature_1D(2) +{ + fe_collection.push_back(FE_Q(1)); +} + + + +template +void +Test::setup_discrete_level_set() +{ + dof_handler.distribute_dofs(fe_collection); + level_set.reinit(dof_handler.n_dofs()); + + Point point_on_zero_contour; + Tensor<1, dim> plane_normal; + for (unsigned int i = 1; i < dim; i++) + { + point_on_zero_contour(i) = 0.25; + plane_normal[i] = 1; + } + + const Functions::SignedDistance::Plane analytical_levelset( + point_on_zero_contour, plane_normal); + + VectorTools::interpolate(dof_handler, analytical_levelset, level_set); +} + + + +template +void +Test::print_quadrature_point_difference_on_both_sides( + const dealii::FEInterfaceValues &fe_interface_values) +{ + for (unsigned int q : fe_interface_values.quadrature_point_indices()) + { + const Point &point_on_side_0 = + fe_interface_values.get_fe_face_values(0).quadrature_point(q); + const Point &point_on_side_1 = + fe_interface_values.get_fe_face_values(1).quadrature_point(q); + + deallog << point_on_side_0.distance(point_on_side_1) << std::endl; + } +} + + + +template +void +Test::run() +{ + setup_discrete_level_set(); + mesh_classifier.reclassify(); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = update_quadrature_points; + region_update_flags.outside = update_quadrature_points; + + NonMatching::FEInterfaceValues fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + dof_handler, + level_set); + + const auto cell = dof_handler.begin_active(); + + // Test that the points agree only on the interface between the two cells, + // which face this is depends on the incoming triangulation. + for (unsigned int face : cell->face_indices()) + if (!cell->face(face)->at_boundary()) + { + const unsigned int invalid_subface = + dealii::numbers::invalid_unsigned_int; + + fe_values.reinit(cell, + face, + invalid_subface, + cell->neighbor(face), + cell->neighbor_of_neighbor(face), + invalid_subface); + + Assert(fe_values.get_inside_fe_values(), ExcInternalError()); + Assert(fe_values.get_outside_fe_values(), ExcInternalError()); + + print_quadrature_point_difference_on_both_sides( + *fe_values.get_inside_fe_values()); + print_quadrature_point_difference_on_both_sides( + *fe_values.get_outside_fe_values()); + } +} + + + +int +main() +{ + initlog(); + const int dim = 3; + + // non_standard_orientation_mesh takes 4 boolean parameters. Test all + // combinations. + const unsigned int n_arguments = 4; + const unsigned int n_combinations = std::pow(2, 4); + for (unsigned int i = 0; i < n_combinations; i++) + { + const std::bitset bits(i); + + const bool face_orientation = bits[0]; + const bool face_flip = bits[1]; + const bool face_rotation = bits[2]; + const bool manipulate_left_cube = bits[3]; + + Triangulation triangulation; + GridGenerator::non_standard_orientation_mesh(triangulation, + face_orientation, + face_flip, + face_rotation, + manipulate_left_cube); + Test test(triangulation); + test.run(); + } +} diff --git a/tests/non_matching/fe_interface_values_non_standard_meshes.output b/tests/non_matching/fe_interface_values_non_standard_meshes.output new file mode 100644 index 0000000000..f8b01367a1 --- /dev/null +++ b/tests/non_matching/fe_interface_values_non_standard_meshes.output @@ -0,0 +1,193 @@ + +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::0 -- 2.39.5