]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement paths for FE_Q_iso_Q1 in NonMatching 15803/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 18 Jul 2023 08:04:54 +0000 (10:04 +0200)
committerMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Wed, 9 Aug 2023 07:14:02 +0000 (09:14 +0200)
include/deal.II/non_matching/immersed_surface_quadrature.h
include/deal.II/non_matching/quadrature_generator.h
source/non_matching/immersed_surface_quadrature.cc
source/non_matching/mesh_classifier.cc
source/non_matching/quadrature_generator.cc
tests/non_matching/discrete_quadrature_generator.cc
tests/non_matching/discrete_quadrature_generator.output

index 764a7ede6f1a16c6fd9af23d6974d07a3670913d..9c315567537b1373a326dba17ef27eba21113ca7 100644 (file)
@@ -122,6 +122,12 @@ namespace NonMatching
                               const std::vector<double> &             weights,
                               const std::vector<Tensor<1, spacedim>> &normals);
 
+    /**
+     * Clears weights, points and normals vectors.
+     */
+    void
+    clear();
+
     /**
      * Extend the given formula by an additional quadrature point.
      * The point, weight and normal should be with respect to reference space,
index e9f5e7b00d6c107159afa8cbf77c5977fd05f45d..4612d95f5fbaf3ef5da732f02948e4af0c13ebe5 100644 (file)
@@ -201,6 +201,12 @@ namespace NonMatching
       const hp::QCollection<1> &quadratures1D,
       const AdditionalData &    additional_data = AdditionalData());
 
+    /**
+     * Clears the inside, outside and surface quadratures.
+     */
+    void
+    clear_quadratures();
+
     /**
      * Construct immersed quadratures rules for the incoming level set
      * function over the BoundingBox.
@@ -216,6 +222,14 @@ namespace NonMatching
     void
     generate(const Function<dim> &level_set, const BoundingBox<dim> &box);
 
+    /**
+     * Same as above but does not clear quadratures and appends it to the
+     * existing quadrature instead.
+     */
+    void
+    generate_append(const Function<dim> &   level_set,
+                    const BoundingBox<dim> &box);
+
     /**
      * Return the quadrature rule for the region
      * $\{x \in B : \psi(x) < 0 \}$
@@ -305,6 +319,12 @@ namespace NonMatching
       const hp::QCollection<1> &quadratures1D,
       const AdditionalData &    additional_data = AdditionalData());
 
+    /**
+     * Clears the inside, outside and surface quadratures.
+     */
+    void
+    clear_quadratures();
+
     /**
      * Construct immersed quadratures rules for the incoming level set
      * function on a given face of the BoundingBox.
@@ -322,6 +342,15 @@ namespace NonMatching
              const BoundingBox<dim> &box,
              const unsigned int      face_index);
 
+    /**
+     * Same as above but does not clear quadratures and appends it to the
+     * existing quadrature instead.
+     */
+    void
+    generate_append(const Function<dim> &   level_set,
+                    const BoundingBox<dim> &box,
+                    const unsigned int      face_index);
+
     /**
      * Return the quadrature rule for the region
      * $\{x \in F : \psi(x) < 0 \}$
@@ -403,6 +432,12 @@ namespace NonMatching
       const hp::QCollection<1> &quadratures1D,
       const AdditionalData &    additional_data = AdditionalData());
 
+    /**
+     * Does nothing. Exists for compatibility reasons.
+     */
+    void
+    clear_quadratures();
+
     /**
      * Construct immersed quadratures rules for the incoming level set
      * function on a given face of the BoundingBox.
@@ -417,6 +452,15 @@ namespace NonMatching
              const BoundingBox<1> &box,
              const unsigned int    face_index);
 
+    /**
+     * Same as above but does not clear quadratures and appends it to the
+     * existing quadrature instead.
+     */
+    void
+    generate_append(const Function<1> &   level_set,
+                    const BoundingBox<1> &box,
+                    const unsigned int    face_index);
+
     /**
      * @copydoc FaceQuadratureGenerator<dim>::get_inside_quadrature()
      */
@@ -518,6 +562,12 @@ namespace NonMatching
     generate(const typename Triangulation<dim>::active_cell_iterator &cell);
 
   private:
+    /**
+     * Construct immersed quadratures for FE_Q_iso_Q1.
+     */
+    void
+    generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box);
+
     /**
      * Function that describes our level set function in reference space.
      */
@@ -575,6 +625,13 @@ namespace NonMatching
              const unsigned int face_index);
 
   private:
+    /**
+     * Construct immersed quadratures for FE_Q_iso_Q1.
+     */
+    void
+    generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box,
+                         unsigned int            face_index);
+
     /**
      * Function that describes our level set function in reference space.
      */
@@ -701,6 +758,12 @@ namespace NonMatching
          */
         ExtendableQuadrature(const Quadrature<dim> &quadrature);
 
+        /**
+         * Clears weights and points vectors.
+         */
+        void
+        clear();
+
         /**
          * Add a point with an associated weight to the quadrature.
          */
@@ -760,6 +823,12 @@ namespace NonMatching
         ExtendableQuadrature<dim> &
         quadrature_by_definiteness(const Definiteness definiteness);
 
+        /**
+         * Clears all quadratures.
+         */
+        void
+        clear();
+
         /**
          * Quadrature for the region $\{x \in B : \psi_i(x) < 0 \forall i \}$ of
          * the box, $B$.
@@ -1396,6 +1465,26 @@ namespace NonMatching
         virtual void
         set_active_cell(
           const typename Triangulation<dim>::active_cell_iterator &cell) = 0;
+
+        /**
+         * Set the dof values and the bounding box of the subcell the
+         * function should be evaluated on. Relevant for FE_Q_iso_Q1.
+         */
+        virtual void
+        set_subcell(const std::vector<unsigned int> &mask,
+                    const BoundingBox<dim> &         subcell_box) = 0;
+
+        /**
+         * Returns flag indicating if the finite element is FE_Q_iso_Q1.
+         */
+        virtual bool
+        is_fe_q_iso_q1() const = 0;
+
+        /**
+         * Number of subdivisions of the FE_Q_iso_Q1 element.
+         */
+        virtual unsigned int
+        n_subdivisions() const = 0;
       };
 
     } // namespace DiscreteQuadratureGeneratorImplementation
index 55b6deef1506a84458c7476dac7a3ce1ae76a7bc..7e71e2b83437f5fe5bdf05f98b1fc6c11cd4930b 100644 (file)
@@ -38,6 +38,17 @@ namespace NonMatching
 
 
 
+  template <int dim, int spacedim>
+  inline void
+  ImmersedSurfaceQuadrature<dim, spacedim>::clear()
+  {
+    this->quadrature_points.clear();
+    this->weights.clear();
+    this->normals.clear();
+  }
+
+
+
   template <int dim, int spacedim>
   void
   ImmersedSurfaceQuadrature<dim, spacedim>::push_back(
index 743342bad1ccbfec877bd289580d066b3301a025..6134b58b9c935ad231bbb0007532021d0e29dbd1 100644 (file)
@@ -17,6 +17,7 @@
 
 #include <deal.II/dofs/dof_accessor.h>
 
+#include "deal.II/fe/fe_q_iso_q1.h"
 #include <deal.II/fe/fe_bernstein.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_values.h>
@@ -410,6 +411,24 @@ namespace NonMatching
                                                       face_index,
                                                       local_levelset_values);
 
+    const FiniteElement<dim> &fe =
+      level_set_description->get_fe_collection()[fe_index];
+
+    const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+      dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
+
+    const FE_Poly<dim> *fe_poly = dynamic_cast<const FE_Poly<dim> *>(&fe);
+
+    const bool is_linear = fe_q_iso_q1 != nullptr ||
+                           (fe_poly != nullptr && fe_poly->get_degree() == 1);
+
+    // shortcut for linear elements
+    if (is_linear)
+      {
+        return internal::MeshClassifierImplementation::location_from_dof_signs(
+          local_levelset_values);
+      }
+
     lagrange_to_bernstein_face[fe_index][face_index].solve(
       local_levelset_values);
 
@@ -466,7 +485,8 @@ namespace NonMatching
     for (unsigned int i = 0; i < fe_collection.size(); i++)
       {
         const FiniteElement<dim> &element = fe_collection[i];
-        const FE_Q<dim> *fe_q = dynamic_cast<const FE_Q<dim> *>(&element);
+        const FE_Q_Base<dim> *    fe_q =
+          dynamic_cast<const FE_Q_Base<dim> *>(&element);
         Assert(fe_q != nullptr, ExcNotImplemented());
 
         const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
index df5e9bb56cc8212fd46180e3cfeb63615bad5005..6e908042fe58220718cf482b70d2a7f934a642c7 100644 (file)
@@ -15,6 +15,9 @@
 
 #include <deal.II/base/function_tools.h>
 
+#include "deal.II/fe/fe_q.h"
+#include "deal.II/fe/fe_q_iso_q1.h"
+
 #include <deal.II/grid/reference_cell.h>
 
 #include <deal.II/lac/block_vector.h>
@@ -641,6 +644,16 @@ namespace NonMatching
 
 
 
+      template <int dim>
+      inline void
+      ExtendableQuadrature<dim>::clear()
+      {
+        this->quadrature_points.clear();
+        this->weights.clear();
+      }
+
+
+
       template <int dim>
       void
       ExtendableQuadrature<dim>::push_back(const Point<dim> &point,
@@ -670,6 +683,18 @@ namespace NonMatching
 
 
 
+      template <int dim>
+      void
+      QPartitioning<dim>::clear()
+      {
+        positive.clear();
+        negative.clear();
+        indefinite.clear();
+        surface.clear();
+      }
+
+
+
       /**
        * Takes a (dim-1)-dimensional point from the cross-section (orthogonal
        * to direction) of the box. Creates the two dim-dimensional points, which
@@ -870,7 +895,7 @@ namespace NonMatching
       void
       QGeneratorBase<dim, spacedim>::clear_quadratures()
       {
-        q_partitioning = QPartitioning<dim>();
+        q_partitioning.clear();
       }
 
 
@@ -1328,6 +1353,25 @@ namespace NonMatching
         set_active_cell(const typename Triangulation<dim>::active_cell_iterator
                           &cell) override;
 
+        /**
+         * @copydoc CellWiseFunction::set_subcell()
+         */
+        void
+        set_subcell(const std::vector<unsigned int> &mask,
+                    const BoundingBox<dim> &         subcell_box) override;
+
+        /**
+         * @copydoc CellWiseFunction::is_fe_q_iso_q1()
+         */
+        bool
+        is_fe_q_iso_q1() const override;
+
+        /**
+         * @copydoc CellWiseFunction::n_subdivisions()
+         */
+        unsigned int
+        n_subdivisions() const override;
+
         /**
          * @copydoc Function::value()
          *
@@ -1386,16 +1430,33 @@ namespace NonMatching
         SmartPointer<const FiniteElement<dim>> element;
 
         /**
-         * DOF-indices of the cell in the last call to set_active_cell()
+         * 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()
+         * set_active_cell().
          */
         std::vector<typename VectorType::value_type> local_dof_values;
 
+        /**
+         * Local solution values of the subcell after the last call to
+         * set_subcell().
+         */
+        std::vector<typename VectorType::value_type> local_dof_values_subcell;
+
+        /**
+         * Bounding box of the subcell after the last call to set_subcell().
+         */
+        BoundingBox<dim> subcell_box;
+
+        /**
+         * Number of subdivisions per line of the FE_Q_iso_Q1 element. Set to
+         * numbers::invalid_unsigned_int for other elements.
+         */
+        unsigned int n_subdivisions_per_line;
+
         /**
          * Description of the 1d polynomial basis for tensor product elements
          * used for the fast path of this class using tensor product
@@ -1422,6 +1483,7 @@ namespace NonMatching
         const VectorType &     dof_values)
         : dof_handler(&dof_handler)
         , global_dof_values(&dof_values)
+        , n_subdivisions_per_line(numbers::invalid_unsigned_int)
       {
         Assert(dof_handler.n_dofs() == dof_values.size(),
                ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
@@ -1443,6 +1505,8 @@ namespace NonMatching
         const auto dof_handler_cell =
           cell->as_dof_handler_iterator(*dof_handler);
 
+        const FE_Q<dim> fe_q1(1);
+
         // Save the element and the local dof values, since this is what we need
         // to evaluate the function.
 
@@ -1454,9 +1518,22 @@ namespace NonMatching
             poly.clear();
             element = &dof_handler_cell->get_fe();
 
-            if (element->n_base_elements() == 1 &&
-                dealii::internal::FEPointEvaluation::is_fast_path_supported(
-                  *element, 0))
+            const FiniteElement<dim> *fe = &dof_handler_cell->get_fe();
+
+            if (const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+                  dynamic_cast<const FE_Q_iso_Q1<dim> *>(
+                    &dof_handler_cell->get_fe()))
+              {
+                this->n_subdivisions_per_line = fe_q_iso_q1->get_degree();
+                fe                            = &fe_q1;
+                local_dof_values_subcell.resize(fe_q1.n_dofs_per_cell());
+              }
+            else
+              this->n_subdivisions_per_line = numbers::invalid_unsigned_int;
+
+            if (fe->n_base_elements() == 1 &&
+                dealii::internal::FEPointEvaluation::is_fast_path_supported(*fe,
+                                                                            0))
               {
                 dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
                   shape_info;
@@ -1465,7 +1542,7 @@ namespace NonMatching
                 renumber = shape_info.lexicographic_numbering;
                 poly =
                   dealii::internal::FEPointEvaluation::get_polynomial_space(
-                    element->base_element(0));
+                    fe->base_element(0));
 
                 polynomials_are_hat_functions =
                   (poly.size() == 2 && poly[0].value(0.) == 1. &&
@@ -1489,6 +1566,38 @@ namespace NonMatching
 
 
 
+      template <int dim, typename VectorType>
+      void
+      RefSpaceFEFieldFunction<dim, VectorType>::set_subcell(
+        const std::vector<unsigned int> &mask,
+        const BoundingBox<dim> &         subcell_box)
+      {
+        for (unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
+          local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]];
+
+        this->subcell_box = subcell_box;
+      }
+
+
+
+      template <int dim, typename VectorType>
+      bool
+      RefSpaceFEFieldFunction<dim, VectorType>::is_fe_q_iso_q1() const
+      {
+        return n_subdivisions_per_line != numbers::invalid_unsigned_int;
+      }
+
+
+
+      template <int dim, typename VectorType>
+      unsigned int
+      RefSpaceFEFieldFunction<dim, VectorType>::n_subdivisions() const
+      {
+        return n_subdivisions_per_line;
+      }
+
+
+
       template <int dim, typename VectorType>
       bool
       RefSpaceFEFieldFunction<dim, VectorType>::cell_is_set() const
@@ -1514,10 +1623,11 @@ namespace NonMatching
             // TODO: this could be extended to a component that is not zero
             return dealii::internal::evaluate_tensor_product_value(
               poly,
-              local_dof_values,
-              point,
+              this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+                                       local_dof_values,
+              this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point,
               polynomials_are_hat_functions,
-              renumber);
+              this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
           }
         else
           {
@@ -1546,10 +1656,13 @@ namespace NonMatching
             // 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,
+                     this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+                                              local_dof_values,
+                     this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) :
+                                              point,
                      polynomials_are_hat_functions,
-                     renumber)
+                     this->is_fe_q_iso_q1() ? std::vector<unsigned int>() :
+                                              renumber)
               .second;
           }
         else
@@ -1578,7 +1691,11 @@ namespace NonMatching
           {
             // 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);
+              poly,
+              this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+                                       local_dof_values,
+              this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point,
+              this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
           }
         else
           {
@@ -1591,6 +1708,57 @@ namespace NonMatching
             return symmetrize(hessian);
           }
       }
+
+
+
+      template <int dim>
+      BoundingBox<dim>
+      create_subcell_box(const BoundingBox<dim> &             unit_box,
+                         const std::array<unsigned int, dim> &subcell_indices,
+                         const unsigned int                   n_subdivisions)
+      {
+        std::pair<Point<dim>, Point<dim>> corners;
+        for (unsigned int d = 0; d < dim; ++d)
+          {
+            AssertIndexRange(subcell_indices[d], n_subdivisions + 1);
+
+            const double split_box_side_length =
+              (1. / n_subdivisions) * unit_box.side_length(d);
+            corners.first[d]  = subcell_indices[d] * split_box_side_length;
+            corners.second[d] = corners.first[d] + split_box_side_length;
+          }
+
+        return BoundingBox<dim>(corners);
+      }
+
+
+
+      template <int dim>
+      std::vector<unsigned int>
+      setup_lexicographic_mask(
+        const std::array<unsigned int, dim> &subcell_indices,
+        unsigned int                         dofs_per_line)
+      {
+        // This code expands the tensor product space, with `mask` increasing
+        // in size by a factor `2` for each iteration in `d`. Note how the
+        // indices in the next line are appended at the end of the vector.
+        std::vector<unsigned int> mask;
+        mask.reserve(1 << dim);
+        mask.push_back(0);
+        unsigned int stride = 1;
+        for (unsigned int d = 0; d < dim; ++d)
+          {
+            const unsigned int previous_mask_size = 1U << d;
+            AssertDimension(mask.size(), previous_mask_size);
+            for (unsigned int i = 0; i < previous_mask_size; ++i)
+              {
+                mask[i] += subcell_indices[d] * stride;
+                mask.push_back(stride + mask[i]);
+              }
+            stride *= dofs_per_line;
+          }
+        return mask;
+      }
     } // namespace DiscreteQuadratureGeneratorImplementation
   }   // namespace internal
 
@@ -1627,10 +1795,30 @@ namespace NonMatching
 
 
 
+  template <int dim>
+  void
+  QuadratureGenerator<dim>::clear_quadratures()
+  {
+    q_generator.clear_quadratures();
+  }
+
+
+
   template <int dim>
   void
   QuadratureGenerator<dim>::generate(const Function<dim> &   level_set,
                                      const BoundingBox<dim> &box)
+  {
+    clear_quadratures();
+    generate_append(level_set, box);
+  }
+
+
+
+  template <int dim>
+  void
+  QuadratureGenerator<dim>::generate_append(const Function<dim> &   level_set,
+                                            const BoundingBox<dim> &box)
   {
     Assert(level_set.n_components == 1,
            ExcMessage(
@@ -1638,8 +1826,6 @@ namespace NonMatching
              " it should have one component."));
     Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
 
-    q_generator.clear_quadratures();
-
     std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
     level_sets.push_back(level_set);
 
@@ -1703,11 +1889,33 @@ namespace NonMatching
 
 
 
+  template <int dim>
+  void
+  FaceQuadratureGenerator<dim>::clear_quadratures()
+  {
+    quadrature_generator.clear_quadratures();
+    surface_quadrature.clear();
+  }
+
+
+
   template <int dim>
   void
   FaceQuadratureGenerator<dim>::generate(const Function<dim> &   level_set,
                                          const BoundingBox<dim> &box,
                                          const unsigned int      face_index)
+  {
+    clear_quadratures();
+    generate_append(level_set, box, face_index);
+  }
+
+
+
+  template <int dim>
+  void
+  FaceQuadratureGenerator<dim>::generate_append(const Function<dim> &level_set,
+                                                const BoundingBox<dim> &box,
+                                                const unsigned int face_index)
   {
     AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
 
@@ -1727,7 +1935,8 @@ namespace NonMatching
     // Reuse the lower dimensional QuadratureGenerator on the face.
     const BoundingBox<dim - 1> cross_section =
       box.cross_section(face_normal_direction);
-    quadrature_generator.generate(face_restriction, cross_section);
+
+    quadrature_generator.generate_append(face_restriction, cross_section);
 
     // We need the dim-dimensional normals of the zero-contour.
     // Recompute these.
@@ -1801,10 +2010,29 @@ namespace NonMatching
 
 
 
+  void
+  FaceQuadratureGenerator<1>::clear_quadratures()
+  {
+    // do nothing
+  }
+
+
+
   void
   FaceQuadratureGenerator<1>::generate(const Function<1> &   level_set,
                                        const BoundingBox<1> &box,
                                        const unsigned int    face_index)
+  {
+    clear_quadratures();
+    generate_append(level_set, box, face_index);
+  }
+
+
+
+  void
+  FaceQuadratureGenerator<1>::generate_append(const Function<1> &   level_set,
+                                              const BoundingBox<1> &box,
+                                              const unsigned int    face_index)
   {
     AssertIndexRange(face_index, GeometryInfo<1>::faces_per_cell);
 
@@ -1885,6 +2113,63 @@ namespace NonMatching
 
 
 
+  template <int dim>
+  void
+  DiscreteQuadratureGenerator<dim>::generate_fe_q_iso_q1(
+    const BoundingBox<dim> &unit_box)
+  {
+    const unsigned int n_subdivisions =
+      reference_space_level_set->n_subdivisions();
+
+    const unsigned int dofs_per_line = n_subdivisions + 1;
+
+    this->clear_quadratures();
+
+    // The triple nested loop below expands the tensor product of the
+    // subdivisions, increasing the size of 'all_subcell_indices' by a factor
+    // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
+    // copies of the data are added at the end of the vector in each
+    // iteration over d. This code avoids a dimension-dependent
+    // number of nested loops over each direction, yet having the same
+    // effect.
+    std::vector<std::array<unsigned int, dim>> all_subcell_indices;
+    all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
+    all_subcell_indices.emplace_back();
+    for (unsigned d = 0; d < dim; ++d)
+      {
+        const unsigned int old_size = all_subcell_indices.size();
+        for (unsigned int i = 1; i < n_subdivisions; ++i)
+          {
+            for (unsigned int j = 0; j < old_size; ++j)
+              {
+                std::array<unsigned int, dim> next_entry =
+                  all_subcell_indices[j];
+                next_entry[d] = i;
+                all_subcell_indices.push_back(next_entry);
+              }
+          }
+      }
+
+
+    for (const auto &subcell_indices : all_subcell_indices)
+      {
+        const std::vector<unsigned int> lexicographic_mask =
+          internal::DiscreteQuadratureGeneratorImplementation::
+            setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
+
+        const auto subcell_box =
+          internal::DiscreteQuadratureGeneratorImplementation::
+            create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
+
+        reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
+
+        QuadratureGenerator<dim>::generate_append(*reference_space_level_set,
+                                                  subcell_box);
+      }
+  }
+
+
+
   template <int dim>
   void
   DiscreteQuadratureGenerator<dim>::generate(
@@ -1896,7 +2181,10 @@ namespace NonMatching
 
     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);
+    if (reference_space_level_set->is_fe_q_iso_q1())
+      generate_fe_q_iso_q1(unit_box);
+    else
+      QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
   }
 
 
@@ -1918,6 +2206,71 @@ namespace NonMatching
 
 
 
+  template <int dim>
+  void
+  DiscreteFaceQuadratureGenerator<dim>::generate_fe_q_iso_q1(
+    const BoundingBox<dim> &unit_box,
+    unsigned int            face_index)
+  {
+    const unsigned int n_subdivisions =
+      reference_space_level_set->n_subdivisions();
+
+    const unsigned int dofs_per_line = n_subdivisions + 1;
+
+    this->clear_quadratures();
+
+    const unsigned int coordinate_direction_face   = face_index / 2;
+    const unsigned int coordinate_orientation_face = face_index % 2;
+
+    // The triple nested loop below expands the tensor product of the
+    // subdivisions, increasing the size of 'all_subcell_indices' by a factor
+    // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
+    // copies of the data are added at the end of the vector in each
+    // iteration over d. This code avoids a dimension-dependent
+    // number of nested loops over each direction, yet having the same
+    // effect. To account for the face all subcell indices in the coordinate of
+    // the face direction are filled with the appropriate index depending on the
+    // coordinate orientation. The other two indices (on the face) are
+    // set accordingly in the nested loops (with a warp around).
+    std::vector<std::array<unsigned int, dim>> all_subcell_indices;
+    all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
+    all_subcell_indices.emplace_back();
+    all_subcell_indices[0][coordinate_direction_face] =
+      coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1;
+    for (unsigned d = 0; d < dim - 1; ++d)
+      {
+        const unsigned int old_size = all_subcell_indices.size();
+        for (unsigned int i = 1; i < n_subdivisions; ++i)
+          {
+            for (unsigned int j = 0; j < old_size; ++j)
+              {
+                std::array<unsigned int, dim> next_entry =
+                  all_subcell_indices[j];
+                next_entry[(coordinate_direction_face + d + 1) % dim] = i;
+                all_subcell_indices.push_back(next_entry);
+              }
+          }
+      }
+
+    for (const auto &subcell_indices : all_subcell_indices)
+      {
+        const std::vector<unsigned int> lexicographic_mask =
+          internal::DiscreteQuadratureGeneratorImplementation::
+            setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
+
+        const auto subcell_box =
+          internal::DiscreteQuadratureGeneratorImplementation::
+            create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
+
+        reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
+
+        FaceQuadratureGenerator<dim>::generate_append(
+          *reference_space_level_set, subcell_box, face_index);
+      }
+  }
+
+
+
   template <int dim>
   void
   DiscreteFaceQuadratureGenerator<dim>::generate(
@@ -1930,9 +2283,12 @@ namespace NonMatching
 
     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);
+    if (reference_space_level_set->is_fe_q_iso_q1())
+      generate_fe_q_iso_q1(unit_box, face_index);
+    else
+      FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
+                                             unit_box,
+                                             face_index);
   }
 } // namespace NonMatching
 #include "quadrature_generator.inst"
index 87d1be03ce93415686427c18913b669509040e11..09b54133721715c39fc8ce46466e6917782aca2e 100644 (file)
@@ -18,6 +18,7 @@
 #include <deal.II/base/quadrature_lib.h>
 
 #include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
 
 #include <deal.II/grid/grid_generator.h>
 
@@ -35,7 +36,7 @@ template <int dim>
 class Test
 {
 public:
-  Test();
+  Test(const FE_Poly<dim> &fe);
 
   void
   run();
@@ -67,11 +68,15 @@ private:
 
 
 template <int dim>
-Test<dim>::Test()
+Test<dim>::Test(const FE_Poly<dim> &fe)
   : dof_handler(triangulation)
 {
-  fe_collection.push_back(FE_Q<dim>(1));
-  const unsigned int n_quadrature_points = 1;
+  const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+    dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
+  const bool is_fe_q_iso_q1 = fe_q_iso_q1 != nullptr;
+  fe_collection.push_back(fe);
+
+  const unsigned int n_quadrature_points = is_fe_q_iso_q1 ? 1 : fe.get_degree();
   q_collection1D.push_back(QGauss<1>(n_quadrature_points));
 }
 
@@ -131,8 +136,11 @@ Test<dim>::test_discrete_quadrature_generator()
     q_collection1D, dof_handler, level_set);
   quadrature_generator.generate(triangulation.begin_active());
 
+  deallog << "inside" << std::endl;
   print_quadrature(quadrature_generator.get_inside_quadrature());
+  deallog << "outside" << std::endl;
   print_quadrature(quadrature_generator.get_outside_quadrature());
+  deallog << "surface" << std::endl;
   print_surface_quadrature(quadrature_generator.get_surface_quadrature());
 }
 
@@ -149,9 +157,12 @@ Test<dim>::test_discrete_face_quadrature_generator()
   for (const auto f : cell->face_indices())
     {
       face_quadrature_generator.generate(cell, f);
-
+      deallog << "face_index = " << f << std::endl;
+      deallog << "inside" << std::endl;
       print_quadrature(face_quadrature_generator.get_inside_quadrature());
+      deallog << "outside" << std::endl;
       print_quadrature(face_quadrature_generator.get_outside_quadrature());
+      deallog << "surface" << std::endl;
       print_surface_quadrature(
         face_quadrature_generator.get_surface_quadrature());
     }
@@ -161,10 +172,10 @@ Test<dim>::test_discrete_face_quadrature_generator()
 
 template <int dim>
 void
-run_test()
+run_test(const FE_Poly<dim> &fe)
 {
   deallog << "dim = " << dim << std::endl;
-  Test<dim> test;
+  Test<dim> test(fe);
   test.run();
   deallog << std::endl;
 }
@@ -176,7 +187,19 @@ main()
 {
   initlog();
 
-  run_test<1>();
-  run_test<2>();
-  run_test<3>();
+  std::shared_ptr<FE_Poly<1>> fe_1;
+  std::shared_ptr<FE_Poly<2>> fe_2;
+  std::shared_ptr<FE_Poly<3>> fe_3;
+  fe_1 = std::make_shared<FE_Q<1>>(1);
+  fe_2 = std::make_shared<FE_Q<2>>(1);
+  fe_3 = std::make_shared<FE_Q<3>>(1);
+  run_test<1>(*fe_1);
+  run_test<2>(*fe_2);
+  run_test<3>(*fe_3);
+  fe_1 = std::make_shared<FE_Q_iso_Q1<1>>(2);
+  fe_2 = std::make_shared<FE_Q_iso_Q1<2>>(2);
+  fe_3 = std::make_shared<FE_Q_iso_Q1<3>>(2);
+  run_test<1>(*fe_1);
+  run_test<2>(*fe_2);
+  run_test<3>(*fe_3);
 }
index 4dfaa1998f76d8be3d4ec938ce3489194069eec2..59b2760ca4ff0f88f81729679c8c554ad81ef665 100644 (file)
 
 DEAL::dim = 1
+DEAL::inside
 DEAL::0.375000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.250000
+DEAL::surface
 DEAL::0.750000, 1.00000, 1.00000
+DEAL::face_index = 0
+DEAL::inside
 DEAL::1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
 DEAL::1.00000
+DEAL::surface
 DEAL::
 DEAL::dim = 2
+DEAL::inside
 DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
 DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
 DEAL::0.500000, 1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
 DEAL::0.500000, 1.00000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
 DEAL::0.375000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.250000
+DEAL::surface
 DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
 DEAL::0.375000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.250000
+DEAL::surface
 DEAL::0.750000, 1.00000, 1.00000, 0.00000
 DEAL::
 DEAL::dim = 3
+DEAL::inside
 DEAL::0.375000, 0.500000, 0.500000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.500000, 0.500000, 0.250000
+DEAL::surface
 DEAL::0.750000, 0.500000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
 DEAL::0.500000, 0.500000, 1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
 DEAL::0.500000, 0.500000, 1.00000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
 DEAL::0.500000, 0.375000, 0.750000
+DEAL::outside
 DEAL::0.500000, 0.875000, 0.250000
+DEAL::surface
 DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
 DEAL::0.500000, 0.375000, 0.750000
+DEAL::outside
 DEAL::0.500000, 0.875000, 0.250000
+DEAL::surface
 DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 4
+DEAL::inside
 DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
 DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 5
+DEAL::inside
 DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
 DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
 DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
 DEAL::
+DEAL::dim = 1
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::1.00000
+DEAL::surface
+DEAL::
+DEAL::dim = 2
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.750000, 0.500000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::0.250000, 0.500000
+DEAL::0.750000, 0.500000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::
+DEAL::dim = 3
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000, 0.125000
+DEAL::0.625000, 0.250000, 0.250000, 0.0625000
+DEAL::0.250000, 0.750000, 0.250000, 0.125000
+DEAL::0.625000, 0.750000, 0.250000, 0.0625000
+DEAL::0.250000, 0.250000, 0.750000, 0.125000
+DEAL::0.625000, 0.250000, 0.750000, 0.0625000
+DEAL::0.250000, 0.750000, 0.750000, 0.125000
+DEAL::0.625000, 0.750000, 0.750000, 0.0625000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.250000, 0.0625000
+DEAL::0.875000, 0.750000, 0.250000, 0.0625000
+DEAL::0.875000, 0.250000, 0.750000, 0.0625000
+DEAL::0.875000, 0.750000, 0.750000, 0.0625000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.250000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.750000, 0.750000, 0.250000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.750000, 0.750000, 0.250000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.625000, 0.125000
+DEAL::0.750000, 0.625000, 0.125000
+DEAL::outside
+DEAL::0.250000, 0.875000, 0.125000
+DEAL::0.750000, 0.875000, 0.125000
+DEAL::surface
+DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.625000, 0.125000
+DEAL::0.750000, 0.625000, 0.125000
+DEAL::outside
+DEAL::0.250000, 0.875000, 0.125000
+DEAL::0.750000, 0.875000, 0.125000
+DEAL::surface
+DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 4
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 5
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 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.