--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/bounding_box.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/std_cxx17/optional.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/non_matching/fe_immersed_values.h>
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <deque>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace NonMatching
+{
+ namespace internal
+ {
+ namespace FEValuesImplementation
+ {
+ template <int dim>
+ 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<dim> fe_values(...);
+ *
+ * for (const auto &cell : dof_handler.active_cell_iterators())
+ * {
+ * fe_values.reinit(cell);
+ *
+ * const std_cxx17::optional<FEValues<dim>> &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 <int dim>
+ class FEValues
+ {
+ public:
+ using AdditionalData = typename QuadratureGenerator<dim>::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 <class VectorType>
+ FEValues(const hp::FECollection<dim> &fe_collection,
+ const Quadrature<1> & quadrature,
+ const RegionUpdateFlags region_update_flags,
+ const MeshClassifier<dim> & mesh_classifier,
+ const DoFHandler<dim> & 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 <class VectorType>
+ FEValues(const hp::MappingCollection<dim> &mapping_collection,
+ const hp::FECollection<dim> & fe_collection,
+ const hp::QCollection<dim> & q_collection,
+ const hp::QCollection<1> & q_collection_1D,
+ const RegionUpdateFlags region_update_flags,
+ const MeshClassifier<dim> & mesh_classifier,
+ const DoFHandler<dim> & 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 <bool level_dof_access>
+ void
+ reinit(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &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<dealii::FEValues<dim>> &
+ 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<dealii::FEValues<dim>> &
+ 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<FEImmersedSurfaceValues<dim>> &
+ 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<dim> &q_collection);
+
+ /**
+ * A pointer to the collection of mappings to be used.
+ */
+ const SmartPointer<const hp::MappingCollection<dim>> mapping_collection;
+
+ /**
+ * A pointer to the collection of finite elements to be used.
+ */
+ const SmartPointer<const hp::FECollection<dim>> 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<const MeshClassifier<dim>> 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<std_cxx17::optional<dealii::FEValues<dim>>>
+ 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<std_cxx17::optional<dealii::FEValues<dim>>>
+ 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<dealii::FEValues<dim>> 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<dealii::FEValues<dim>> 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<NonMatching::FEImmersedSurfaceValues<dim>>
+ fe_values_surface;
+
+ /**
+ * Object that generates the immersed quadrature rules.
+ */
+ QuadratureGenerator<dim> quadrature_generator;
+
+ /**
+ * Function that describes our level set function in reference space.
+ */
+ const std::unique_ptr<
+ internal::FEValuesImplementation::CellWiseFunction<dim>>
+ 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 <int dim>
+ class CellWiseFunction : public Function<dim>
+ {
+ 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<dim>::active_cell_iterator &cell) = 0;
+ };
+ } // namespace FEValuesImplementation
+ } // namespace internal
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
SET(_src
fe_immersed_values.cc
+ fe_values.cc
mesh_classifier.cc
quadrature_generator.cc
coupling.cc
SET(_inst
fe_immersed_values.inst.in
+ fe_values.inst.in
mesh_classifier.inst.in
quadrature_generator.inst.in
coupling.inst.in
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/exceptions.h>
+#include <deal.II/base/numbers.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/non_matching/fe_values.h>
+
+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 <int dim, class VectorType = Vector<double>>
+ class RefSpaceFEFieldFunction : public CellWiseFunction<dim>
+ {
+ 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<dim> &dof_handler,
+ const VectorType & dof_values);
+
+ /**
+ * @copydoc CellWiseFunction::set_active_cell()
+ */
+ void
+ set_active_cell(const typename Triangulation<dim>::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<dim> & 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<dim> & 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<dim> & 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<const DoFHandler<dim>> dof_handler;
+
+ /**
+ * Pointer to the vector of solution coefficients passed to the
+ * constructor.
+ */
+ const SmartPointer<const VectorType> global_dof_values;
+
+ /**
+ * Pointer to the element associated with the cell in the last call to
+ * set_active_cell().
+ */
+ SmartPointer<const FiniteElement<dim>> element;
+
+ /**
+ * DOF-indices of the cell in the last call to set_active_cell()
+ */
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ /**
+ * Local solution values of the cell in the last call to
+ * set_active_cell()
+ */
+ std::vector<typename VectorType::value_type> local_dof_values;
+ };
+
+
+
+ template <int dim, class VectorType>
+ RefSpaceFEFieldFunction<dim, VectorType>::RefSpaceFEFieldFunction(
+ const DoFHandler<dim> &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 <int dim, class VectorType>
+ void
+ RefSpaceFEFieldFunction<dim, VectorType>::set_active_cell(
+ const typename Triangulation<dim>::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<dim>::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<VectorType>::get(
+ *global_dof_values, local_dof_indices[i]);
+ }
+
+
+
+ template <int dim, class VectorType>
+ bool
+ RefSpaceFEFieldFunction<dim, VectorType>::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 <int dim, class VectorType>
+ double
+ RefSpaceFEFieldFunction<dim, VectorType>::value(
+ const Point<dim> & 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 <int dim, class VectorType>
+ Tensor<1, dim>
+ RefSpaceFEFieldFunction<dim, VectorType>::gradient(
+ const Point<dim> & 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 <int dim, class VectorType>
+ SymmetricTensor<2, dim>
+ RefSpaceFEFieldFunction<dim, VectorType>::hessian(
+ const Point<dim> & 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 <int dim>
+ template <class VectorType>
+ FEValues<dim>::FEValues(const hp::FECollection<dim> &fe_collection,
+ const Quadrature<1> & quadrature,
+ const RegionUpdateFlags region_update_flags,
+ const MeshClassifier<dim> & mesh_classifier,
+ const DoFHandler<dim> & dof_handler,
+ const VectorType & level_set,
+ const AdditionalData & additional_data)
+ : mapping_collection(&dealii::hp::StaticMappingQ1<dim>::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<internal::FEValuesImplementation::
+ RefSpaceFEFieldFunction<dim, VectorType>>(
+ dof_handler,
+ level_set))
+ {
+ // Tensor products of each quadrature in q_collection_1D. Used on the
+ // non-intersected cells.
+ hp::QCollection<dim> q_collection;
+ for (unsigned int i = 0; i < q_collection_1D.size(); ++i)
+ q_collection.push_back(Quadrature<dim>(q_collection_1D[i]));
+
+ initialize(q_collection);
+ }
+
+
+
+ template <int dim>
+ template <class VectorType>
+ FEValues<dim>::FEValues(const hp::MappingCollection<dim> &mapping_collection,
+ const hp::FECollection<dim> & fe_collection,
+ const hp::QCollection<dim> & q_collection,
+ const hp::QCollection<1> & q_collection_1D,
+ const RegionUpdateFlags region_update_flags,
+ const MeshClassifier<dim> & mesh_classifier,
+ const DoFHandler<dim> & 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<internal::FEValuesImplementation::
+ RefSpaceFEFieldFunction<dim, VectorType>>(
+ dof_handler,
+ level_set))
+ {
+ initialize(q_collection);
+ }
+
+
+
+ template <int dim>
+ void
+ FEValues<dim>::initialize(const hp::QCollection<dim> &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<dim> 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 <int dim>
+ template <bool level_dof_access>
+ void
+ FEValues<dim>::reinit(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &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<dim> unit_box = create_unit_bounding_box<dim>();
+ quadrature_generator.generate(*reference_space_level_set, unit_box);
+
+ const Quadrature<dim> &inside_quadrature =
+ quadrature_generator.get_inside_quadrature();
+ const Quadrature<dim> &outside_quadrature =
+ quadrature_generator.get_outside_quadrature();
+ const ImmersedSurfaceQuadrature<dim> &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 <int dim>
+ const std_cxx17::optional<dealii::FEValues<dim>> &
+ FEValues<dim>::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 <int dim>
+ const std_cxx17::optional<dealii::FEValues<dim>> &
+ FEValues<dim>::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 <int dim>
+ const std_cxx17::optional<FEImmersedSurfaceValues<dim>> &
+ FEValues<dim>::get_surface_fe_values() const
+ {
+ return fe_values_surface;
+ }
+
+
+#include "fe_values.inst"
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<deal_II_dimension>;
+ }
+
+for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
+ {
+ template class internal::FEValuesImplementation::
+ RefSpaceFEFieldFunction<deal_II_dimension, VEC>;
+
+ template FEValues<deal_II_dimension>::FEValues(
+ const hp::MappingCollection<deal_II_dimension> &,
+ const hp::FECollection<deal_II_dimension> &,
+ const hp::QCollection<deal_II_dimension> &,
+ const hp::QCollection<1> &,
+ const RegionUpdateFlags,
+ const MeshClassifier<deal_II_dimension> &,
+ const DoFHandler<deal_II_dimension> &,
+ const VEC &,
+ const typename FEValues<deal_II_dimension>::AdditionalData &);
+
+
+ template FEValues<deal_II_dimension>::FEValues(
+ const hp::FECollection<deal_II_dimension> &,
+ const Quadrature<1> &,
+ const RegionUpdateFlags,
+ const MeshClassifier<deal_II_dimension> &,
+ const DoFHandler<deal_II_dimension> &,
+ const VEC &,
+ const typename FEValues<deal_II_dimension>::AdditionalData &);
+ }
+
+// Template reinit function
+for (deal_II_dimension : DIMENSIONS; lda : BOOL)
+ {
+ template void FEValues<deal_II_dimension>::reinit(
+ const TriaIterator<
+ DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/exceptions.h>
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx17/optional.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/non_matching/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+// Assert that the two incoming cells are the same. Throw an exception if not.
+template <int dim>
+void
+assert_cells_are_the_same(
+ const typename Triangulation<dim>::cell_iterator &expected,
+ const typename Triangulation<dim>::cell_iterator &cell)
+{
+ AssertThrow(expected == cell, ExcInternalError());
+ deallog << "OK" << std::endl;
+}
+
+
+
+// Assert that the incoming optional does not have a value.
+template <class FEVALUES>
+void
+assert_doesnt_have_value(const std_cxx17::optional<FEVALUES> &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 <int dim>
+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<dim> &fe_values) const;
+
+ Triangulation<dim> triangulation;
+ hp::FECollection<dim> fe_collection;
+ hp::DoFHandler<dim> dof_handler;
+
+ hp::MappingCollection<dim> mapping_collection;
+ hp::QCollection<dim> q_collection;
+ hp::QCollection<1> q_collection1D;
+
+ Vector<double> level_set;
+ NonMatching::MeshClassifier<dim> mesh_classifier;
+};
+
+
+
+template <int dim>
+Test<dim>::Test()
+ : dof_handler(triangulation)
+ , mesh_classifier(dof_handler, level_set)
+{
+ fe_collection.push_back(FE_Q<dim>(1));
+ mapping_collection.push_back(MappingCartesian<dim>());
+ const unsigned int n_quadrature_points = 1;
+ q_collection.push_back(QGauss<dim>(n_quadrature_points));
+ q_collection1D.push_back(QGauss<1>(n_quadrature_points));
+}
+
+
+
+template <int dim>
+void
+Test<dim>::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<dim> 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<dim> 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 <int dim>
+void
+Test<dim>::setup_mesh()
+{
+ const Point<dim> lower_left;
+ Point<dim> upper_right;
+
+ std::vector<unsigned int> 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 <int dim>
+void
+Test<dim>::setup_discrete_level_set()
+{
+ Point<dim> point_on_zero_contour;
+ point_on_zero_contour[0] = 1.5;
+ const Functions::LevelSet::Plane<dim> analytical_levelset(
+ point_on_zero_contour, Point<dim>::unit_vector(0));
+
+ level_set.reinit(dof_handler.n_dofs());
+ VectorTools::interpolate(dof_handler, analytical_levelset, level_set);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_fe_values_reinitializes_correctly(
+ NonMatching::FEValues<dim> &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<dim>(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<dim>(cell,
+ fe_values.get_inside_fe_values()->get_cell());
+ assert_cells_are_the_same<dim>(cell,
+ fe_values.get_outside_fe_values()->get_cell());
+ assert_cells_are_the_same<dim>(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<dim>(cell,
+ fe_values.get_outside_fe_values()->get_cell());
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ deallog << "dim = " << dim << std::endl;
+ Test<dim> test;
+ test.run();
+ deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff