]> https://gitweb.dealii.org/ - dealii.git/commitdiff
QuadratureGenerator with discrete levelset 13615/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 12 Apr 2022 11:33:45 +0000 (13:33 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Thu, 21 Apr 2022 11:17:30 +0000 (13:17 +0200)
include/deal.II/non_matching/fe_values.h
include/deal.II/non_matching/quadrature_generator.h
source/non_matching/fe_values.cc
source/non_matching/fe_values.inst.in
source/non_matching/quadrature_generator.cc
source/non_matching/quadrature_generator.inst.in
tests/non_matching/discrete_quadrature_generator.cc [new file with mode: 0644]
tests/non_matching/discrete_quadrature_generator.output [new file with mode: 0644]

index 1225cf66dafb05fed2c4c5cc2b34ad368d53599c..dfa7b79b081631ee92ab522f594ef6e908e58515 100644 (file)
@@ -41,16 +41,6 @@ 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$:
@@ -375,45 +365,9 @@ namespace NonMatching
     /**
      * 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;
+    DiscreteQuadratureGenerator<dim> quadrature_generator;
   };
 
-
-  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
 
index 97653b097a26bf4c68fee58e2cafc854466e3718..38dfae22e0529df82ae7a070a16a1a74cd5b00ea 100644 (file)
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/std_cxx17/optional.h>
 
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/tria.h>
+
 #include <deal.II/hp/q_collection.h>
 
 #include <deal.II/non_matching/immersed_surface_quadrature.h>
@@ -39,8 +43,15 @@ namespace NonMatching
     {
       template <int dim, int spacedim>
       class QGenerator;
-    }
-  } // namespace internal
+    } // namespace QuadratureGeneratorImplementation
+
+
+    namespace DiscreteQuadratureGeneratorImplementation
+    {
+      template <int dim>
+      class CellWiseFunction;
+    } // namespace DiscreteQuadratureGeneratorImplementation
+  }   // namespace internal
 
 
   /**
@@ -459,6 +470,118 @@ namespace NonMatching
   };
 
 
+  /**
+   * This class generates the same type of immersed quadrature rules as those
+   * described in the QuadratureGenerator class. The difference is that this
+   * class handles the case when the the domain is a discrete level set
+   * function, i.e., when the level set function is described as a
+   * (DoFHandler, Vector) pair. The generate()-function of this class takes a
+   * cell in real space and constructs the immersed quadrature rules in
+   * reference space over this cell. These quadrature rules can then be obtained
+   * with one of the functions:
+   * get_inside_quadrature(),
+   * get_outside_quadrature(), and
+   * get_surface_quadrature().
+   *
+   * Internally, the quadrature generation is done by transforming the discrete
+   * level set function from real space to reference space and using the same
+   * algorithm as in the QuadratureGenerator class.
+   */
+  template <int dim>
+  class DiscreteQuadratureGenerator : public QuadratureGenerator<dim>
+  {
+  public:
+    using AdditionalData = AdditionalQGeneratorData;
+
+    /**
+     * Constructor, the discrete level set function is described by the
+     * incoming DoFHandler and Vector. Pointers to these are stored
+     * internally, so they must have a longer lifetime than the created this
+     * class. The hp::QCollection<1> and AdditionalData is passed to the
+     * QuadratureGenerator class.
+     */
+    template <class VectorType>
+    DiscreteQuadratureGenerator(
+      const hp::QCollection<1> &quadratures1D,
+      const DoFHandler<dim> &   dof_handler,
+      const VectorType &        level_set,
+      const AdditionalData &    additional_data = AdditionalData());
+
+    /**
+     * Construct immersed quadratures rules based on the discrete level
+     * set vector over the incoming cell.
+     */
+    template <bool level_dof_access>
+    void
+    generate(
+      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
+
+  private:
+    /**
+     * Function that describes our level set function in reference space.
+     */
+    std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
+                      CellWiseFunction<dim>>
+      reference_space_level_set;
+  };
+
+  /**
+   * This class generates the same type of immersed quadrature rules as those
+   * described in the FaceQuadratureGenerator class. The difference is that this
+   * class handles the case when the the domain is a discrete level set
+   * function, i.e., when the level set function is described as a
+   * (DoFHandler, Vector) pair. The generate()-function of this class takes a
+   * cell in real space plus the respective face index and constructs the
+   * immersed quadrature rules in reference space over this face. These
+   * quadrature rules can then be obtained with one of the functions:
+   * get_inside_quadrature(),
+   * get_outside_quadrature(), and
+   * get_surface_quadrature().
+   *
+   * Internally, the quadrature generation is done by transforming the discrete
+   * level set function from real space to reference space and using the same
+   * algorithm as in the FaceQuadratureGenerator class.
+   */
+  template <int dim>
+  class DiscreteFaceQuadratureGenerator : public FaceQuadratureGenerator<dim>
+  {
+  public:
+    using AdditionalData = AdditionalQGeneratorData;
+
+    /**
+     * Constructor, the discrete level set function is described by the
+     * incoming DoFHandler and Vector. Pointers to these are stored
+     * internally, so they must have a longer lifetime than the created this
+     * class. The hp::QCollection<1> and AdditionalData is passed to the
+     * QuadratureGenerator class.
+     */
+    template <class VectorType>
+    DiscreteFaceQuadratureGenerator(
+      const hp::QCollection<1> &quadratures1D,
+      const DoFHandler<dim> &   dof_handler,
+      const VectorType &        level_set,
+      const AdditionalData &    additional_data = AdditionalData());
+
+    /**
+     * Construct immersed quadratures rules based on the discrete level
+     * set vector over the incoming face described by cell and face index.
+     */
+    template <bool level_dof_access>
+    void
+    generate(
+      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+      const unsigned int face_index);
+
+  private:
+    /**
+     * Function that describes our level set function in reference space.
+     */
+    std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
+                      CellWiseFunction<dim>>
+      reference_space_level_set;
+  };
+
+
   namespace internal
   {
     namespace QuadratureGeneratorImplementation
@@ -1247,6 +1370,33 @@ namespace NonMatching
         const std::vector<FunctionBounds<dim>> &all_function_bounds);
 
     } // namespace QuadratureGeneratorImplementation
+
+
+    namespace DiscreteQuadratureGeneratorImplementation
+    {
+      /**
+       * 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 DiscreteQuadratureGeneratorImplementation
   }   // namespace internal
 
 } // namespace NonMatching
index d4e0bcd40a4fabd01dc828cd6fc4cc8d4c64c90d..75d82d252669b4fcf694badb9bee145cc49445f5 100644 (file)
@@ -43,332 +43,6 @@ namespace NonMatching
   {}
 
 
-  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;
-
-        /**
-         * Description of the 1D polynomial basis for tensor product elements
-         * used for the fast path of this class using tensor product
-         * evaluators.
-         */
-        std::vector<Polynomials::Polynomial<double>> poly;
-
-        /**
-         * Renumbering for the tensor-product evaluator in the fast path.
-         */
-        std::vector<unsigned int> renumber;
-
-        /**
-         * Check whether the shape functions are linear.
-         */
-        bool polynomials_are_hat_functions;
-      };
-
-
-
-      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.
-
-        // Check if we can use the fast path. In case we have a different
-        // element from the one used before we need to set up the data
-        // structures again.
-        if (element != &dof_handler_cell->get_fe())
-          {
-            poly.clear();
-            element = &dof_handler_cell->get_fe();
-
-            if (element->n_base_elements() == 1 &&
-                dealii::internal::FEPointEvaluation::is_fast_path_supported(
-                  *element, 0))
-              {
-                dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
-                  shape_info;
-
-                shape_info.reinit(QMidpoint<1>(), *element, 0);
-                renumber = shape_info.lexicographic_numbering;
-                poly =
-                  dealii::internal::FEPointEvaluation::get_polynomial_space(
-                    element->base_element(0));
-
-                polynomials_are_hat_functions =
-                  (poly.size() == 2 && poly[0].value(0.) == 1. &&
-                   poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
-                   poly[1].value(1.) == 1.);
-              }
-          }
-        else
-          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());
-
-        if (!poly.empty() && component == 0)
-          {
-            // TODO: this could be extended to a component that is not zero
-            return dealii::internal::evaluate_tensor_product_value_and_gradient(
-                     poly,
-                     local_dof_values,
-                     point,
-                     polynomials_are_hat_functions,
-                     renumber)
-              .first;
-          }
-        else
-          {
-            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());
-
-        if (!poly.empty() && component == 0)
-          {
-            // TODO: this could be extended to a component that is not zero
-            return dealii::internal::evaluate_tensor_product_value_and_gradient(
-                     poly,
-                     local_dof_values,
-                     point,
-                     polynomials_are_hat_functions,
-                     renumber)
-              .second;
-          }
-        else
-          {
-            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());
-
-        if (!poly.empty() && component == 0)
-          {
-            // TODO: this could be extended to a component that is not zero
-            return dealii::internal::evaluate_tensor_product_hessian(
-              poly, local_dof_values, point, renumber);
-          }
-        else
-          {
-            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>
@@ -384,12 +58,10 @@ namespace NonMatching
     , 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))
+    , 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.
@@ -418,12 +90,10 @@ namespace NonMatching
     , 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))
+    , quadrature_generator(q_collection_1D,
+                           dof_handler,
+                           level_set,
+                           additional_data)
   {
     initialize(q_collection);
   }
@@ -513,10 +183,7 @@ namespace NonMatching
             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);
+            quadrature_generator.generate(cell);
 
             const Quadrature<dim> &inside_quadrature =
               quadrature_generator.get_inside_quadrature();
index e94c54ec6f0d92b120f37cb039c52324bff3652e..e9cdf571c47bc521ecf8cb79e84e64e30dd6d395 100644 (file)
@@ -21,9 +21,6 @@ for (deal_II_dimension : DIMENSIONS)
 
 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> &,
index da551e41a1212648eff450fcec5e86b7673de4fe..1f5ba03a164f647ec675f1aae2c9e87495cfc774 100644 (file)
 
 #include <deal.II/grid/reference_cell.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/matrix_free/fe_point_evaluation.h>
+
 #include <deal.II/non_matching/quadrature_generator.h>
 
 #include <boost/math/special_functions/relative_difference.hpp>
@@ -1251,7 +1265,334 @@ namespace NonMatching
         AssertIndexRange(q_index, this->q_collection1D->size());
         this->q_index = q_index;
       }
+
+
+
     } // namespace QuadratureGeneratorImplementation
+
+
+
+    namespace DiscreteQuadratureGeneratorImplementation
+    {
+      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;
+
+        /**
+         * Description of the 1D polynomial basis for tensor product elements
+         * used for the fast path of this class using tensor product
+         * evaluators.
+         */
+        std::vector<Polynomials::Polynomial<double>> poly;
+
+        /**
+         * Renumbering for the tensor-product evaluator in the fast path.
+         */
+        std::vector<unsigned int> renumber;
+
+        /**
+         * Check whether the shape functions are linear.
+         */
+        bool polynomials_are_hat_functions;
+      };
+
+
+
+      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.
+
+        // Check if we can use the fast path. In case we have a different
+        // element from the one used before we need to set up the data
+        // structures again.
+        if (element != &dof_handler_cell->get_fe())
+          {
+            poly.clear();
+            element = &dof_handler_cell->get_fe();
+
+            if (element->n_base_elements() == 1 &&
+                dealii::internal::FEPointEvaluation::is_fast_path_supported(
+                  *element, 0))
+              {
+                dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
+                  shape_info;
+
+                shape_info.reinit(QMidpoint<1>(), *element, 0);
+                renumber = shape_info.lexicographic_numbering;
+                poly =
+                  dealii::internal::FEPointEvaluation::get_polynomial_space(
+                    element->base_element(0));
+
+                polynomials_are_hat_functions =
+                  (poly.size() == 2 && poly[0].value(0.) == 1. &&
+                   poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
+                   poly[1].value(1.) == 1.);
+              }
+          }
+        else
+          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());
+
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_value_and_gradient(
+                     poly,
+                     local_dof_values,
+                     point,
+                     polynomials_are_hat_functions,
+                     renumber)
+              .first;
+          }
+        else
+          {
+            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());
+
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_value_and_gradient(
+                     poly,
+                     local_dof_values,
+                     point,
+                     polynomials_are_hat_functions,
+                     renumber)
+              .second;
+          }
+        else
+          {
+            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());
+
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_hessian(
+              poly, local_dof_values, point, renumber);
+          }
+        else
+          {
+            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 DiscreteQuadratureGeneratorImplementation
   }   // namespace internal
 
 
@@ -1525,6 +1866,69 @@ namespace NonMatching
   {
     return surface_quadrature;
   }
+
+
+
+  template <int dim>
+  template <class VectorType>
+  DiscreteQuadratureGenerator<dim>::DiscreteQuadratureGenerator(
+    const hp::QCollection<1> &quadratures1D,
+    const DoFHandler<dim> &   dof_handler,
+    const VectorType &        level_set,
+    const AdditionalData &    additional_data)
+    : QuadratureGenerator<dim>(quadratures1D, additional_data)
+    , reference_space_level_set(
+        std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
+                           RefSpaceFEFieldFunction<dim, VectorType>>(
+          dof_handler,
+          level_set))
+  {}
+
+
+
+  template <int dim>
+  template <bool level_dof_access>
+  void
+  DiscreteQuadratureGenerator<dim>::generate(
+    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
+  {
+    reference_space_level_set->set_active_cell(cell);
+    const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
+    QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
+  }
+
+
+
+  template <int dim>
+  template <class VectorType>
+  DiscreteFaceQuadratureGenerator<dim>::DiscreteFaceQuadratureGenerator(
+    const hp::QCollection<1> &quadratures1D,
+    const DoFHandler<dim> &   dof_handler,
+    const VectorType &        level_set,
+    const AdditionalData &    additional_data)
+    : FaceQuadratureGenerator<dim>(quadratures1D, additional_data)
+    , reference_space_level_set(
+        std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
+                           RefSpaceFEFieldFunction<dim, VectorType>>(
+          dof_handler,
+          level_set))
+  {}
+
+
+
+  template <int dim>
+  template <bool level_dof_access>
+  void
+  DiscreteFaceQuadratureGenerator<dim>::generate(
+    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+    const unsigned int                                               face_index)
+  {
+    reference_space_level_set->set_active_cell(cell);
+    const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
+    FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
+                                           unit_box,
+                                           face_index);
+  }
 } // namespace NonMatching
 #include "quadrature_generator.inst"
 DEAL_II_NAMESPACE_CLOSE
index 3b6d03c08aaec5ac439bdc95dde487f1334bf4c3..d956a98bd16893dc2d2f01150ecd96ee03020720 100644 (file)
@@ -19,10 +19,12 @@ for (deal_II_dimension : DIMENSIONS)
     namespace NonMatching
     \{
       template class QuadratureGenerator<deal_II_dimension>;
+      template class DiscreteQuadratureGenerator<deal_II_dimension>;
 
 #if 1 < deal_II_dimension
       template class FaceQuadratureGenerator<deal_II_dimension>;
 #endif
+      template class DiscreteFaceQuadratureGenerator<deal_II_dimension>;
 
       namespace internal
       \{
@@ -76,3 +78,33 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
       UpThroughDimensionCreator<deal_II_dimension, deal_II_space_dimension>;
 #endif
   }
+
+for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
+  {
+    template NonMatching::DiscreteQuadratureGenerator<deal_II_dimension>::
+      DiscreteQuadratureGenerator(const hp::QCollection<1> &,
+                                  const DoFHandler<deal_II_dimension> &,
+                                  const VEC &,
+                                  const AdditionalData &);
+
+    template NonMatching::DiscreteFaceQuadratureGenerator<deal_II_dimension>::
+      DiscreteFaceQuadratureGenerator(const hp::QCollection<1> &,
+                                      const DoFHandler<deal_II_dimension> &,
+                                      const VEC &,
+                                      const AdditionalData &);
+  }
+
+// Template generate function
+for (deal_II_dimension : DIMENSIONS; lda : BOOL)
+  {
+    template void
+    NonMatching::DiscreteQuadratureGenerator<deal_II_dimension>::generate(
+      const TriaIterator<
+        DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+
+    template void
+    NonMatching::DiscreteFaceQuadratureGenerator<deal_II_dimension>::generate(
+      const TriaIterator<
+        DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
+      const unsigned int);
+  }
diff --git a/tests/non_matching/discrete_quadrature_generator.cc b/tests/non_matching/discrete_quadrature_generator.cc
new file mode 100644 (file)
index 0000000..8501821
--- /dev/null
@@ -0,0 +1,191 @@
+// ---------------------------------------------------------------------
+//
+// 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/fe/fe_q.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "quadrature_printing.h"
+
+using namespace dealii;
+
+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_discrete_quadrature_generator();
+
+  void
+  test_discrete_face_quadrature_generator();
+
+  Triangulation<dim>    triangulation;
+  hp::FECollection<dim> fe_collection;
+  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();
+  test_discrete_quadrature_generator();
+  test_discrete_face_quadrature_generator();
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_mesh()
+{
+  const Point<dim> lower_left;
+  Point<dim>       upper_right;
+  upper_right[0] = 2.;
+
+  for (unsigned int d = 1; d < dim; ++d)
+    {
+      upper_right[d] = 1.;
+    }
+
+  GridGenerator::hyper_rectangle(triangulation, 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_discrete_quadrature_generator()
+{
+  NonMatching::DiscreteQuadratureGenerator<dim> quadrature_generator(
+    q_collection1D, dof_handler, level_set);
+  quadrature_generator.generate(dof_handler.begin_active());
+
+  print_quadrature(quadrature_generator.get_inside_quadrature());
+  print_quadrature(quadrature_generator.get_outside_quadrature());
+  print_surface_quadrature(quadrature_generator.get_surface_quadrature());
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_discrete_face_quadrature_generator()
+{
+  NonMatching::DiscreteFaceQuadratureGenerator<dim> face_quadrature_generator(
+    q_collection1D, dof_handler, level_set);
+
+  const auto &cell = dof_handler.begin_active();
+  for (const auto f : cell->face_indices())
+    {
+      face_quadrature_generator.generate(cell, f);
+
+      print_quadrature(face_quadrature_generator.get_inside_quadrature());
+      print_quadrature(face_quadrature_generator.get_outside_quadrature());
+      print_surface_quadrature(
+        face_quadrature_generator.get_surface_quadrature());
+    }
+}
+
+
+
+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>();
+}
diff --git a/tests/non_matching/discrete_quadrature_generator.output b/tests/non_matching/discrete_quadrature_generator.output
new file mode 100644 (file)
index 0000000..4dfaa19
--- /dev/null
@@ -0,0 +1,40 @@
+
+DEAL::dim = 1
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::
+DEAL::dim = 2
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000
+DEAL::0.500000, 1.00000
+DEAL::0.500000, 1.00000
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::
+DEAL::dim = 3
+DEAL::0.375000, 0.500000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.500000, 1.00000
+DEAL::0.500000, 0.500000, 1.00000
+DEAL::0.500000, 0.375000, 0.750000
+DEAL::0.500000, 0.875000, 0.250000
+DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.375000, 0.750000
+DEAL::0.500000, 0.875000, 0.250000
+DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.