Merge pull request #16932 from peterrum/nm_fevalues_cellaccessor

NonMatching::FEValues: accept CellAccessor

include/deal.II/non_matching/fe_values.h

@@ -216,14 +216,57 @@ namespace NonMatching
      * 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_inside_fe_values(), get_outside_fe_values(), or
      * get_surface_fe_values().
+     *
+     * If the @p q_index argument is left at its default value, then we use
+     * that quadrature formula within the hp::QCollection passed to the
+     * constructor of this class with index given by
+     * <code>cell-@>active_fe_index()</code>, i.e. the same index as that of
+     * the finite element. In this case, there should be a corresponding
+     * quadrature formula for each finite element in the hp::FECollection. As
+     * a special case, if the quadrature collection contains only a single
+     * element (a frequent case if one wants to use the same quadrature object
+     * for all finite elements in an hp-discretization, even if that may not
+     * be the most efficient), then this single quadrature is used unless a
+     * different value for this argument is specified. On the other hand, if a
+     * value is given for this argument, it overrides the choice of
+     * <code>cell-@>active_fe_index()</code> or the choice for the single
+     * quadrature.
+     *
+     * If the @p mapping_index argument is left at its default value, then we
+     * use that mapping object within the hp::MappingCollection passed to the
+     * constructor of this class with index given by
+     * <code>cell-@>active_fe_index()</code>, i.e. the same index as that of
+     * the finite element. As above, if the mapping collection contains only a
+     * single element (a frequent case if one wants to use a $Q_1$ mapping for
+     * all finite elements in an hp-discretization), then this single mapping
+     * is used unless a different value for this argument is specified.
      */
     template <bool level_dof_access>
     void
     reinit(
-      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
+      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+      const unsigned int q_index       = numbers::invalid_unsigned_int,
+      const unsigned int mapping_index = numbers::invalid_unsigned_int);
+
+    /**
+     * Like the previous function, but for non-DoFHandler iterators. The reason
+     * this function exists is so that one can use NonMatching::FEValues for
+     * Triangulation objects too.
+     *
+     * Since <code>cell-@>active_fe_index()</code> doesn't make sense for
+     * triangulation iterators, this function chooses the zero-th finite
+     * element, mapping, and quadrature object from the relevant constructions
+     * passed to the constructor of this object. The only exception is if you
+     * specify a value different from the default value for any of these last
+     * three arguments.
+     */
+    void
+    reinit(const TriaIterator<CellAccessor<dim, dim>> &cell,
+           const unsigned int q_index       = numbers::invalid_unsigned_int,
+           const unsigned int mapping_index = numbers::invalid_unsigned_int,
+           const unsigned int fe_index      = numbers::invalid_unsigned_int);
 
     /**
      * Return an dealii::FEValues object reinitialized with a quadrature for the
@@ -258,6 +301,16 @@ namespace NonMatching
     get_surface_fe_values() const;
 
   private:
+    /**
+     * Internal function called by the reinit() functions.
+     */
+    template <typename CellIteratorType>
+    void
+    reinit_internal(const CellIteratorType &cell,
+                    const unsigned int      q_index,
+                    const unsigned int      mapping_index,
+                    const unsigned int      fe_index);
+
     /**
      * Do work common to the constructors. The incoming QCollection should be
      * quadratures integrating over $[0, 1]^{dim}$. These will be used on the

source/non_matching/fe_values.cc

@@ -148,10 +148,72 @@ namespace NonMatching
   template <bool level_dof_access>
   void
   FEValues<dim>::reinit(
-    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
+    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+    const unsigned int                                               q_index,
+    const unsigned int mapping_index)
+  {
+    this->reinit_internal(cell,
+                          q_index,
+                          mapping_index,
+                          cell->active_fe_index());
+  }
+
+
+
+  template <int dim>
+  void
+  FEValues<dim>::reinit(const TriaIterator<CellAccessor<dim, dim>> &cell,
+                        const unsigned int                          q_index,
+                        const unsigned int mapping_index,
+                        const unsigned int fe_index)
+  {
+    this->reinit_internal(cell, q_index, mapping_index, fe_index);
+  }
+
+
+
+  template <int dim>
+  template <typename CellIteratorType>
+  void
+  FEValues<dim>::reinit_internal(const CellIteratorType &cell,
+                                 const unsigned int      q_index_in,
+                                 const unsigned int      mapping_index_in,
+                                 const unsigned int      fe_index_in)
   {
     current_cell_location = mesh_classifier->location_to_level_set(cell);
-    active_fe_index       = cell->active_fe_index();
+
+    if (fe_index_in == numbers::invalid_unsigned_int)
+      this->active_fe_index = 0;
+    else
+      this->active_fe_index = fe_index_in;
+
+    unsigned int mapping_index = mapping_index_in;
+    unsigned int q_index       = q_index_in;
+    unsigned int q_index_1D    = q_index_in;
+
+    if (mapping_index == numbers::invalid_unsigned_int)
+      {
+        if (mapping_collection->size() > 1)
+          mapping_index = active_fe_index;
+        else
+          mapping_index = 0;
+      }
+
+    if (q_index == numbers::invalid_unsigned_int)
+      {
+        if (fe_values_inside_full_quadrature.size() > 1)
+          q_index = active_fe_index;
+        else
+          q_index = 0;
+      }
+
+    if (q_index_1D == numbers::invalid_unsigned_int)
+      {
+        if (q_collection_1D.size() > 1)
+          q_index_1D = active_fe_index;
+        else
+          q_index_1D = 0;
+      }
 
     // These objects were created with a quadrature based on the previous cell
     // and are thus no longer valid.
@@ -163,22 +225,29 @@ namespace NonMatching
       {
         case LocationToLevelSet::inside:
           {
-            fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell);
+            Assert((active_fe_index == mapping_index) ||
+                     ((mapping_collection->size() == 1) &&
+                      (mapping_index == 0)),
+                   ExcNotImplemented());
+            Assert(active_fe_index == q_index, ExcNotImplemented());
+
+            fe_values_inside_full_quadrature.at(q_index)->reinit(cell);
             break;
           }
         case LocationToLevelSet::outside:
           {
-            fe_values_outside_full_quadrature.at(active_fe_index)->reinit(cell);
+            Assert((active_fe_index == mapping_index) ||
+                     ((mapping_collection->size() == 1) &&
+                      (mapping_index == 0)),
+                   ExcNotImplemented());
+            Assert(active_fe_index == q_index, ExcNotImplemented());
+
+            fe_values_outside_full_quadrature.at(q_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);
+            quadrature_generator.set_1D_quadrature(q_index_1D);
             quadrature_generator.generate(cell);
 
             const Quadrature<dim> &inside_quadrature =

source/non_matching/fe_values.inst.in

@@ -69,7 +69,9 @@ for (deal_II_dimension : DIMENSIONS; lda : BOOL)
   {
     template void FEValues<deal_II_dimension>::reinit(
       const TriaIterator<
-        DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+        DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
+      const unsigned int,
+      const unsigned int);
 
     template void FEInterfaceValues<deal_II_dimension>::do_reinit(
       const TriaIterator<

tests/non_matching/fe_values.cc

@@ -88,9 +88,10 @@ class Test
   void
   setup_discrete_level_set();
 
+  template <typename IteratorType>
   void
-  test_fe_values_reinitializes_correctly(
-    NonMatching::FEValues<dim> &fe_values) const;
+  test_fe_values_reinitializes_correctly(NonMatching::FEValues<dim> &fe_values,
+                                         IteratorType cell) const;
 
   Triangulation<dim>    triangulation;
   hp::FECollection<dim> fe_collection;
@@ -139,7 +140,10 @@ Test<dim>::run()
                                          mesh_classifier,
                                          dof_handler,
                                          level_set);
-    test_fe_values_reinitializes_correctly(fe_values);
+    test_fe_values_reinitializes_correctly(fe_values,
+                                           triangulation.begin_active());
+    test_fe_values_reinitializes_correctly(fe_values,
+                                           dof_handler.begin_active());
   }
   {
     // Test with the "more advanced" constructor.
@@ -151,7 +155,10 @@ Test<dim>::run()
                                          mesh_classifier,
                                          dof_handler,
                                          level_set);
-    test_fe_values_reinitializes_correctly(fe_values);
+    test_fe_values_reinitializes_correctly(fe_values,
+                                           triangulation.begin_active());
+    test_fe_values_reinitializes_correctly(fe_values,
+                                           dof_handler.begin_active());
   }
 }
 
@@ -197,13 +204,14 @@ Test<dim>::setup_discrete_level_set()
 
 
 template <int dim>
+template <typename IteratorType>
 void
 Test<dim>::test_fe_values_reinitializes_correctly(
-  NonMatching::FEValues<dim> &fe_values) const
+  NonMatching::FEValues<dim> &fe_values,
+  IteratorType                cell) 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,

tests/non_matching/fe_values.with_lapack=on.output

@@ -18,6 +18,24 @@ 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::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
 DEAL::
 DEAL::dim = 2
 DEAL::OK
@@ -38,6 +56,24 @@ 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::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
 DEAL::
 DEAL::dim = 3
 DEAL::OK
@@ -58,4 +94,22 @@ 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::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
 DEAL::