]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Changes to FEInterfaceValues and NM::FEInterfaceValues 17045/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 19 May 2024 07:47:05 +0000 (09:47 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 20 Jun 2024 07:49:24 +0000 (09:49 +0200)
include/deal.II/fe/fe_interface_values.h
include/deal.II/non_matching/fe_values.h
source/non_matching/fe_values.cc
source/non_matching/fe_values.inst.in
tests/non_matching/fe_interface_values.cc
tests/non_matching/fe_interface_values.with_lapack=on.output

index 3206af44b454271e0fe1637e64c038edfea6d85b..359e99555c76a78b65c2e75446a262437207bf6c 100644 (file)
@@ -1441,6 +1441,8 @@ public:
    *   finite element across the interface (only used if the FEInterface object
    *   is initialized with an hp::FECollection, an hp::QCollection, and possibly
    *   an hp::MappingCollection).
+   * @param[in] fe_index_neighbor Active fe index of neighboring cell. Useful
+   *   if hp capabilities are used and for non-DoFHandler iterators.
    */
   template <typename CellIteratorType, typename CellNeighborIteratorType>
   void
@@ -1450,9 +1452,10 @@ public:
          const CellNeighborIteratorType &cell_neighbor,
          const unsigned int              face_no_neighbor,
          const unsigned int              sub_face_no_neighbor,
-         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);
+         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,
+         const unsigned int fe_index_neighbor = numbers::invalid_unsigned_int);
 
   /**
    * Re-initialize this object to be used on an interface given by a single face
@@ -2318,11 +2321,24 @@ FEInterfaceValues<dim, spacedim>::reinit(
   const unsigned int              sub_face_no_neighbor,
   const unsigned int              q_index,
   const unsigned int              mapping_index,
-  const unsigned int              fe_index)
+  const unsigned int              fe_index_in,
+  const unsigned int              fe_index_neighbor_in)
 {
   Assert(internal_fe_face_values || internal_hp_fe_face_values,
          ExcNotInitialized());
 
+  constexpr bool is_dof_cell_accessor =
+    std::is_same_v<DoFCellAccessor<dim, spacedim, true>,
+                   typename CellIteratorType::AccessorType> ||
+    std::is_same_v<DoFCellAccessor<dim, spacedim, false>,
+                   typename CellIteratorType::AccessorType>;
+
+  constexpr bool is_dof_cell_accessor_neighbor =
+    std::is_same_v<DoFCellAccessor<dim, spacedim, true>,
+                   typename CellNeighborIteratorType::AccessorType> ||
+    std::is_same_v<DoFCellAccessor<dim, spacedim, false>,
+                   typename CellNeighborIteratorType::AccessorType>;
+
   if (internal_fe_face_values)
     {
       if (sub_face_no == numbers::invalid_unsigned_int)
@@ -2357,6 +2373,28 @@ FEInterfaceValues<dim, spacedim>::reinit(
     }
   else if (internal_hp_fe_face_values)
     {
+      unsigned int active_fe_index = fe_index_in;
+      unsigned int active_fe_index_neighbor =
+        (fe_index_neighbor_in != numbers::invalid_unsigned_int) ?
+          fe_index_neighbor_in :
+          fe_index_in;
+
+      if (active_fe_index == numbers::invalid_unsigned_int)
+        {
+          if constexpr (is_dof_cell_accessor)
+            active_fe_index = cell->active_fe_index();
+          else
+            active_fe_index = 0;
+        }
+
+      if (active_fe_index_neighbor == numbers::invalid_unsigned_int)
+        {
+          if constexpr (is_dof_cell_accessor_neighbor)
+            active_fe_index_neighbor = cell_neighbor->active_fe_index();
+          else
+            active_fe_index_neighbor = 0;
+        }
+
       unsigned int used_q_index       = q_index;
       unsigned int used_mapping_index = mapping_index;
 
@@ -2376,10 +2414,10 @@ FEInterfaceValues<dim, spacedim>::reinit(
       // same :-(
       if (used_q_index == numbers::invalid_unsigned_int)
         if (internal_hp_fe_face_values
-              ->get_quadrature_collection()[cell->active_fe_index()] ==
+              ->get_quadrature_collection()[active_fe_index] ==
             internal_hp_fe_face_values
-              ->get_quadrature_collection()[cell_neighbor->active_fe_index()])
-          used_q_index = cell->active_fe_index();
+              ->get_quadrature_collection()[active_fe_index_neighbor])
+          used_q_index = active_fe_index;
 
       // Third check, if the above did not already suffice. We see if we
       // can get somewhere via the dominated's finite element index.
@@ -2387,7 +2425,7 @@ FEInterfaceValues<dim, spacedim>::reinit(
         ((used_q_index == numbers::invalid_unsigned_int) ||
              (used_mapping_index == numbers::invalid_unsigned_int) ?
            internal_hp_fe_face_values->get_fe_collection().find_dominated_fe(
-             {cell->active_fe_index(), cell_neighbor->active_fe_index()}) :
+             {active_fe_index, active_fe_index_neighbor}) :
            numbers::invalid_unsigned_int);
 
       if (used_q_index == numbers::invalid_unsigned_int)
@@ -2420,14 +2458,18 @@ FEInterfaceValues<dim, spacedim>::reinit(
       if (sub_face_no == numbers::invalid_unsigned_int)
         {
           internal_hp_fe_face_values->reinit(
-            cell, face_no, used_q_index, used_mapping_index, fe_index);
+            cell, face_no, used_q_index, used_mapping_index, active_fe_index);
           fe_face_values = &const_cast<FEFaceValues<dim, spacedim> &>(
             internal_hp_fe_face_values->get_present_fe_values());
         }
       else
         {
-          internal_hp_fe_subface_values->reinit(
-            cell, face_no, sub_face_no, used_q_index, used_mapping_index);
+          internal_hp_fe_subface_values->reinit(cell,
+                                                face_no,
+                                                sub_face_no,
+                                                used_q_index,
+                                                used_mapping_index,
+                                                active_fe_index);
 
           fe_face_values = &const_cast<FESubfaceValues<dim, spacedim> &>(
             internal_hp_fe_subface_values->get_present_fe_values());
@@ -2437,18 +2479,21 @@ FEInterfaceValues<dim, spacedim>::reinit(
           internal_hp_fe_face_values_neighbor->reinit(cell_neighbor,
                                                       face_no_neighbor,
                                                       used_q_index,
-                                                      used_mapping_index);
+                                                      used_mapping_index,
+                                                      active_fe_index_neighbor);
 
           fe_face_values_neighbor = &const_cast<FEFaceValues<dim, spacedim> &>(
             internal_hp_fe_face_values_neighbor->get_present_fe_values());
         }
       else
         {
-          internal_hp_fe_subface_values_neighbor->reinit(cell_neighbor,
-                                                         face_no_neighbor,
-                                                         sub_face_no_neighbor,
-                                                         used_q_index,
-                                                         used_mapping_index);
+          internal_hp_fe_subface_values_neighbor->reinit(
+            cell_neighbor,
+            face_no_neighbor,
+            sub_face_no_neighbor,
+            used_q_index,
+            used_mapping_index,
+            active_fe_index_neighbor);
 
           fe_face_values_neighbor =
             &const_cast<FESubfaceValues<dim, spacedim> &>(
@@ -2463,47 +2508,48 @@ FEInterfaceValues<dim, spacedim>::reinit(
     }
 
   // Set up dof mapping and remove duplicates (for continuous elements).
-  {
-    // Get dof indices first:
-    std::vector<types::global_dof_index> v(
-      fe_face_values->get_fe().n_dofs_per_cell());
-    cell->get_active_or_mg_dof_indices(v);
-    std::vector<types::global_dof_index> v2(
-      fe_face_values_neighbor->get_fe().n_dofs_per_cell());
-    cell_neighbor->get_active_or_mg_dof_indices(v2);
-
-    // Fill a map from the global dof index to the left and right
-    // local index.
-    std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
-                                          tempmap;
-    std::pair<unsigned int, unsigned int> invalid_entry(
-      numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);
-
-    for (unsigned int i = 0; i < v.size(); ++i)
-      {
-        // If not already existing, add an invalid entry:
-        auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
-        result.first->second.first = i;
-      }
-
-    for (unsigned int i = 0; i < v2.size(); ++i)
-      {
-        // If not already existing, add an invalid entry:
-        auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
-        result.first->second.second = i;
-      }
-
-    // Transfer from the map to the sorted std::vectors.
-    dofmap.resize(tempmap.size());
-    interface_dof_indices.resize(tempmap.size());
-    unsigned int idx = 0;
-    for (auto &x : tempmap)
-      {
-        interface_dof_indices[idx] = x.first;
-        dofmap[idx]                = {{x.second.first, x.second.second}};
-        ++idx;
-      }
-  }
+  if constexpr (is_dof_cell_accessor_neighbor && is_dof_cell_accessor)
+    {
+      // Get dof indices first:
+      std::vector<types::global_dof_index> v(
+        fe_face_values->get_fe().n_dofs_per_cell());
+      cell->get_active_or_mg_dof_indices(v);
+      std::vector<types::global_dof_index> v2(
+        fe_face_values_neighbor->get_fe().n_dofs_per_cell());
+      cell_neighbor->get_active_or_mg_dof_indices(v2);
+
+      // Fill a map from the global dof index to the left and right
+      // local index.
+      std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
+                                            tempmap;
+      std::pair<unsigned int, unsigned int> invalid_entry(
+        numbers::invalid_unsigned_int, numbers::invalid_unsigned_int);
+
+      for (unsigned int i = 0; i < v.size(); ++i)
+        {
+          // If not already existing, add an invalid entry:
+          auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
+          result.first->second.first = i;
+        }
+
+      for (unsigned int i = 0; i < v2.size(); ++i)
+        {
+          // If not already existing, add an invalid entry:
+          auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
+          result.first->second.second = i;
+        }
+
+      // Transfer from the map to the sorted std::vectors.
+      dofmap.resize(tempmap.size());
+      interface_dof_indices.resize(tempmap.size());
+      unsigned int idx = 0;
+      for (auto &x : tempmap)
+        {
+          interface_dof_indices[idx] = x.first;
+          dofmap[idx]                = {{x.second.first, x.second.second}};
+          ++idx;
+        }
+    }
 }
 
 
@@ -2522,12 +2568,17 @@ FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
 
   if (internal_fe_face_values)
     {
+      Assert((q_index == 0 || q_index == numbers::invalid_unsigned_int),
+             ExcNotImplemented());
+      Assert((mapping_index == 0 ||
+              mapping_index == numbers::invalid_unsigned_int),
+             ExcNotImplemented());
+      Assert((fe_index == 0 || fe_index == numbers::invalid_unsigned_int),
+             ExcNotImplemented());
+
       internal_fe_face_values->reinit(cell, face_no);
       fe_face_values          = internal_fe_face_values.get();
       fe_face_values_neighbor = nullptr;
-
-      interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
-      cell->get_active_or_mg_dof_indices(interface_dof_indices);
     }
   else if (internal_hp_fe_face_values)
     {
@@ -2536,15 +2587,31 @@ FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
       fe_face_values = &const_cast<FEFaceValues<dim> &>(
         internal_hp_fe_face_values->get_present_fe_values());
       fe_face_values_neighbor = nullptr;
-
-      interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
-      cell->get_active_or_mg_dof_indices(interface_dof_indices);
     }
 
-  dofmap.resize(interface_dof_indices.size());
-  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
+  if constexpr (std::is_same_v<typename CellIteratorType::AccessorType,
+                               DoFCellAccessor<dim, spacedim, true>> ||
+                std::is_same_v<typename CellIteratorType::AccessorType,
+                               DoFCellAccessor<dim, spacedim, false>>)
     {
-      dofmap[i] = {{i, numbers::invalid_unsigned_int}};
+      if (internal_fe_face_values)
+        {
+          interface_dof_indices.resize(
+            fe_face_values->get_fe().n_dofs_per_cell());
+          cell->get_active_or_mg_dof_indices(interface_dof_indices);
+        }
+      else if (internal_hp_fe_face_values)
+        {
+          interface_dof_indices.resize(
+            fe_face_values->get_fe().n_dofs_per_cell());
+          cell->get_active_or_mg_dof_indices(interface_dof_indices);
+        }
+
+      dofmap.resize(interface_dof_indices.size());
+      for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
+        {
+          dofmap[i] = {{i, numbers::invalid_unsigned_int}};
+        }
     }
 }
 
index 6a014fba8ad060ef5ecca7a48d0d070d99c6ab5e..7d66598757bcb7f658fa8150366e805a3ab9997d 100644 (file)
@@ -578,7 +578,10 @@ namespace NonMatching
            const unsigned int              sub_face_no,
            const CellNeighborIteratorType &cell_neighbor,
            const unsigned int              face_no_neighbor,
-           const unsigned int              sub_face_no_neighbor);
+           const unsigned int              sub_face_no_neighbor,
+           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);
 
 
     /**
@@ -590,7 +593,11 @@ namespace NonMatching
      */
     template <typename CellIteratorType>
     void
-    reinit(const CellIteratorType &cell, const unsigned int face_no);
+    reinit(const CellIteratorType &cell,
+           const unsigned int      face_no,
+           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::FEInterfaceValues object reinitialized with a
@@ -631,12 +638,14 @@ namespace NonMatching
      * reinit on a single dealii::FEInterfaceValues object, which is what
      * differs between the two reinit functions.
      */
-    template <bool level_dof_access>
+    template <typename CellAccessorType>
     void
-    do_reinit(
-      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
-      const unsigned int                                               face_no,
-      const std::function<void(dealii::FEInterfaceValues<dim> &)> &call_reinit);
+    do_reinit(const TriaIterator<CellAccessorType>          &cell,
+              const unsigned int                             face_no,
+              const unsigned int                             q_index,
+              const unsigned int                             active_fe_index,
+              const std::function<void(dealii::FEInterfaceValues<dim> &,
+                                       const unsigned int)> &call_reinit);
 
     /**
      * A pointer to the collection of mappings to be used.
@@ -659,11 +668,6 @@ namespace NonMatching
      */
     LocationToLevelSet current_face_location;
 
-    /**
-     * Active fe index of the last cell that reinit was called with.
-     */
-    unsigned int active_fe_index;
-
     /**
      * The update flags passed to the constructor.
      */
@@ -675,38 +679,17 @@ namespace NonMatching
     const SmartPointer<const MeshClassifier<dim>> mesh_classifier;
 
     /**
-     * For each element in the FECollection passed to the constructor,
-     * this object contains an dealii::FEInterfaceValues object created with a
-     * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and
-     * UpdateFlags for the inside region. Thus, these optionals should always
-     * contain a value.
-     *
-     * When LocationToLevelSet of the cell is INSIDE (and we do not need
-     * to generate an immersed quadrature), we return the
-     * dealii::FEInterfaceValues object in this container corresponding to the
-     * cell's active_fe_index.
-     *
-     * This container is a std::deque, which is compatible with the
-     * `FEInterfaceValues` class that does not have a copy-constructor.
+     * FEInterfaceValues corresponding to cells with LocationToLevelSet
+     * INSIDE (not needed to generate an immersed quadrature).
      */
-    std::deque<std::optional<dealii::FEInterfaceValues<dim>>>
+    std::optional<dealii::FEInterfaceValues<dim>>
       fe_values_inside_full_quadrature;
 
     /**
-     * For each element in the FECollection passed to the constructor,
-     * this object contains an dealii::FEInterfaceValues object created with a
-     * quadrature rule over the full reference cell: $[0, 1]^{dim-1}$ and
-     * UpdateFlags for the outside region. Thus, these optionals should always
-     * contain a value.
-     *
-     * When LocationToLevelSet of the cell is OUTSIDE (and we do not need
-     * to generate an immersed quadrature), we return the dealii::FEValues
-     * object in this container corresponding to the cell's active_fe_index.
-     *
-     * This container is a std::deque, which is compatible with the
-     * `FEInterfaceValues` class that does not have a copy-constructor.
+     * FEInterfaceValues corresponding to cells with LocationToLevelSet
+     * OUTSIDE (not needed to generate an immersed quadrature).
      */
-    std::deque<std::optional<dealii::FEInterfaceValues<dim>>>
+    std::optional<dealii::FEInterfaceValues<dim>>
       fe_values_outside_full_quadrature;
 
     /**
@@ -742,16 +725,21 @@ namespace NonMatching
   template <typename CellIteratorType>
   inline void
   FEInterfaceValues<dim>::reinit(const CellIteratorType &cell,
-                                 const unsigned int      face_no)
+                                 const unsigned int      face_no,
+                                 const unsigned int      q_index,
+                                 const unsigned int      mapping_index,
+                                 const unsigned int      active_fe_index)
   {
     // Lambda describing how we should call reinit on a single
     // dealii::FEInterfaceValues object.
     const auto reinit_operation =
-      [&cell, face_no](dealii::FEInterfaceValues<dim> &fe_interface_values) {
-        fe_interface_values.reinit(cell, face_no);
+      [&](dealii::FEInterfaceValues<dim> &fe_interface_values,
+          const unsigned int              q_index) {
+        fe_interface_values.reinit(
+          cell, face_no, q_index, mapping_index, active_fe_index);
       };
 
-    do_reinit(cell, face_no, reinit_operation);
+    do_reinit(cell, face_no, q_index, active_fe_index, reinit_operation);
   }
 
 
@@ -764,7 +752,10 @@ namespace NonMatching
                                  const unsigned int              sub_face_no,
                                  const CellNeighborIteratorType &cell_neighbor,
                                  const unsigned int face_no_neighbor,
-                                 const unsigned int sub_face_no_neighbor)
+                                 const unsigned int sub_face_no_neighbor,
+                                 const unsigned int q_index,
+                                 const unsigned int mapping_index,
+                                 const unsigned int active_fe_index)
   {
     Assert(sub_face_no == numbers::invalid_unsigned_int, ExcNotImplemented());
     Assert(sub_face_no_neighbor == numbers::invalid_unsigned_int,
@@ -773,22 +764,20 @@ namespace NonMatching
     // Lambda describing how we should call reinit on a single
     // dealii::FEInterfaceValues object.
     const auto reinit_operation =
-      [&cell,
-       face_no,
-       sub_face_no,
-       &cell_neighbor,
-       face_no_neighbor,
-       sub_face_no_neighbor](
-        dealii::FEInterfaceValues<dim> &fe_interface_values) {
+      [&](dealii::FEInterfaceValues<dim> &fe_interface_values,
+          const unsigned int              q_index) {
         fe_interface_values.reinit(cell,
                                    face_no,
                                    sub_face_no,
                                    cell_neighbor,
                                    face_no_neighbor,
-                                   sub_face_no_neighbor);
+                                   sub_face_no_neighbor,
+                                   q_index,
+                                   mapping_index,
+                                   active_fe_index);
       };
 
-    do_reinit(cell, face_no, reinit_operation);
+    do_reinit(cell, face_no, q_index, active_fe_index, reinit_operation);
   }
 
 #endif
index 56675f15ce4165721d10d8cfa6f2a5ddd5b2a680..7e190e334a2a6a942961924a5b8edbb3842b4af3 100644 (file)
@@ -398,7 +398,6 @@ namespace NonMatching
     const hp::QCollection<dim - 1> &q_collection)
   {
     current_face_location = LocationToLevelSet::unassigned;
-    active_fe_index       = numbers::invalid_unsigned_int;
 
     Assert(fe_collection->size() > 0,
            ExcMessage("Incoming hp::FECollection can not be empty."));
@@ -417,43 +416,31 @@ namespace NonMatching
       ExcMessage(
         "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
 
-    // For each element in fe_collection, create dealii::FEInterfaceValues
-    // objects to use on the non-intersected cells.
-    fe_values_inside_full_quadrature.resize(fe_collection->size());
-    fe_values_outside_full_quadrature.resize(fe_collection->size());
-    for (unsigned int fe_index = 0; fe_index < fe_collection->size();
-         ++fe_index)
-      {
-        const unsigned int mapping_index =
-          mapping_collection->size() > 1 ? fe_index : 0;
-        const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
-
-        fe_values_inside_full_quadrature[fe_index].emplace(
-          (*mapping_collection)[mapping_index],
-          (*fe_collection)[fe_index],
-          q_collection[q_index],
-          region_update_flags.inside);
-        fe_values_outside_full_quadrature[fe_index].emplace(
-          (*mapping_collection)[mapping_index],
-          (*fe_collection)[fe_index],
-          q_collection[q_index],
-          region_update_flags.outside);
-      }
+    fe_values_inside_full_quadrature.emplace(*mapping_collection,
+                                             *fe_collection,
+                                             q_collection,
+                                             region_update_flags.inside);
+    fe_values_outside_full_quadrature.emplace(*mapping_collection,
+                                              *fe_collection,
+                                              q_collection,
+                                              region_update_flags.outside);
   }
 
 
 
   template <int dim>
-  template <bool level_dof_access>
+  template <typename CellAccessorType>
   void
   FEInterfaceValues<dim>::do_reinit(
-    const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
-    const unsigned int                                               face_no,
-    const std::function<void(dealii::FEInterfaceValues<dim> &)> &call_reinit)
+    const TriaIterator<CellAccessorType>          &cell,
+    const unsigned int                             face_no,
+    const unsigned int                             q_index_in,
+    const unsigned int                             active_fe_index_in,
+    const std::function<void(dealii::FEInterfaceValues<dim> &,
+                             const unsigned int)> &call_reinit)
   {
     current_face_location =
       mesh_classifier->location_to_level_set(cell, face_no);
-    active_fe_index = cell->active_fe_index();
 
     // These objects were created with a quadrature based on the previous cell
     // and are thus no longer valid.
@@ -464,22 +451,44 @@ namespace NonMatching
       {
         case LocationToLevelSet::inside:
           {
-            call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index));
+            call_reinit(*fe_values_inside_full_quadrature, q_index_in);
             break;
           }
         case LocationToLevelSet::outside:
           {
-            call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index));
+            call_reinit(*fe_values_outside_full_quadrature, q_index_in);
             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;
+            unsigned int q_index = q_index_in;
+
+            if (q_index == numbers::invalid_unsigned_int)
+              {
+                unsigned int active_fe_index = active_fe_index_in;
+
+                if (active_fe_index == numbers::invalid_unsigned_int)
+                  {
+                    if constexpr (std::is_same_v<
+                                    DoFCellAccessor<dim, dim, true>,
+                                    CellAccessorType> ||
+                                  std::is_same_v<
+                                    DoFCellAccessor<dim, dim, false>,
+                                    CellAccessorType>)
+                      active_fe_index = cell->active_fe_index();
+                    else
+                      active_fe_index = 0;
+                  }
+
+                if (q_collection_1D.size() > 1)
+                  q_index = active_fe_index;
+                else
+                  q_index = 0;
+              }
+
+            AssertIndexRange(q_index, q_collection_1D.size());
 
-            face_quadrature_generator.set_1D_quadrature(q1D_index);
+            face_quadrature_generator.set_1D_quadrature(q_index);
             face_quadrature_generator.generate(cell, face_no);
 
             const Quadrature<dim - 1> &inside_quadrature =
@@ -492,22 +501,24 @@ namespace NonMatching
             // object if that is the case.
             if (inside_quadrature.size() > 0)
               {
-                fe_values_inside.emplace((*mapping_collection)[mapping_index],
-                                         (*fe_collection)[active_fe_index],
-                                         inside_quadrature,
+                fe_values_inside.emplace(*mapping_collection,
+                                         *fe_collection,
+                                         hp::QCollection<dim - 1>(
+                                           inside_quadrature),
                                          region_update_flags.inside);
 
-                call_reinit(*fe_values_inside);
+                call_reinit(*fe_values_inside, /*q_index=*/0);
               }
 
             if (outside_quadrature.size() > 0)
               {
-                fe_values_outside.emplace((*mapping_collection)[mapping_index],
-                                          (*fe_collection)[active_fe_index],
-                                          outside_quadrature,
+                fe_values_outside.emplace(*mapping_collection,
+                                          *fe_collection,
+                                          hp::QCollection<dim - 1>(
+                                            outside_quadrature),
                                           region_update_flags.outside);
 
-                call_reinit(*fe_values_outside);
+                call_reinit(*fe_values_outside, /*q_index=*/0);
               }
             break;
           }
@@ -526,7 +537,7 @@ namespace NonMatching
   FEInterfaceValues<dim>::get_inside_fe_values() const
   {
     if (current_face_location == LocationToLevelSet::inside)
-      return fe_values_inside_full_quadrature.at(active_fe_index);
+      return fe_values_inside_full_quadrature;
     else
       return fe_values_inside;
   }
@@ -538,7 +549,7 @@ namespace NonMatching
   FEInterfaceValues<dim>::get_outside_fe_values() const
   {
     if (current_face_location == LocationToLevelSet::outside)
-      return fe_values_outside_full_quadrature.at(active_fe_index);
+      return fe_values_outside_full_quadrature;
     else
       return fe_values_outside;
   }
index 31a3b1a3f5823d0171772160bb0c79f0b5c70447..ff168fba4d41ea92a0fa59450b33b8b53078c948 100644 (file)
@@ -65,6 +65,17 @@ for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
   }
 
 // Template reinit functions
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template void FEInterfaceValues<deal_II_dimension>::do_reinit(
+      const TriaIterator<CellAccessor<deal_II_dimension, deal_II_dimension>> &,
+      const unsigned int,
+      const unsigned int,
+      const unsigned int,
+      const std::function<void(dealii::FEInterfaceValues<deal_II_dimension> &,
+                               const unsigned int)> &);
+  }
+
 for (deal_II_dimension : DIMENSIONS; lda : BOOL)
   {
     template void FEValues<deal_II_dimension>::reinit(
@@ -77,6 +88,8 @@ for (deal_II_dimension : DIMENSIONS; lda : BOOL)
       const TriaIterator<
         DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
       const unsigned int,
-      const std::function<void(dealii::FEInterfaceValues<deal_II_dimension> &)>
-        &);
+      const unsigned int,
+      const unsigned int,
+      const std::function<void(dealii::FEInterfaceValues<deal_II_dimension> &,
+                               const unsigned int)> &);
   }
index 1f100aa288be7216c58381584e605de1b4a391ac..6c6d8f40c57c27948860aa30bdc6633118601ad8 100644 (file)
@@ -55,9 +55,11 @@ private:
   void
   setup_discrete_level_set();
 
+  template <typename IteratorType>
   void
   print_which_optionals_have_values_on_cell_0(
-    NonMatching::FEInterfaceValues<dim> &fe_values);
+    NonMatching::FEInterfaceValues<dim> &fe_values,
+    IteratorType                         cell);
 
   Triangulation<dim>    triangulation;
   hp::FECollection<dim> fe_collection;
@@ -106,7 +108,10 @@ Test<dim>::run()
                                                   mesh_classifier,
                                                   dof_handler,
                                                   level_set);
-    print_which_optionals_have_values_on_cell_0(fe_values);
+    print_which_optionals_have_values_on_cell_0(fe_values,
+                                                triangulation.begin_active());
+    print_which_optionals_have_values_on_cell_0(fe_values,
+                                                dof_handler.begin_active());
   }
   {
     deallog << "Advanced constructor:" << std::endl;
@@ -118,7 +123,10 @@ Test<dim>::run()
                                                   mesh_classifier,
                                                   dof_handler,
                                                   level_set);
-    print_which_optionals_have_values_on_cell_0(fe_values);
+    print_which_optionals_have_values_on_cell_0(fe_values,
+                                                triangulation.begin_active());
+    print_which_optionals_have_values_on_cell_0(fe_values,
+                                                dof_handler.begin_active());
   }
 }
 
@@ -166,12 +174,12 @@ Test<dim>::setup_discrete_level_set()
 
 
 template <int dim>
+template <typename IteratorType>
 void
 Test<dim>::print_which_optionals_have_values_on_cell_0(
-  NonMatching::FEInterfaceValues<dim> &fe_values)
+  NonMatching::FEInterfaceValues<dim> &fe_values,
+  IteratorType                         cell)
 {
-  const auto cell = dof_handler.begin_active();
-
   for (const unsigned int face_index : cell->face_indices())
     {
       deallog << "face " << face_index << std::endl;
index 52bca2f64357dae2e53cad64e2fc2526334407be..84666ddd52c65162fdb71e28327001c3069d0a1d 100644 (file)
@@ -5,11 +5,19 @@ DEAL::face 0
 DEAL::inside has value
 DEAL::face 1
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::face 1
+DEAL::outside has value
 DEAL::Advanced constructor:
 DEAL::face 0
 DEAL::inside has value
 DEAL::face 1
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::face 1
+DEAL::outside has value
 DEAL::
 DEAL::dim = 2
 DEAL::Simple constructor:
@@ -23,6 +31,16 @@ DEAL::face 2
 DEAL::inside has value
 DEAL::face 3
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::face 3
+DEAL::outside has value
 DEAL::Advanced constructor:
 DEAL::face 0
 DEAL::inside has value
@@ -34,6 +52,16 @@ DEAL::face 2
 DEAL::inside has value
 DEAL::face 3
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::face 3
+DEAL::outside has value
 DEAL::
 DEAL::dim = 3
 DEAL::Simple constructor:
@@ -53,6 +81,22 @@ DEAL::face 4
 DEAL::inside has value
 DEAL::face 5
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 3
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 4
+DEAL::inside has value
+DEAL::face 5
+DEAL::outside has value
 DEAL::Advanced constructor:
 DEAL::face 0
 DEAL::inside has value
@@ -70,4 +114,20 @@ DEAL::face 4
 DEAL::inside has value
 DEAL::face 5
 DEAL::outside has value
+DEAL::face 0
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 1
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 2
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 3
+DEAL::inside has value
+DEAL::outside has value
+DEAL::face 4
+DEAL::inside has value
+DEAL::face 5
+DEAL::outside has value
 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.