Merge pull request #13834 from peterrum/multiply_dof_numbers

Generalize FETools::Compositing::multiply_dof_numbers()

include/deal.II/fe/fe_base.h

@@ -980,9 +980,9 @@ internal::GenericDoFsPerObject::generate(const FiniteElementData<dim> &fe)
 
   internal::GenericDoFsPerObject result;
 
-  result.dofs_per_object_exclusive.resize(dim + 1);
-  result.dofs_per_object_inclusive.resize(dim + 1);
-  result.object_index.resize(dim + 1);
+  result.dofs_per_object_exclusive.resize(4);
+  result.dofs_per_object_inclusive.resize(4);
+  result.object_index.resize(4);
 
   unsigned int counter = 0;
 
@@ -1023,24 +1023,32 @@ internal::GenericDoFsPerObject::generate(const FiniteElementData<dim> &fe)
       }
 
   {
-    result.dofs_per_object_exclusive[dim].emplace_back(
-      fe.template n_dofs_per_object<dim>());
+    const auto c = fe.template n_dofs_per_object<dim>();
+
+    result.dofs_per_object_exclusive[dim].emplace_back(c);
     result.dofs_per_object_inclusive[dim].emplace_back(fe.n_dofs_per_cell());
     result.object_index[dim].emplace_back(counter);
+
+    counter += c;
   }
 
-  result.first_object_index_on_face.resize(dim);
+  for (unsigned int d = dim + 1; d <= 3; ++d)
+    {
+      result.dofs_per_object_exclusive[d].emplace_back(0);
+      result.dofs_per_object_inclusive[d].emplace_back(0);
+      result.object_index[d].emplace_back(counter);
+    }
+
+  result.first_object_index_on_face.resize(3);
   for (unsigned int face_no : reference_cell.face_indices())
     {
       result.first_object_index_on_face[0].emplace_back(0);
 
-      if (dim >= 2)
-        result.first_object_index_on_face[1].emplace_back(
-          fe.get_first_face_line_index(face_no));
+      result.first_object_index_on_face[1].emplace_back(
+        fe.get_first_face_line_index(face_no));
 
-      if (dim == 3)
-        result.first_object_index_on_face[2].emplace_back(
-          fe.get_first_face_quad_index(face_no));
+      result.first_object_index_on_face[2].emplace_back(
+        fe.get_first_face_quad_index(face_no));
     }
 
   return result;

include/deal.II/fe/fe_tools.templates.h

@@ -85,11 +85,6 @@ namespace FETools
     {
       AssertDimension(fes.size(), multiplicities.size());
 
-      unsigned int multiplied_dofs_per_vertex = 0;
-      unsigned int multiplied_dofs_per_line   = 0;
-      unsigned int multiplied_dofs_per_quad   = 0;
-      unsigned int multiplied_dofs_per_hex    = 0;
-
       unsigned int multiplied_n_components = 0;
 
       unsigned int degree = 0; // degree is the maximal degree of the components
@@ -103,23 +98,69 @@ namespace FETools
             break;
           }
 
+      dealii::internal::GenericDoFsPerObject dpo;
+
+      std::vector<dealii::internal::GenericDoFsPerObject> dpos_in(fes.size());
+
       for (unsigned int i = 0; i < fes.size(); ++i)
         if (multiplicities[i] > 0)
-          {
-            // TODO: the implementation makes the assumption that all faces have
-            // the same number of dofs -> don't construct DPO but
-            // PrecomputedFiniteElementData
-            AssertDimension(fes[i]->n_unique_quads(), 1);
-
-            multiplied_dofs_per_vertex +=
-              fes[i]->n_dofs_per_vertex() * multiplicities[i];
-            multiplied_dofs_per_line +=
-              fes[i]->n_dofs_per_line() * multiplicities[i];
-            multiplied_dofs_per_quad +=
-              fes[i]->n_dofs_per_quad(0) * multiplicities[i];
-            multiplied_dofs_per_hex +=
-              fes[i]->n_dofs_per_hex() * multiplicities[i];
+          dpos_in[i] =
+            dealii::internal::GenericDoFsPerObject::generate(*fes[i]);
+
+      // helper function to fill a vector of GenericDoFsPerObject
+      // according to multiplicities
+      const auto fill_dpo_vector =
+        [&](const std::function<std::vector<std::vector<unsigned int>> &(
+              dealii::internal::GenericDoFsPerObject &)> &get_vector) {
+          auto &vector_dst = get_vector(dpo);
+
+          // allocate memory
+          for (unsigned int i = 0; i < fes.size(); ++i)
+            if (multiplicities[i] > 0)
+              {
+                const auto &vector_src = get_vector(dpos_in[i]);
+
+                vector_dst.resize(vector_src.size());
+
+                for (unsigned int j = 0; j < vector_src.size(); ++j)
+                  vector_dst[j].assign(vector_src[j].size(), 0);
+
+                break;
+              }
+
+          // fill vector according to multiplicities
+          for (unsigned int i = 0; i < fes.size(); ++i)
+            if (multiplicities[i] > 0)
+              {
+                const auto &vector_src = get_vector(dpos_in[i]);
 
+                for (unsigned int j = 0; j < vector_src.size(); ++j)
+                  for (unsigned int k = 0; k < vector_src[j].size(); ++k)
+                    vector_dst[j][k] += vector_src[j][k] * multiplicities[i];
+              }
+        };
+
+      // go through each field of GenericDoFsPerObject
+      fill_dpo_vector(
+        [](auto &dpo) -> std::vector<std::vector<unsigned int>> & {
+          return dpo.dofs_per_object_exclusive;
+        });
+      fill_dpo_vector(
+        [](auto &dpo) -> std::vector<std::vector<unsigned int>> & {
+          return dpo.dofs_per_object_inclusive;
+        });
+      fill_dpo_vector(
+        [](auto &dpo) -> std::vector<std::vector<unsigned int>> & {
+          return dpo.object_index;
+        });
+      fill_dpo_vector(
+        [](auto &dpo) -> std::vector<std::vector<unsigned int>> & {
+          return dpo.first_object_index_on_face;
+        });
+
+      for (unsigned int i = 0; i < fes.size(); ++i)
+        if (multiplicities[i] > 0)
+          {
             multiplied_n_components +=
               fes[i]->n_components() * multiplicities[i];
 
@@ -150,14 +191,6 @@ namespace FETools
               total_conformity & fes[index]->conforming_space);
       }
 
-      std::vector<unsigned int> dpo;
-      dpo.push_back(multiplied_dofs_per_vertex);
-      dpo.push_back(multiplied_dofs_per_line);
-      if (dim > 1)
-        dpo.push_back(multiplied_dofs_per_quad);
-      if (dim > 2)
-        dpo.push_back(multiplied_dofs_per_hex);
-
       BlockIndices block_indices(0, 0);
 
       for (unsigned int base = 0; base < fes.size(); ++base)

source/fe/fe_data.cc

@@ -94,6 +94,19 @@ namespace internal
 
     return result;
   }
+
+  unsigned int
+  number_unique_entries(const std::vector<unsigned int> &vector)
+  {
+    if (std::all_of(vector.begin(), vector.end(), [&](const auto &e) {
+          return e == vector.front();
+        }))
+      {
+        return 1;
+      }
+    else
+      return vector.size();
+  }
 } // namespace internal
 
 
@@ -133,6 +146,8 @@ FiniteElementData<dim>::FiniteElementData(
                       block_indices)
 {}
 
+
+
 template <int dim>
 FiniteElementData<dim>::FiniteElementData(
   const internal::GenericDoFsPerObject &data,
@@ -142,16 +157,17 @@ FiniteElementData<dim>::FiniteElementData(
   const Conformity                      conformity,
   const BlockIndices &                  block_indices)
   : reference_cell_kind(reference_cell)
-  , number_unique_quads(data.dofs_per_object_inclusive[2].size())
-  , number_unique_faces(data.dofs_per_object_inclusive[dim - 1].size())
+  , number_unique_quads(
+      internal::number_unique_entries(data.dofs_per_object_inclusive[2]))
+  , number_unique_faces(
+      internal::number_unique_entries(data.dofs_per_object_inclusive[dim - 1]))
   , dofs_per_vertex(data.dofs_per_object_exclusive[0][0])
   , dofs_per_line(data.dofs_per_object_exclusive[1][0])
-  , n_dofs_on_quad(dim > 1 ? data.dofs_per_object_exclusive[2] :
-                             std::vector<unsigned int>{0})
+  , n_dofs_on_quad(data.dofs_per_object_exclusive[2])
   , dofs_per_quad(n_dofs_on_quad[0])
   , dofs_per_quad_max(
       *max_element(n_dofs_on_quad.begin(), n_dofs_on_quad.end()))
-  , dofs_per_hex(dim > 2 ? data.dofs_per_object_exclusive[3][0] : 0)
+  , dofs_per_hex(data.dofs_per_object_exclusive[3][0])
   , first_line_index(data.object_index[1][0])
   , first_index_of_quads(data.object_index[2])
   , first_quad_index(first_index_of_quads[0])

tests/fe/generic_dofs_per_object_01.output

@@ -1,27 +1,27 @@
 
 DEAL::FE_Q<1>(1)
-DEAL::1 1 0 
-DEAL::1 1 2 
-DEAL::0 1 2 
-DEAL::0 0 
+DEAL::1 1 0 0 0 
+DEAL::1 1 2 0 0 
+DEAL::0 1 2 2 2 
+DEAL::0 0 1 1 2 2 
 DEAL::
 DEAL::FE_Q<1>(2)
-DEAL::1 1 1 
-DEAL::1 1 3 
-DEAL::0 1 2 
-DEAL::0 0 
+DEAL::1 1 1 0 0 
+DEAL::1 1 3 0 0 
+DEAL::0 1 2 3 3 
+DEAL::0 0 1 1 2 2 
 DEAL::
 DEAL::FE_Q<2>(1)
-DEAL::1 1 1 1 0 0 0 0 0 
-DEAL::1 1 1 1 2 2 2 2 4 
-DEAL::0 1 2 3 4 4 4 4 4 
-DEAL::0 0 0 0 2 2 2 2 
+DEAL::1 1 1 1 0 0 0 0 0 0 
+DEAL::1 1 1 1 2 2 2 2 4 0 
+DEAL::0 1 2 3 4 4 4 4 4 4 
+DEAL::0 0 0 0 2 2 2 2 4 4 4 4 
 DEAL::
 DEAL::FE_Q<2>(2)
-DEAL::1 1 1 1 1 1 1 1 1 
-DEAL::1 1 1 1 3 3 3 3 9 
-DEAL::0 1 2 3 4 5 6 7 8 
-DEAL::0 0 0 0 2 2 2 2 
+DEAL::1 1 1 1 1 1 1 1 1 0 
+DEAL::1 1 1 1 3 3 3 3 9 0 
+DEAL::0 1 2 3 4 5 6 7 8 9 
+DEAL::0 0 0 0 2 2 2 2 5 5 5 5 
 DEAL::
 DEAL::FE_Q<3>(1)
 DEAL::1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
@@ -36,16 +36,16 @@ DEAL::0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
 DEAL::0 0 0 0 0 0 4 4 4 4 4 4 8 8 8 8 8 8 
 DEAL::
 DEAL::FE_SimplexP<2>(1)
-DEAL::1 1 1 0 0 0 0 
-DEAL::1 1 1 2 2 2 3 
-DEAL::0 1 2 3 3 3 3 
-DEAL::0 0 0 2 2 2 
+DEAL::1 1 1 0 0 0 0 0 
+DEAL::1 1 1 2 2 2 3 0 
+DEAL::0 1 2 3 3 3 3 3 
+DEAL::0 0 0 2 2 2 3 3 3 
 DEAL::
 DEAL::FE_SimplexP<2>(2)
-DEAL::1 1 1 1 1 1 0 
-DEAL::1 1 1 3 3 3 6 
-DEAL::0 1 2 3 4 5 6 
-DEAL::0 0 0 2 2 2 
+DEAL::1 1 1 1 1 1 0 0 
+DEAL::1 1 1 3 3 3 6 0 
+DEAL::0 1 2 3 4 5 6 6 
+DEAL::0 0 0 2 2 2 4 4 4 
 DEAL::
 DEAL::FE_SimplexP<3>(1)
 DEAL::1 1 1 1 0 0 0 0 0 0 0 0 0 0 0