From: Simon Sticko Date: Thu, 18 Nov 2021 21:42:03 +0000 (+0100) Subject: Add class NonMatching::FEValues X-Git-Tag: v9.4.0-rc1~780^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1a555f5cf879f7b53b541ab6d2e76ba5021b8b4e;p=dealii.git Add class NonMatching::FEValues This class makes it easier to assemble cut finite element methods, when the domain is described by a level set function. In the same way as hp::FEValues, NonMatching::FEValues generates dealii::FEValues like objects, but where these contain quadrature rules to integrate over the 3 different regions of the cell that are defined by the sign of the level set function. --- diff --git a/include/deal.II/non_matching/fe_values.h b/include/deal.II/non_matching/fe_values.h new file mode 100644 index 0000000000..1225cf66da --- /dev/null +++ b/include/deal.II/non_matching/fe_values.h @@ -0,0 +1,420 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_non_matching_fe_values +#define dealii_non_matching_fe_values + +#include +#include +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace NonMatching +{ + namespace internal + { + namespace FEValuesImplementation + { + template + class CellWiseFunction; + } + } // namespace internal + + + /** + * Struct storing UpdateFlags for the 3 regions of a cell, $K$, that is + * defined by the sign of a level set function, $\psi$: + * @f[ + * N = \{x \in K : \psi(x) < 0 \}, \\ + * P = \{x \in K : \psi(x) > 0 \}, \\ + * S = \{x \in K : \psi(x) = 0 \}. + * @f] + * As in the QuadratureGenerator class, we refer to $N$, $P$ and $S$ as the + * inside, outside, and surface region. RegionUpdateFlags is used to describe + * how the FEValues objects, which are created by NonMatching::FEValues, + * should be updated. + */ + struct RegionUpdateFlags + { + /** + * Constructor, sets the UpdateFlags for each region to update_default. + */ + RegionUpdateFlags(); + + /** + * Flags for the region $\{x \in K : \psi(x) < 0 \}$ + */ + UpdateFlags inside; + + /** + * Flags for the region $\{x \in K : \psi(x) > 0 \}$ + */ + UpdateFlags outside; + + /** + * Flags for the region $\{x \in K : \psi(x) = 0 \}$ + */ + UpdateFlags surface; + }; + + + /** + * This class is intended to facilitate assembling in immersed (in the sense + * of cut) finite element methods when the domain is described by a level set + * function, $\psi : \mathbb{R}^{dim} \to \mathbb{R}$. In this type of + * method, we typically need to integrate over 3 different regions of each + * cell, $K$: + * @f[ + * N = \{x \in K : \psi(x) < 0 \}, \\ + * P = \{x \in K : \psi(x) > 0 \}, \\ + * S = \{x \in K : \psi(x) = 0 \}. + * @f] + * Thus we need quadrature rules for these 3 regions: + * @image html immersed_quadratures.svg + * As in the QuadratureGenerator class, we refer to $N$, $P$, and $S$ as the + * inside, outside, and surface regions. The constructor of this class takes a + * discrete level set function described by a DoFHandler and a Vector. When + * the reinit() function is called, the QuadratureGenerator will be called in + * the background to create these immersed quadrature rules. This class then + * creates dealii::FEValues objects for the inside/outside regions and an + * FEImmersedSurfaceValues object for the surface region. These objects can + * then be accessed through one of the functions: get_inside_fe_values(), + * get_outside_fe_values(), or + * get_surface_fe_values(). + * Since a cut between a cell and the domain can be arbitrarily small, the + * underlying algorithm may generate a quadrature rule with 0 points. This + * can, for example, happen if the relative size of the cut is similar to the + * floating-point accuracy. Since the FEValues-like objects are not allowed to + * contain 0 points, the object that get_inside/outside/surface_fe_values() + * returns is wrapped in a std_cxx17::optional. This requires us to check if + * the returned FEValues-like object contains a value before we use it: + * @code + * NonMatching::FEValues fe_values(...); + * + * for (const auto &cell : dof_handler.active_cell_iterators()) + * { + * fe_values.reinit(cell); + * + * const std_cxx17::optional> &fe_values_inside = + * fe_values.get_inside_fe_values(); + * + * if (fe_values_inside) + * { + * // Assemble locally + * for (const unsigned int q_index : + * fe_values_inside->quadrature_point_indices()) + * { + * // ... + * } + * } + * } + * @endcode + * + * Of course, it is somewhat expensive to generate the immersed quadrature + * rules and create FEValues objects with the generated quadratures. 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. It only generates the immersed + * quadrature rules if the cell is intersected. If the cell is completely + * inside or outside, it returns a cached FEValues object created with a + * quadrature over the reference cell: $[0, 1]^{dim}$. + */ + template + class FEValues + { + public: + using AdditionalData = typename QuadratureGenerator::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 elements 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/surface 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 + FEValues(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}$ + * that should be used when a cell 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/surface 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 + FEValues(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 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_surface_fe_values(). + */ + template + void + reinit( + const TriaIterator> &cell); + + /** + * Return an dealii::FEValues 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::FEValues 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; + + /** + * Return an dealii::FEValues object reinitialized with a quadrature for the + * surface 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 not intersected, the returned optional will not contain a value. + */ + const std_cxx17::optional> & + get_surface_fe_values() const; + + private: + /** + * Do work common to the constructors. The incoming QCollection should be + * quadratures integrating over $[0, 1]^{dim}$. These will be used on the + * non-intersected cells. + */ + void + initialize(const hp::QCollection &q_collection); + + /** + * 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 cell that reinit was called with. + */ + LocationToLevelSet current_cell_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::FEValues object created with a quadrature + * rule over the full reference cell: $[0, 1]^{dim}$ 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::FEValues + * object in this container corresponding to the cell's active_fe_index. + * + * This container is a std::deque, which is compatible with the `FEValues` + * 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::FEValues object created with a quadrature + * rule over the full reference cell: $[0, 1]^{dim}$ 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 `FEValues` + * class that does not have a copy-constructor. + */ + std::deque>> + fe_values_outside_full_quadrature; + + /** + * FEValues 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; + + /** + * FEValues 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_outside; + + /** + * FEImmersedSurfaceValues object created with a quadrature rule integrating + * over the surface 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_surface; + + /** + * Object that generates the immersed quadrature rules. + */ + QuadratureGenerator quadrature_generator; + + /** + * Function that describes our level set function in reference space. + */ + const std::unique_ptr< + internal::FEValuesImplementation::CellWiseFunction> + reference_space_level_set; + }; + + + namespace internal + { + namespace FEValuesImplementation + { + /** + * Interface for a scalar Function which has a + * set_active_cell(..)-function. That is, a function which we in some way + * need to associate with a given cell in order to evaluate. + */ + template + class CellWiseFunction : public Function + { + public: + /** + * Destructor. Declared to make it virtual. + */ + virtual ~CellWiseFunction() = default; + + /** + * Set the cell that the function should be evaluated on. + */ + virtual void + set_active_cell( + const typename Triangulation::active_cell_iterator &cell) = 0; + }; + } // namespace FEValuesImplementation + } // namespace internal + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt index aa3d565b8b..f4b2d1d0f8 100644 --- a/source/non_matching/CMakeLists.txt +++ b/source/non_matching/CMakeLists.txt @@ -17,6 +17,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src fe_immersed_values.cc + fe_values.cc mesh_classifier.cc quadrature_generator.cc coupling.cc @@ -25,6 +26,7 @@ SET(_src SET(_inst fe_immersed_values.inst.in + fe_values.inst.in mesh_classifier.inst.in quadrature_generator.inst.in coupling.inst.in diff --git a/source/non_matching/fe_values.cc b/source/non_matching/fe_values.cc new file mode 100644 index 0000000000..2b26e6aa85 --- /dev/null +++ b/source/non_matching/fe_values.cc @@ -0,0 +1,523 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace NonMatching +{ + RegionUpdateFlags::RegionUpdateFlags() + : inside(update_default) + , outside(update_default) + , surface(update_default) + {} + + + namespace internal + { + namespace FEValuesImplementation + { + DeclExceptionMsg( + ExcCellNotSet, + "The set_active_cell function has to be called before calling this function."); + + + /** + * This class evaluates a function defined by a solution vector and a + * DoFHandler transformed to reference space. To be precise, if we let + * $\hat{x}$ be a point on the reference cell, this class implements the + * function + * + * $\hat{f}(\hat{x}) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(\hat{x})$, + * + * where $f_j$ are the local solution values and $\hat{\phi}_j(\hat(x))$ + * are the local reference space shape functions. The gradient and Hessian + * of this function are thus derivatives with respect to the reference + * space coordinates, $\hat{x}_0, \hat{x}_1, \ldots$. + * + * Note that this class is similar to FEFieldFunction, but that + * FEFieldFunction implements the following function on a given cell, $K$, + * + * $f(x) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(F_K^{-1}(x))$, + * + * which has the same coefficients but uses real space basis functions. + * Here, $F_K$ is the mapping from the reference cell to the real cell. + * + * Before calling the value/gradient/hessian function, the set_active_cell + * function must be called to specify which cell the function should be + * evaluated on. + */ + template > + class RefSpaceFEFieldFunction : public CellWiseFunction + { + public: + /** + * Constructor. Takes the solution vector and the associated DoFHandler. + * + * Pointers to the input arguments are stored internally, so they must + * have a longer lifetime than the created RefSpaceFEFieldFunction + * object. + */ + RefSpaceFEFieldFunction(const DoFHandler &dof_handler, + const VectorType & dof_values); + + /** + * @copydoc CellWiseFunction::set_active_cell() + */ + void + set_active_cell(const typename Triangulation::active_cell_iterator + &cell) override; + + /** + * @copydoc Function::value() + * + * @note The set_active_cell function must be called before this function. + * The incoming point should be on the reference cell, but this is not + * checked. + */ + double + value(const Point & point, + const unsigned int component = 0) const override; + + /** + * @copydoc Function::gradient() + * + * @note The set_active_cell function must be called before this function. + * The incoming point should be on the reference cell, but this is not + * checked. + */ + Tensor<1, dim> + gradient(const Point & point, + const unsigned int component = 0) const override; + + /** + * @copydoc Function::hessian() + * + * @note The set_active_cell function must be called before this function. + * The incoming point should be on the reference cell, but this is not + * checked. + */ + SymmetricTensor<2, dim> + hessian(const Point & point, + const unsigned int component = 0) const override; + + private: + /** + * Return whether the set_active_cell function has been called. + */ + bool + cell_is_set() const; + + /** + * Pointer to the DoFHandler passed to the constructor. + */ + const SmartPointer> dof_handler; + + /** + * Pointer to the vector of solution coefficients passed to the + * constructor. + */ + const SmartPointer global_dof_values; + + /** + * Pointer to the element associated with the cell in the last call to + * set_active_cell(). + */ + SmartPointer> element; + + /** + * DOF-indices of the cell in the last call to set_active_cell() + */ + std::vector local_dof_indices; + + /** + * Local solution values of the cell in the last call to + * set_active_cell() + */ + std::vector local_dof_values; + }; + + + + template + RefSpaceFEFieldFunction::RefSpaceFEFieldFunction( + const DoFHandler &dof_handler, + const VectorType & dof_values) + : dof_handler(&dof_handler) + , global_dof_values(&dof_values) + { + Assert(dof_handler.n_dofs() == dof_values.size(), + ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size())); + } + + + + template + void + RefSpaceFEFieldFunction::set_active_cell( + const typename Triangulation::active_cell_iterator &cell) + { + Assert( + &cell->get_triangulation() == &dof_handler->get_triangulation(), + ExcMessage( + "The incoming cell must belong to the triangulation associated with " + "the DoFHandler passed to the constructor.")); + + const typename DoFHandler::active_cell_iterator dof_handler_cell( + &dof_handler->get_triangulation(), + cell->level(), + cell->index(), + dof_handler); + + // Save the element and the local dof values, since this is what we need + // to evaluate the function. + element = &dof_handler_cell->get_fe(); + + local_dof_indices.resize(element->dofs_per_cell); + dof_handler_cell->get_dof_indices(local_dof_indices); + + local_dof_values.resize(element->dofs_per_cell); + + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + local_dof_values[i] = + dealii::internal::ElementAccess::get( + *global_dof_values, local_dof_indices[i]); + } + + + + template + bool + RefSpaceFEFieldFunction::cell_is_set() const + { + // If set cell hasn't been called the size of local_dof_values will be + // zero. + return local_dof_values.size() > 0; + } + + + + template + double + RefSpaceFEFieldFunction::value( + const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + Assert(cell_is_set(), ExcCellNotSet()); + + double value = 0; + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + value += local_dof_values[i] * + element->shape_value_component(i, point, component); + + return value; + } + + + + template + Tensor<1, dim> + RefSpaceFEFieldFunction::gradient( + const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + Assert(cell_is_set(), ExcCellNotSet()); + + Tensor<1, dim> gradient; + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + gradient += local_dof_values[i] * + element->shape_grad_component(i, point, component); + + return gradient; + } + + + + template + SymmetricTensor<2, dim> + RefSpaceFEFieldFunction::hessian( + const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + Assert(cell_is_set(), ExcCellNotSet()); + + Tensor<2, dim> hessian; + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + hessian += local_dof_values[i] * + element->shape_grad_grad_component(i, point, component); + + return symmetrize(hessian); + } + } // namespace FEValuesImplementation + } // namespace internal + + + + template + template + FEValues::FEValues(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) + , quadrature_generator(q_collection_1D, additional_data) + , reference_space_level_set( + std::make_unique>( + dof_handler, + level_set)) + { + // 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 + FEValues::FEValues(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) + , quadrature_generator(q_collection_1D, additional_data) + , reference_space_level_set( + std::make_unique>( + dof_handler, + level_set)) + { + initialize(q_collection); + } + + + + template + void + FEValues::initialize(const hp::QCollection &q_collection) + { + current_cell_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::FEValues 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 + FEValues::reinit( + const TriaIterator> &cell) + { + current_cell_location = mesh_classifier->location_to_level_set(cell); + 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_surface.reset(); + fe_values_outside.reset(); + + switch (current_cell_location) + { + case LocationToLevelSet::inside: + { + fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell); + break; + } + case LocationToLevelSet::outside: + { + fe_values_outside_full_quadrature.at(active_fe_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); + + reference_space_level_set->set_active_cell(cell); + const BoundingBox unit_box = create_unit_bounding_box(); + quadrature_generator.generate(*reference_space_level_set, unit_box); + + const Quadrature &inside_quadrature = + quadrature_generator.get_inside_quadrature(); + const Quadrature &outside_quadrature = + quadrature_generator.get_outside_quadrature(); + const ImmersedSurfaceQuadrature &surface_quadrature = + quadrature_generator.get_surface_quadrature(); + + // Even if a cell is formally intersected the number of created + // quadrature points can be 0. Avoid creating an FEValues 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); + + fe_values_inside->reinit(cell); + } + + if (outside_quadrature.size() > 0) + { + fe_values_outside.emplace((*mapping_collection)[mapping_index], + (*fe_collection)[active_fe_index], + outside_quadrature, + region_update_flags.outside); + + fe_values_outside->reinit(cell); + } + + if (surface_quadrature.size() > 0) + { + fe_values_surface.emplace((*mapping_collection)[mapping_index], + (*fe_collection)[active_fe_index], + surface_quadrature, + region_update_flags.surface); + fe_values_surface->reinit(cell); + } + + break; + } + default: + { + Assert(false, ExcInternalError()); + break; + } + } + } + + + + template + const std_cxx17::optional> & + FEValues::get_inside_fe_values() const + { + if (current_cell_location == LocationToLevelSet::inside) + return fe_values_inside_full_quadrature.at(active_fe_index); + else + return fe_values_inside; + } + + + + template + const std_cxx17::optional> & + FEValues::get_outside_fe_values() const + { + if (current_cell_location == LocationToLevelSet::outside) + return fe_values_outside_full_quadrature.at(active_fe_index); + else + return fe_values_outside; + } + + + + template + const std_cxx17::optional> & + FEValues::get_surface_fe_values() const + { + return fe_values_surface; + } + + +#include "fe_values.inst" + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/fe_values.inst.in b/source/non_matching/fe_values.inst.in new file mode 100644 index 0000000000..e94c54ec6f --- /dev/null +++ b/source/non_matching/fe_values.inst.in @@ -0,0 +1,55 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2014 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 (deal_II_dimension : DIMENSIONS) + { + template class FEValues; + } + +for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) + { + template class internal::FEValuesImplementation:: + RefSpaceFEFieldFunction; + + template FEValues::FEValues( + const hp::MappingCollection &, + const hp::FECollection &, + const hp::QCollection &, + const hp::QCollection<1> &, + const RegionUpdateFlags, + const MeshClassifier &, + const DoFHandler &, + const VEC &, + const typename FEValues::AdditionalData &); + + + template FEValues::FEValues( + const hp::FECollection &, + const Quadrature<1> &, + const RegionUpdateFlags, + const MeshClassifier &, + const DoFHandler &, + const VEC &, + const typename FEValues::AdditionalData &); + } + +// Template reinit function +for (deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template void FEValues::reinit( + const TriaIterator< + DoFCellAccessor> &); + } diff --git a/tests/non_matching/fe_values.cc b/tests/non_matching/fe_values.cc new file mode 100644 index 0000000000..0e15ab086f --- /dev/null +++ b/tests/non_matching/fe_values.cc @@ -0,0 +1,260 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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 + +#include "../tests.h" + +using namespace dealii; + +// Assert that the two incoming cells are the same. Throw an exception if not. +template +void +assert_cells_are_the_same( + const typename Triangulation::cell_iterator &expected, + const typename Triangulation::cell_iterator &cell) +{ + AssertThrow(expected == cell, ExcInternalError()); + deallog << "OK" << std::endl; +} + + + +// Assert that the incoming optional does not have a value. +template +void +assert_doesnt_have_value(const std_cxx17::optional &fe_values) +{ + AssertThrow(!fe_values, ExcInternalError()); + deallog << "OK" << std::endl; +} + + + +// Set up a triangulation with 3 elements in a row |-0-|-1-|-2-|. +// and a level set function such that +// +// element LocationToLevelSet +// 0 inside +// 1 intersected +// 2 outside +// +// Test the following: +// 1. That we can construct a NonMatching::FEValues object. +// 2. That when we call reinit on each cell in the triangulation, +// the boost::optionals that we get from get_inside/outside/surface_fe_values +// are set up correctly. That is, the 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(); + + // Setup a discrete level set function corresponding to + // $\psi(x) = (x_0 - 1.5) = 0$ + void + setup_discrete_level_set(); + + void + test_fe_values_reinitializes_correctly( + NonMatching::FEValues &fe_values) const; + + Triangulation triangulation; + hp::FECollection fe_collection; + hp::DoFHandler dof_handler; + + hp::MappingCollection mapping_collection; + hp::QCollection q_collection; + hp::QCollection<1> q_collection1D; + + 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(); + dof_handler.distribute_dofs(fe_collection); + setup_discrete_level_set(); + mesh_classifier.reclassify(); + + const NonMatching::RegionUpdateFlags region_update_flags; + + { + // Test with the "simple" constructor. + NonMatching::FEValues fe_values(fe_collection, + q_collection1D[0], + region_update_flags, + mesh_classifier, + dof_handler, + level_set); + test_fe_values_reinitializes_correctly(fe_values); + } + { + // Test with the "more advanced" constructor. + NonMatching::FEValues fe_values(mapping_collection, + fe_collection, + q_collection, + q_collection1D, + region_update_flags, + mesh_classifier, + dof_handler, + level_set); + test_fe_values_reinitializes_correctly(fe_values); + } +} + + + +template +void +Test::setup_mesh() +{ + const Point lower_left; + Point upper_right; + + std::vector repetitions; + upper_right[0] = 3; + repetitions.push_back(3); + 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() +{ + Point point_on_zero_contour; + point_on_zero_contour[0] = 1.5; + const Functions::LevelSet::Plane analytical_levelset( + point_on_zero_contour, Point::unit_vector(0)); + + level_set.reinit(dof_handler.n_dofs()); + VectorTools::interpolate(dof_handler, analytical_levelset, level_set); +} + + + +template +void +Test::test_fe_values_reinitializes_correctly( + NonMatching::FEValues &fe_values) 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, + fe_values.get_inside_fe_values()->get_cell()); + assert_doesnt_have_value(fe_values.get_surface_fe_values()); + assert_doesnt_have_value(fe_values.get_outside_fe_values()); + + + // The second is intersected so all FEValues object should be + // initialized. + cell++; + fe_values.reinit(cell); + + assert_cells_are_the_same(cell, + fe_values.get_inside_fe_values()->get_cell()); + assert_cells_are_the_same(cell, + fe_values.get_outside_fe_values()->get_cell()); + assert_cells_are_the_same(cell, + fe_values.get_surface_fe_values()->get_cell()); + + // The third is outside so only the outside FEValues object should be + // initialized. + cell++; + fe_values.reinit(cell); + + assert_doesnt_have_value(fe_values.get_inside_fe_values()); + assert_doesnt_have_value(fe_values.get_surface_fe_values()); + assert_cells_are_the_same(cell, + fe_values.get_outside_fe_values()->get_cell()); +} + + + +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_values.output b/tests/non_matching/fe_values.output new file mode 100644 index 0000000000..bd38077635 --- /dev/null +++ b/tests/non_matching/fe_values.output @@ -0,0 +1,61 @@ + +DEAL::dim = 1 +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 +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 +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::