]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Pass face_no to FE functions 10783/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 1 Aug 2020 12:43:56 +0000 (14:43 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 13 Aug 2020 17:51:29 +0000 (19:51 +0200)
15 files changed:
include/deal.II/dofs/dof_accessor.templates.h
include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/multigrid/mg_constrained_dofs.h
include/deal.II/numerics/matrix_creator.templates.h
include/deal.II/numerics/vector_tools_boundary.templates.h
include/deal.II/numerics/vector_tools_constraints.templates.h
include/deal.II/numerics/vector_tools_project.templates.h
source/dofs/dof_handler.cc
source/dofs/dof_handler_policy.cc
source/dofs/dof_tools.cc
source/dofs/dof_tools_constraints.cc
source/dofs/dof_tools_sparsity.cc
source/grid/grid_tools_dof_handlers.cc
source/multigrid/mg_tools.cc

index c7c9cc567f480ce54e67be27ab716f91a6b1f34c..e84ffdbda94264adb1646c963d5b032d962d7585 100644 (file)
@@ -351,7 +351,7 @@ namespace internal
               dof_handler.object_dof_ptr[obj_level][d][obj_index];
             const unsigned int ptr_1 =
               ptr_0 +
-              dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>();
+              dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(0);
 
             return {ptr_0, ptr_1};
           }
@@ -819,11 +819,11 @@ namespace internal
           {
             const auto &fe = accessor.get_fe(fe_index_);
 
-            const unsigned int                          //
-              dofs_per_vertex = fe.n_dofs_per_vertex(), //
-              dofs_per_line   = fe.n_dofs_per_line(),   //
-              dofs_per_quad   = fe.n_dofs_per_quad(),   //
-              dofs_per_hex    = fe.n_dofs_per_hex();    //
+            const unsigned int                                   //
+              dofs_per_vertex = fe.n_dofs_per_vertex(),          //
+              dofs_per_line   = fe.n_dofs_per_line(),            //
+              dofs_per_quad   = fe.n_dofs_per_quad(0 /*dummy*/), //
+              dofs_per_hex    = fe.n_dofs_per_hex();             //
 
             const unsigned int inner_dofs =
               structdim == 1 ? dofs_per_line :
@@ -1211,7 +1211,7 @@ namespace internal
           (void)fe_index;
 
           for (unsigned int d = 0;
-               d < fe.template n_dofs_per_object<structdim_>();
+               d < fe.template n_dofs_per_object<structdim_>(0);
                ++d, ++index)
             index_value[index] = accessor.mg_dof_index(level, mapping(d));
         }
@@ -1291,7 +1291,7 @@ namespace internal
           (void)fe_index;
 
           for (unsigned int d = 0;
-               d < fe.template n_dofs_per_object<structdim_>();
+               d < fe.template n_dofs_per_object<structdim_>(0);
                ++d, ++index)
             accessor.set_mg_dof_index(level, mapping(d), index_value[index]);
         }
index 7e993a6d11f328e71abf5aa6f951a95adb2e50d6..ab14ffdbb2009df4fa2365ff901057cb5d6f4096 100644 (file)
@@ -140,25 +140,31 @@ namespace CUDAWrappers
       void
       setup_constraint_weigths(unsigned int fe_degree)
       {
+        const unsigned int face_no =
+          0; // we assume that all faces have the same number of dofs
+
         FE_Q<2>            fe_q(fe_degree);
-        FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(),
-                                                fe_q.n_dofs_per_face());
-        fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix);
+        FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(face_no),
+                                                fe_q.n_dofs_per_face(face_no));
+        fe_q.get_subface_interpolation_matrix(fe_q,
+                                              0,
+                                              interpolation_matrix,
+                                              face_no);
 
         std::vector<unsigned int> mapping =
           FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree);
 
-        FullMatrix<double> mapped_matrix(fe_q.n_dofs_per_face(),
-                                         fe_q.n_dofs_per_face());
-        for (unsigned int i = 0; i < fe_q.n_dofs_per_face(); ++i)
-          for (unsigned int j = 0; j < fe_q.n_dofs_per_face(); ++j)
+        FullMatrix<double> mapped_matrix(fe_q.n_dofs_per_face(face_no),
+                                         fe_q.n_dofs_per_face(face_no));
+        for (unsigned int i = 0; i < fe_q.n_dofs_per_face(face_no); ++i)
+          for (unsigned int j = 0; j < fe_q.n_dofs_per_face(face_no); ++j)
             mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]);
 
         cudaError_t error_code =
           cudaMemcpyToSymbol(internal::constraint_weights,
                              &mapped_matrix[0][0],
-                             sizeof(double) * fe_q.n_dofs_per_face() *
-                               fe_q.n_dofs_per_face());
+                             sizeof(double) * fe_q.n_dofs_per_face(face_no) *
+                               fe_q.n_dofs_per_face(face_no));
         AssertCuda(error_code);
       }
     } // namespace internal
index 19e6710dac844e77d3d463e1cd446174fcaa3ea8..8519cb123cd6a69cfedb0f650e9af6d46fbf0341 100644 (file)
@@ -760,7 +760,8 @@ namespace internal
             const FiniteElement<dim> &fe = fes[fe_index];
             // cache number of finite elements and dofs_per_cell
             dof_info[no].dofs_per_cell.push_back(fe.n_dofs_per_cell());
-            dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face());
+            dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face(
+              0)); // we assume that all faces have the same number of dofs
             dof_info[no].dimension       = dim;
             dof_info[no].n_base_elements = fe.n_base_elements();
             dof_info[no].n_components.resize(dof_info[no].n_base_elements);
index 7baf0b58b6bec5a90a9070399412c3325299d712..53e9d88dc758b3410ea453d0af0d7d2c518bef07 100644 (file)
@@ -273,7 +273,7 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim, spacedim> &dof)
                     }
 
                   const unsigned int dofs_per_face =
-                    cell->face(f)->get_fe(0).n_dofs_per_face();
+                    dof.get_fe(0).n_dofs_per_face(f);
                   std::vector<types::global_dof_index> dofs_1(dofs_per_face);
                   std::vector<types::global_dof_index> dofs_2(dofs_per_face);
 
index 3b37e5b399cac4fba131f0cd389716d6caaf95da..57f067889ed39c639623d674a29eae6e022ced37 100644 (file)
@@ -948,8 +948,6 @@ namespace MatrixCreator
       const bool fe_is_system    = (n_components != 1);
       const bool fe_is_primitive = fe.is_primitive();
 
-      const unsigned int dofs_per_face = fe.n_dofs_per_face();
-
       copy_data.cell          = cell;
       copy_data.dofs_per_cell = fe.n_dofs_per_cell();
 
@@ -973,7 +971,7 @@ namespace MatrixCreator
       copy_data.dofs.resize(copy_data.dofs_per_cell);
       cell->get_dof_indices(copy_data.dofs);
 
-      std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face);
+      std::vector<types::global_dof_index> dofs_on_face_vector;
 
       // Because CopyData objects are reused and emplace_back is
       // used, dof_is_on_face, cell_matrix, and cell_vector must be
@@ -1147,6 +1145,7 @@ namespace MatrixCreator
               }
 
 
+            dofs_on_face_vector.resize(fe.n_dofs_per_face(face));
             cell->face(face)->get_dof_indices(dofs_on_face_vector);
             // for each dof on the cell, have a flag whether it is on
             // the face
@@ -1379,7 +1378,6 @@ namespace MatrixCreator
       const FiniteElement<dim, spacedim> &fe              = cell->get_fe();
       const bool                          fe_is_system    = (n_components != 1);
       const bool                          fe_is_primitive = fe.is_primitive();
-      const unsigned int                  dofs_per_face = fe.n_dofs_per_face();
 
       copy_data.cell          = cell;
       copy_data.dofs_per_cell = fe.n_dofs_per_cell();
@@ -1407,7 +1405,7 @@ namespace MatrixCreator
       std::vector<number>         rhs_values_scalar;
       std::vector<Vector<number>> rhs_values_system;
 
-      std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face);
+      std::vector<types::global_dof_index> dofs_on_face_vector;
 
       copy_data.dofs.resize(copy_data.dofs_per_cell);
       cell->get_dof_indices(copy_data.dofs);
@@ -1599,6 +1597,7 @@ namespace MatrixCreator
                   }
               }
 
+            dofs_on_face_vector.resize(fe.n_dofs_per_face(face));
             cell->face(face)->get_dof_indices(dofs_on_face_vector,
                                               cell->active_fe_index());
             // for each dof on the cell, have a
index 0e8c79faf022b6fbef80e6661dc49e2bf7c1fa4b..7801f0d42e5be9f8ba9ff43ff2045f632485a3ce 100644 (file)
@@ -131,12 +131,14 @@ namespace VectorTools
                                                    function_values);
 
                   for (unsigned int i = 0; i < fe.n_dofs_per_vertex(); ++i)
-                    if (component_mask[fe.face_system_to_component_index(i)
+                    if (component_mask[fe.face_system_to_component_index(
+                                           i, direction)
                                          .first])
                       boundary_values[cell->vertex_dof_index(
                         direction, i, cell->active_fe_index())] =
                         function_values(
-                          fe.face_system_to_component_index(i).first);
+                          fe.face_system_to_component_index(i, direction)
+                            .first);
                 }
         }
       else // dim > 1
@@ -157,6 +159,9 @@ namespace VectorTools
           dof_values_system.reserve(
             dof.get_fe_collection().max_dofs_per_face());
 
+          // TODO: get support for each face -> PR #10764
+          const unsigned int face_no = 0;
+
           // before we start with the loop over all cells create an hp::FEValues
           // object that holds the interpolation points of all finite elements
           // that may ever be in use
@@ -173,9 +178,9 @@ namespace VectorTools
               //
               // to do this, we check whether the FE has support points on the
               // face at all:
-              if (fe.has_face_support_points())
-                q_collection.push_back(
-                  Quadrature<dim - 1>(fe.get_unit_face_support_points()));
+              if (fe.has_face_support_points(face_no))
+                q_collection.push_back(Quadrature<dim - 1>(
+                  fe.get_unit_face_support_points(face_no)));
               else
                 {
                   // if not, then we should try a more clever way. the idea is
@@ -194,13 +199,15 @@ namespace VectorTools
                   // components we are interested in are primitive (by the above
                   // check), we can safely put such a check in front
                   std::vector<Point<dim - 1>> unit_support_points(
-                    fe.n_dofs_per_face());
+                    fe.n_dofs_per_face(face_no));
 
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
-                    if (fe.is_primitive(fe.face_to_cell_index(i, 0)))
-                      if (component_mask[fe.face_system_to_component_index(i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i)
+                    if (fe.is_primitive(fe.face_to_cell_index(i, face_no)))
+                      if (component_mask[fe.face_system_to_component_index(
+                                             i, face_no)
                                            .first] == true)
-                        unit_support_points[i] = fe.unit_face_support_point(i);
+                        unit_support_points[i] =
+                          fe.unit_face_support_point(i, face_no);
 
                   q_collection.push_back(
                     Quadrature<dim - 1>(unit_support_points));
@@ -260,7 +267,7 @@ namespace VectorTools
                   // element in use here has DoFs on the face at all
                   if ((function_map.find(boundary_component) !=
                        function_map.end()) &&
-                      (cell->get_fe().n_dofs_per_face() > 0))
+                      (cell->get_fe().n_dofs_per_face(face_no) > 0))
                     {
                       // face is of the right component
                       x_fe_values.reinit(cell, face_no);
@@ -269,7 +276,7 @@ namespace VectorTools
 
                       // get indices, physical location and boundary values of
                       // dofs on this face
-                      face_dofs.resize(fe.n_dofs_per_face());
+                      face_dofs.resize(fe.n_dofs_per_face(face_no));
                       face->get_dof_indices(face_dofs, cell->active_fe_index());
                       const std::vector<Point<spacedim>> &dof_locations =
                         fe_values.get_quadrature_points();
@@ -278,12 +285,14 @@ namespace VectorTools
                         {
                           // resize array. avoid construction of a memory
                           // allocating temporary if possible
-                          if (dof_values_system.size() < fe.n_dofs_per_face())
-                            dof_values_system.resize(fe.n_dofs_per_face(),
-                                                     Vector<number>(
-                                                       fe.n_components()));
+                          if (dof_values_system.size() <
+                              fe.n_dofs_per_face(face_no))
+                            dof_values_system.resize(
+                              fe.n_dofs_per_face(face_no),
+                              Vector<number>(fe.n_components()));
                           else
-                            dof_values_system.resize(fe.n_dofs_per_face());
+                            dof_values_system.resize(
+                              fe.n_dofs_per_face(face_no));
 
                           function_map.find(boundary_component)
                             ->second->vector_value_list(dof_locations,
@@ -298,7 +307,8 @@ namespace VectorTools
                               unsigned int component;
                               if (fe.is_primitive())
                                 component =
-                                  fe.face_system_to_component_index(i).first;
+                                  fe.face_system_to_component_index(i, face_no)
+                                    .first;
                               else
                                 {
                                   // non-primitive case. make sure that this
@@ -360,7 +370,7 @@ namespace VectorTools
                         // fe has only one component, so save some computations
                         {
                           // get only the one component that this function has
-                          dof_values_scalar.resize(fe.n_dofs_per_face());
+                          dof_values_scalar.resize(fe.n_dofs_per_face(face_no));
                           function_map.find(boundary_component)
                             ->second->value_list(dof_locations,
                                                  dof_values_scalar,
@@ -1026,7 +1036,7 @@ namespace VectorTools
 
           // Compute the degrees of
           // freedom.
-          for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
             if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
                  (fe.system_to_base_index(fe.face_to_cell_index(i, face))
                     .first == base_indices) &&
@@ -1201,7 +1211,7 @@ namespace VectorTools
 
                   // Compute the degrees
                   // of freedom.
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                     if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
                           nullptr) &&
                          (fe.system_to_base_index(
@@ -1271,7 +1281,7 @@ namespace VectorTools
                   for (unsigned int d = 0; d < dim; ++d)
                     tmp[d] = values[q_point](first_vector_component + d);
 
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                     if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
                           nullptr) &&
                          (fe.system_to_base_index(
@@ -1337,7 +1347,7 @@ namespace VectorTools
 
                   unsigned int index = 0;
 
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                     if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
                           nullptr) &&
                          (fe.system_to_base_index(
@@ -1391,7 +1401,7 @@ namespace VectorTools
               {
                 unsigned int index = 0;
 
-                for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                   if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
                        (fe.system_to_base_index(fe.face_to_cell_index(i, face))
                           .first == base_indices) &&
@@ -1431,7 +1441,7 @@ namespace VectorTools
                   for (unsigned int d = 0; d < dim; ++d)
                     tmp[d] = values[q_point](first_vector_component + d);
 
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                     if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
                           nullptr) &&
                          (fe.system_to_base_index(
@@ -1487,7 +1497,7 @@ namespace VectorTools
 
                   unsigned int index = 0;
 
-                  for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                     if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
                           nullptr) &&
                          (fe.system_to_base_index(
@@ -1527,7 +1537,7 @@ namespace VectorTools
 
               unsigned int index = 0;
 
-              for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+              for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                 if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
                      (fe.system_to_base_index(fe.face_to_cell_index(i, face))
                         .first == base_indices) &&
@@ -1583,7 +1593,12 @@ namespace VectorTools
     // corresponding degrees of freedom.
     const unsigned int    superdegree = dof_handler.get_fe().degree;
     const QGauss<dim - 1> reference_face_quadrature(2 * superdegree);
-    const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face();
+
+    // TODO: the implementation makes the assumption that all faces have the
+    // same number of dofs
+    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+    const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face(0);
+
     const hp::FECollection<dim> &fe_collection(dof_handler.get_fe_collection());
     const hp::MappingCollection<dim> mapping_collection(mapping);
     hp::QCollection<dim>             face_quadrature_collection;
@@ -1860,7 +1875,7 @@ namespace VectorTools
                         }
 
                       const unsigned int dofs_per_face =
-                        cell->get_fe().n_dofs_per_face();
+                        cell->get_fe().n_dofs_per_face(face);
 
                       dofs_processed.resize(dofs_per_face);
                       dof_values.resize(dofs_per_face);
@@ -1953,7 +1968,7 @@ namespace VectorTools
                       const unsigned int superdegree = cell->get_fe().degree;
                       const unsigned int degree      = superdegree - 1;
                       const unsigned int dofs_per_face =
-                        cell->get_fe().n_dofs_per_face();
+                        cell->get_fe().n_dofs_per_face(face);
 
                       dofs_processed.resize(dofs_per_face);
                       dof_values.resize(dofs_per_face);
@@ -2428,7 +2443,8 @@ namespace VectorTools
                                                                         1);
 
               unsigned int associated_edge_dof_index = 0;
-              for (unsigned int face_idx = 0; face_idx < fe.n_dofs_per_face();
+              for (unsigned int face_idx = 0;
+                   face_idx < fe.n_dofs_per_face(face);
                    ++face_idx)
                 {
                   const unsigned int cell_idx =
@@ -2663,7 +2679,7 @@ namespace VectorTools
               // Loop over these quad-interior dofs.
               unsigned int associated_face_dof_index = 0;
               for (unsigned int quad_dof_idx = 0;
-                   quad_dof_idx < fe.n_dofs_per_quad();
+                   quad_dof_idx < fe.n_dofs_per_quad(face);
                    ++quad_dof_idx)
                 {
                   const unsigned int face_idx =
@@ -2925,7 +2941,7 @@ namespace VectorTools
                                 }
 
                               const unsigned int dofs_per_face =
-                                cell->get_fe().n_dofs_per_face();
+                                cell->get_fe().n_dofs_per_face(face);
 
                               dofs_processed.resize(dofs_per_face);
                               dof_values.resize(dofs_per_face);
@@ -3051,7 +3067,7 @@ namespace VectorTools
                                 cell->get_fe().degree;
                               const unsigned int degree = superdegree - 1;
                               const unsigned int dofs_per_face =
-                                cell->get_fe().n_dofs_per_face();
+                                cell->get_fe().n_dofs_per_face(face);
 
                               dofs_processed.resize(dofs_per_face);
                               dof_values.resize(dofs_per_face);
@@ -3206,7 +3222,7 @@ namespace VectorTools
                                                                       0};
       std::vector<Vector<double>> values(fe_values.n_quadrature_points,
                                          Vector<double>(2));
-      Vector<double>              dof_values(fe.n_dofs_per_face());
+      Vector<double>              dof_values(fe.n_dofs_per_face(face));
 
       // Get the values of the boundary function at the quadrature points.
       {
@@ -3231,7 +3247,7 @@ namespace VectorTools
                       jacobians[q_point][1][face_coordinate_direction[face]] *
                         jacobians[q_point][1][face_coordinate_direction[face]]);
 
-          for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
             dof_values(i) +=
               tmp * (normals[q_point] *
                      fe_values[vec].value(
@@ -3244,14 +3260,14 @@ namespace VectorTools
         }
 
       std::vector<types::global_dof_index> face_dof_indices(
-        fe.n_dofs_per_face());
+        fe.n_dofs_per_face(face));
 
       cell->face(face)->get_dof_indices(face_dof_indices,
                                         cell->active_fe_index());
 
       // Copy the computed values in the AffineConstraints only, if the degree
       // of freedom is not already constrained.
-      for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
         if (!(constraints.is_constrained(face_dof_indices[i])) &&
             fe.get_nonzero_components(fe.face_to_cell_index(
               i,
@@ -3307,7 +3323,7 @@ namespace VectorTools
           {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
       std::vector<Vector<double>> values(fe_values.n_quadrature_points,
                                          Vector<double>(3));
-      Vector<double>              dof_values_local(fe.n_dofs_per_face());
+      Vector<double>              dof_values_local(fe.n_dofs_per_face(face));
 
       {
         const std::vector<Point<3>> &quadrature_points =
@@ -3340,7 +3356,7 @@ namespace VectorTools
                jacobians[q_point][2][face_coordinate_directions[face][1]] *
                  jacobians[q_point][2][face_coordinate_directions[face][1]]));
 
-          for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
             dof_values_local(i) +=
               tmp * (normals[q_point] *
                      fe_values[vec].value(
@@ -3353,12 +3369,12 @@ namespace VectorTools
         }
 
       std::vector<types::global_dof_index> face_dof_indices(
-        fe.n_dofs_per_face());
+        fe.n_dofs_per_face(face));
 
       cell->face(face)->get_dof_indices(face_dof_indices,
                                         cell->active_fe_index());
 
-      for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
         if (projected_dofs[face_dof_indices[i]] < fe.degree &&
             fe.get_nonzero_components(fe.face_to_cell_index(
               i,
index 410f034c416e0f72c48b0a65f87c27ba9ee58cd5..24477aac63080c4a2a457a178638bd79e6aff38f 100644 (file)
@@ -490,6 +490,11 @@ namespace VectorTools
     for (unsigned int i = 0; i < fe_collection.size(); ++i)
       mapping_collection.push_back(mapping);
 
+    // TODO: the implementation makes the assumption that all faces have the
+    // same number of dofs
+    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+    const unsigned int face_no = 0;
+
     // now also create a quadrature collection for the faces of a cell. fill
     // it with a quadrature formula with the support points on faces for each
     // FE
@@ -497,9 +502,10 @@ namespace VectorTools
     for (unsigned int i = 0; i < fe_collection.size(); ++i)
       {
         const std::vector<Point<dim - 1>> &unit_support_points =
-          fe_collection[i].get_unit_face_support_points();
+          fe_collection[i].get_unit_face_support_points(face_no);
 
-        Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(),
+        Assert(unit_support_points.size() ==
+                 fe_collection[i].n_dofs_per_face(face_no),
                ExcInternalError());
 
         face_quadrature_collection.push_back(
@@ -549,7 +555,7 @@ namespace VectorTools
                 cell->face(face_no);
 
               // get the indices of the dofs on this cell...
-              face_dofs.resize(fe.n_dofs_per_face());
+              face_dofs.resize(fe.n_dofs_per_face(face_no));
               face->get_dof_indices(face_dofs, cell->active_fe_index());
 
               x_fe_face_values.reinit(cell, face_no);
@@ -559,7 +565,7 @@ namespace VectorTools
               // then identify which of them correspond to the selected set of
               // vector components
               for (unsigned int i = 0; i < face_dofs.size(); ++i)
-                if (fe.face_system_to_component_index(i).first ==
+                if (fe.face_system_to_component_index(i, face_no).first ==
                     first_vector_component)
                   {
                     // find corresponding other components of vector
@@ -572,18 +578,19 @@ namespace VectorTools
                         "Error: the finite element does not have enough components "
                         "to define a normal direction."));
 
-                    for (unsigned int k = 0; k < fe.n_dofs_per_face(); ++k)
+                    for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
+                         ++k)
                       if ((k != i) &&
                           (face_quadrature_collection[cell->active_fe_index()]
                              .point(k) ==
                            face_quadrature_collection[cell->active_fe_index()]
                              .point(i)) &&
-                          (fe.face_system_to_component_index(k).first >=
-                           first_vector_component) &&
-                          (fe.face_system_to_component_index(k).first <
+                          (fe.face_system_to_component_index(k, face_no)
+                             .first >= first_vector_component) &&
+                          (fe.face_system_to_component_index(k, face_no).first <
                            first_vector_component + dim))
                         vector_dofs.dof_indices
-                          [fe.face_system_to_component_index(k).first -
+                          [fe.face_system_to_component_index(k, face_no).first -
                            first_vector_component] = face_dofs[k];
 
                     for (unsigned int d = 0; d < dim; ++d)
@@ -1076,6 +1083,11 @@ namespace VectorTools
     for (unsigned int i = 0; i < fe_collection.size(); ++i)
       mapping_collection.push_back(mapping);
 
+    // TODO: the implementation makes the assumption that all faces have the
+    // same number of dofs
+    AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+    const unsigned int face_no = 0;
+
     // now also create a quadrature collection for the faces of a cell. fill
     // it with a quadrature formula with the support points on faces for each
     // FE
@@ -1083,9 +1095,10 @@ namespace VectorTools
     for (unsigned int i = 0; i < fe_collection.size(); ++i)
       {
         const std::vector<Point<dim - 1>> &unit_support_points =
-          fe_collection[i].get_unit_face_support_points();
+          fe_collection[i].get_unit_face_support_points(face_no);
 
-        Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(),
+        Assert(unit_support_points.size() ==
+                 fe_collection[i].n_dofs_per_face(face_no),
                ExcInternalError());
 
         face_quadrature_collection.push_back(
@@ -1123,7 +1136,7 @@ namespace VectorTools
                 cell->face(face_no);
 
               // get the indices of the dofs on this cell...
-              face_dofs.resize(fe.n_dofs_per_face());
+              face_dofs.resize(fe.n_dofs_per_face(face_no));
               face->get_dof_indices(face_dofs, cell->active_fe_index());
 
               x_fe_face_values.reinit(cell, face_no);
@@ -1133,22 +1146,23 @@ namespace VectorTools
               std::map<types::global_dof_index, double> dof_to_b_value;
 
               unsigned int n_scalar_indices = 0;
-              cell_vector_dofs.resize(fe.n_dofs_per_face());
-              for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+              cell_vector_dofs.resize(fe.n_dofs_per_face(face_no));
+              for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i)
                 {
-                  if (fe.face_system_to_component_index(i).first >=
+                  if (fe.face_system_to_component_index(i, face_no).first >=
                         first_vector_component &&
-                      fe.face_system_to_component_index(i).first <
+                      fe.face_system_to_component_index(i, face_no).first <
                         first_vector_component + dim)
                     {
                       const unsigned int component =
-                        fe.face_system_to_component_index(i).first -
+                        fe.face_system_to_component_index(i, face_no).first -
                         first_vector_component;
-                      n_scalar_indices =
-                        std::max(n_scalar_indices,
-                                 fe.face_system_to_component_index(i).second +
-                                   1);
-                      cell_vector_dofs[fe.face_system_to_component_index(i)
+                      n_scalar_indices = std::max(
+                        n_scalar_indices,
+                        fe.face_system_to_component_index(i, face_no).second +
+                          1);
+                      cell_vector_dofs[fe.face_system_to_component_index(
+                                           i, face_no)
                                          .second][component] = face_dofs[i];
 
                       const Point<dim> point = fe_values.quadrature_point(i);
index 6db682e4e1d2fcc1ec250cdb56501de536f3ffe4..56c522bcba8a82e721a22e2de07a378cb58e09fe 100644 (file)
@@ -83,10 +83,10 @@ namespace VectorTools
         for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->at_boundary(f))
             {
-              face_dof_indices.resize(cell->get_fe().n_dofs_per_face());
+              face_dof_indices.resize(cell->get_fe().n_dofs_per_face(f));
               cell->face(f)->get_dof_indices(face_dof_indices,
                                              cell->active_fe_index());
-              for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face();
+              for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face(f);
                    ++i)
                 // enter zero boundary values
                 // for all boundary nodes
index ea798e48a11cc5159420a2b5abc9c7a1f7b0a59a..1335bd6f4cb1ef15ff321e74c73eaaaebfa9891d 100644 (file)
@@ -309,7 +309,8 @@ namespace internal
           {
             dof_handler.object_dof_indices[i][2].resize(
               dof_handler.tria->n_raw_cells(i) *
-                dof_handler.get_fe().n_dofs_per_quad(),
+                dof_handler.get_fe().n_dofs_per_quad(
+                  0 /*note: in 2D there is only one quad*/),
               numbers::invalid_dof_index);
 
             dof_handler.object_dof_ptr[i][2].reserve(
@@ -317,7 +318,8 @@ namespace internal
             for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
                  j++)
               dof_handler.object_dof_ptr[i][2].push_back(
-                j * dof_handler.get_fe().n_dofs_per_quad());
+                j * dof_handler.get_fe().n_dofs_per_quad(
+                      0 /*note: in 2D there is only one quad*/));
 
             dof_handler.cell_dof_cache_indices[i].resize(
               dof_handler.tria->n_raw_cells(i) *
@@ -421,7 +423,7 @@ namespace internal
                   {
                     const auto &face = cell->face(face_index);
                     const auto  n_dofs_per_quad =
-                      dof_handler.get_fe().n_dofs_per_quad(/*face_index*/);
+                      dof_handler.get_fe().n_dofs_per_quad(face_index);
 
                     auto &n_dofs_per_quad_target =
                       dof_handler.object_dof_ptr[0][2][face->index() + 1];
@@ -549,9 +551,10 @@ namespace internal
               std::make_unique<
                 internal::DoFHandlerImplementation::DoFLevel<2>>());
             dof_handler.mg_levels.back()->dof_object.dofs =
-              std::vector<types::global_dof_index>(tria.n_raw_quads(i) *
-                                                     fe.n_dofs_per_quad(),
-                                                   numbers::invalid_dof_index);
+              std::vector<types::global_dof_index>(
+                tria.n_raw_quads(i) *
+                  fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
+                numbers::invalid_dof_index);
           }
 
         dof_handler.mg_faces =
@@ -635,10 +638,13 @@ namespace internal
           std::vector<types::global_dof_index>(tria.n_raw_lines() *
                                                  fe.n_dofs_per_line(),
                                                numbers::invalid_dof_index);
-        dof_handler.mg_faces->quads.dofs =
-          std::vector<types::global_dof_index>(tria.n_raw_quads() *
-                                                 fe.n_dofs_per_quad(),
-                                               numbers::invalid_dof_index);
+
+        // TODO: the implementation makes the assumption that all faces have the
+        // same number of dofs
+        AssertDimension(fe.n_unique_faces(), 1);
+        dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
+          tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
+          numbers::invalid_dof_index);
 
         const unsigned int n_vertices = tria.n_vertices();
 
@@ -1331,17 +1337,18 @@ namespace internal
                           fe_slots_needed = 1;
                           n_face_slots +=
                             dof_handler.get_fe(cell->active_fe_index())
-                              .template n_dofs_per_object<dim - 1>();
+                              .template n_dofs_per_object<dim - 1>(face);
                         }
                       else
                         {
                           fe_slots_needed = 2;
                           n_face_slots +=
                             dof_handler.get_fe(cell->active_fe_index())
-                              .template n_dofs_per_object<dim - 1>() +
+                              .template n_dofs_per_object<dim - 1>(face) +
                             dof_handler
                               .get_fe(cell->neighbor(face)->active_fe_index())
-                              .template n_dofs_per_object<dim - 1>();
+                              .template n_dofs_per_object<dim - 1>(
+                                cell->neighbor_face_no(face));
                         }
 
                       // mark this face as visited
@@ -1412,7 +1419,7 @@ namespace internal
                           const unsigned int fe = cell->active_fe_index();
                           const unsigned int n_dofs =
                             dof_handler.get_fe(fe)
-                              .template n_dofs_per_object<dim - 1>();
+                              .template n_dofs_per_object<dim - 1>(face);
                           const unsigned int offset =
                             dof_handler
                               .hp_object_fe_ptr[d][cell->face(face)->index()];
@@ -1426,20 +1433,25 @@ namespace internal
                         }
                       else
                         {
-                          unsigned int fe_1 = cell->active_fe_index();
+                          unsigned int fe_1      = cell->active_fe_index();
+                          unsigned int face_no_1 = face;
                           unsigned int fe_2 =
                             cell->neighbor(face)->active_fe_index();
+                          unsigned int face_no_2 = cell->neighbor_face_no(face);
 
                           if (fe_2 < fe_1)
-                            std::swap(fe_1, fe_2);
+                            {
+                              std::swap(fe_1, fe_2);
+                              std::swap(face_no_1, face_no_2);
+                            }
 
                           const unsigned int n_dofs_1 =
                             dof_handler.get_fe(fe_1)
-                              .template n_dofs_per_object<dim - 1>();
+                              .template n_dofs_per_object<dim - 1>(face_no_1);
 
                           const unsigned int n_dofs_2 =
                             dof_handler.get_fe(fe_2)
-                              .template n_dofs_per_object<dim - 1>();
+                              .template n_dofs_per_object<dim - 1>(face_no_2);
 
                           const unsigned int offset =
                             dof_handler
@@ -2316,7 +2328,7 @@ DoFHandler<dim, spacedim>::n_boundary_dofs() const
             if (face->at_boundary())
               {
                 const unsigned int dofs_per_face =
-                  cell->get_fe().n_dofs_per_face();
+                  cell->get_fe().n_dofs_per_face(iface);
                 dofs_on_face.resize(dofs_per_face);
 
                 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
@@ -2368,7 +2380,7 @@ DoFHandler<dim, spacedim>::n_boundary_dofs(
                 (boundary_ids.find(boundary_id) != boundary_ids.end()))
               {
                 const unsigned int dofs_per_face =
-                  cell->get_fe().n_dofs_per_face();
+                  cell->get_fe().n_dofs_per_face(iface);
                 dofs_on_face.resize(dofs_per_face);
 
                 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
index e230c033960aad029f422f1c7ef2d854edbcf1a2..c3fae7123800f7927df0d8de300f208b9ed79510 100644 (file)
@@ -145,8 +145,12 @@ namespace internal
         ensure_existence_and_return_dof_identities(
           const FiniteElement<dim, spacedim> &fe1,
           const FiniteElement<dim, spacedim> &fe2,
-          std::unique_ptr<DoFIdentities> &    identities)
+          std::unique_ptr<DoFIdentities> &    identities,
+          const unsigned int face_no = numbers::invalid_unsigned_int)
         {
+          Assert(structdim == 2 || face_no == numbers::invalid_unsigned_int,
+                 ExcInternalError());
+
           // see if we need to fill this entry, or whether it already
           // exists
           if (identities.get() == nullptr)
@@ -170,7 +174,7 @@ namespace internal
                   case 2:
                     {
                       identities = std::make_unique<DoFIdentities>(
-                        fe1.hp_quad_dof_identities(fe2));
+                        fe1.hp_quad_dof_identities(fe2, face_no));
                       break;
                     }
 
@@ -183,10 +187,10 @@ namespace internal
               for (unsigned int i = 0; i < identities->size(); ++i)
                 {
                   Assert((*identities)[i].first <
-                           fe1.template n_dofs_per_object<structdim>(),
+                           fe1.template n_dofs_per_object<structdim>(face_no),
                          ExcInternalError());
                   Assert((*identities)[i].second <
-                           fe2.template n_dofs_per_object<structdim>(),
+                           fe2.template n_dofs_per_object<structdim>(face_no),
                          ExcInternalError());
                 }
             }
@@ -807,6 +811,11 @@ namespace internal
                       fe_indices,
                       /*codim=*/dim - 2);
 
+                  const unsigned int most_dominating_fe_index_face_no =
+                    cell->active_fe_index() == most_dominating_fe_index ?
+                      q :
+                      cell->neighbor_face_no(q);
+
                   // if we found the most dominating element, then use
                   // this to eliminate some of the degrees of freedom
                   // by identification. otherwise, the code that
@@ -829,7 +838,8 @@ namespace internal
                                 quad_dof_identities
                                   [most_dominating_fe_index][other_fe_index]
                                   [cell->quad(q)->reference_cell_type() ==
-                                   ReferenceCell::Type::Quad]);
+                                   ReferenceCell::Type::Quad],
+                                most_dominating_fe_index_face_no);
 
                             for (const auto &identity : identities)
                               {
@@ -1558,6 +1568,11 @@ namespace internal
                       fe_indices,
                       /*codim=*/dim - 2);
 
+                  const unsigned int most_dominating_fe_index_face_no =
+                    cell->active_fe_index() == most_dominating_fe_index ?
+                      q :
+                      cell->neighbor_face_no(q);
+
                   // if we found the most dominating element, then use
                   // this to eliminate some of the degrees of freedom
                   // by identification. otherwise, the code that
@@ -1580,7 +1595,8 @@ namespace internal
                                 quad_dof_identities
                                   [most_dominating_fe_index][other_fe_index]
                                   [cell->quad(q)->reference_cell_type() ==
-                                   ReferenceCell::Type::Quad]);
+                                   ReferenceCell::Type::Quad],
+                                most_dominating_fe_index_face_no);
 
                             for (const auto &identity : identities)
                               {
@@ -2312,7 +2328,7 @@ namespace internal
 
                           for (unsigned int d = 0;
                                d <
-                               dof_handler.get_fe(fe_index).n_dofs_per_quad();
+                               dof_handler.get_fe(fe_index).n_dofs_per_quad(q);
                                ++d)
                             {
                               const types::global_dof_index old_dof_index =
@@ -2674,7 +2690,7 @@ namespace internal
                   if (cell->quad(l)->user_flag_set())
                     {
                       for (unsigned int d = 0;
-                           d < dof_handler.get_fe().n_dofs_per_quad();
+                           d < dof_handler.get_fe().n_dofs_per_quad(l);
                            ++d)
                         {
                           const dealii::types::global_dof_index idx =
index 341edb43965b71a9169851ede2b658728b77af63..fce2cbdfa8e1ba754d2a81a9d23a13efb6881041 100644 (file)
@@ -712,12 +712,12 @@ namespace DoFTools
               {
                 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
 
-                const unsigned int dofs_per_face = fe.n_dofs_per_face();
+                const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
                 face_dof_indices.resize(dofs_per_face);
                 cell->face(face)->get_dof_indices(face_dof_indices,
                                                   cell->active_fe_index());
 
-                for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+                for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                   if (!check_vector_component)
                     selected_dofs.add_index(face_dof_indices[i]);
                   else
@@ -744,9 +744,9 @@ namespace DoFTools
                                           numbers::invalid_unsigned_int)));
                       if (fe.is_primitive(cell_index))
                         {
-                          if (component_mask
-                                [fe.face_system_to_component_index(i).first] ==
-                              true)
+                          if (component_mask[fe.face_system_to_component_index(
+                                                 i, face)
+                                               .first] == true)
                             selected_dofs.add_index(face_dof_indices[i]);
                         }
                       else // not primitive
@@ -1041,7 +1041,7 @@ namespace DoFTools
                         {
                           // simply take all DoFs that live on this subface
                           std::vector<types::global_dof_index> ldi(
-                            fe.n_dofs_per_face());
+                            fe.n_dofs_per_face(f, child));
                           face->child(child)->get_dof_indices(ldi);
                           selected_dofs.add_indices(ldi.begin(), ldi.end());
                         }
@@ -2111,7 +2111,8 @@ namespace DoFTools
       for (const unsigned int f : cell->face_indices())
         if (cell->at_boundary(f))
           {
-            const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+            const unsigned int dofs_per_face =
+              cell->get_fe().n_dofs_per_face(f);
             dofs_on_face.resize(dofs_per_face);
             cell->face(f)->get_dof_indices(dofs_on_face,
                                            cell->active_fe_index());
@@ -2156,7 +2157,8 @@ namespace DoFTools
         if (boundary_ids.find(cell->face(f)->boundary_id()) !=
             boundary_ids.end())
           {
-            const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+            const unsigned int dofs_per_face =
+              cell->get_fe().n_dofs_per_face(f);
             dofs_on_face.resize(dofs_per_face);
             cell->face(f)->get_dof_indices(dofs_on_face,
                                            cell->active_fe_index());
@@ -2533,12 +2535,11 @@ namespace DoFTools
             // Exclude degrees of freedom on faces opposite to the vertex
             exclude.resize(fe.n_dofs_per_cell());
             std::fill(exclude.begin(), exclude.end(), false);
-            const unsigned int dpf = fe.n_dofs_per_face();
 
             for (const unsigned int face : cell->face_indices())
               if (cell->at_boundary(face) ||
                   cell->neighbor(face)->level() != cell->level())
-                for (unsigned int i = 0; i < dpf; ++i)
+                for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                   exclude[fe.face_to_cell_index(i, face)] = true;
             for (types::global_dof_index j = 0; j < indices.size(); ++j)
               if (!exclude[j])
@@ -2594,13 +2595,11 @@ namespace DoFTools
               {
                 // Eliminate dofs on faces of the child which are on faces
                 // of the parent
-                const unsigned int dpf = fe.n_dofs_per_face();
-
                 for (unsigned int d = 0; d < dim; ++d)
                   {
                     const unsigned int face =
                       GeometryInfo<dim>::vertex_to_face[child][d];
-                    for (unsigned int i = 0; i < dpf; ++i)
+                    for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                       exclude[fe.face_to_cell_index(i, face)] = true;
                   }
 
@@ -2610,7 +2609,8 @@ namespace DoFTools
                   for (const unsigned int face :
                        GeometryInfo<dim>::face_indices())
                     if (cell->at_boundary(face))
-                      for (unsigned int i = 0; i < dpf; ++i)
+                      for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
+                           ++i)
                         exclude[fe.face_to_cell_index(i, face)] = false;
               }
 
@@ -2750,7 +2750,6 @@ namespace DoFTools
                 // vertex
                 exclude.resize(fe.n_dofs_per_cell());
                 std::fill(exclude.begin(), exclude.end(), false);
-                const unsigned int dpf = fe.n_dofs_per_face();
 
                 for (unsigned int d = 0; d < dim; ++d)
                   {
@@ -2758,7 +2757,7 @@ namespace DoFTools
                       GeometryInfo<dim>::vertex_to_face[v][d];
                     const unsigned int face =
                       GeometryInfo<dim>::opposite_face[a_face];
-                    for (unsigned int i = 0; i < dpf; ++i)
+                    for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
                       {
                         // For each dof, get the block it is in and decide to
                         // exclude it or not
index 3fd168895cd2531a0b592bba545cf27781e79048..ae4f919a880d04140bdb156b19780a788923f336 100644 (file)
@@ -175,15 +175,23 @@ namespace DoFTools
         const FullMatrix<double> &          face_interpolation_matrix,
         std::vector<bool> &                 primary_dof_mask)
       {
-        Assert(fe1.n_dofs_per_face() >= fe2.n_dofs_per_face(),
+        // TODO: the implementation makes the assumption that all faces have the
+        // same number of dofs
+        AssertDimension(fe1.n_unique_faces(), 1);
+        AssertDimension(fe2.n_unique_faces(), 1);
+        const unsigned int face_no = 0;
+        (void)face_no;
+
+        Assert(fe1.n_dofs_per_face(face_no) >= fe2.n_dofs_per_face(face_no),
                ExcInternalError());
-        AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face());
+        AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face(face_no));
 
         Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(),
                ExcInternalError());
         Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(),
                ExcInternalError());
-        Assert((dim < 3) || (fe2.n_dofs_per_quad() <= fe1.n_dofs_per_quad()),
+        Assert((dim < 3) ||
+                 (fe2.n_dofs_per_quad(face_no) <= fe1.n_dofs_per_quad(face_no)),
                ExcInternalError());
 
         // the idea here is to designate as many DoFs in fe1 per object (vertex,
@@ -268,9 +276,9 @@ namespace DoFTools
             // same algorithm as above
             unsigned int dofs_added = 0;
             unsigned int i          = 0;
-            while (dofs_added < fe2.n_dofs_per_quad())
+            while (dofs_added < fe2.n_dofs_per_quad(q))
               {
-                Assert(i < fe1.n_dofs_per_quad(), ExcInternalError());
+                Assert(i < fe1.n_dofs_per_quad(q), ExcInternalError());
 
                 primary_dof_list.push_back(index + i);
                 if (check_primary_dof_list(face_interpolation_matrix,
@@ -281,11 +289,11 @@ namespace DoFTools
 
                 ++i;
               }
-            index += fe1.n_dofs_per_quad();
+            index += fe1.n_dofs_per_quad(q);
           }
 
-        AssertDimension(index, fe1.n_dofs_per_face());
-        AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face());
+        AssertDimension(index, fe1.n_dofs_per_face(face_no));
+        AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face(face_no));
 
         // finally copy the list into the mask
         std::fill(primary_dof_mask.begin(), primary_dof_mask.end(), false);
@@ -307,10 +315,16 @@ namespace DoFTools
         const FullMatrix<double> &          face_interpolation_matrix,
         std::unique_ptr<std::vector<bool>> &primary_dof_mask)
       {
+        // TODO: the implementation makes the assumption that all faces have the
+        // same number of dofs
+        AssertDimension(fe1.n_unique_faces(), 1);
+        AssertDimension(fe2.n_unique_faces(), 1);
+        const unsigned int face_no = 0;
+
         if (primary_dof_mask == nullptr)
           {
             primary_dof_mask =
-              std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face());
+              std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face(face_no));
             select_primary_dofs_for_face_restriction(fe1,
                                                      fe2,
                                                      face_interpolation_matrix,
@@ -332,12 +346,17 @@ namespace DoFTools
         const FiniteElement<dim, spacedim> & fe2,
         std::unique_ptr<FullMatrix<double>> &matrix)
       {
+        // TODO: the implementation makes the assumption that all faces have the
+        // same number of dofs
+        AssertDimension(fe1.n_unique_faces(), 1);
+        AssertDimension(fe2.n_unique_faces(), 1);
+        const unsigned int face_no = 0;
+
         if (matrix == nullptr)
           {
-            matrix =
-              std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(),
-                                                   fe1.n_dofs_per_face());
-            fe1.get_face_interpolation_matrix(fe2, *matrix);
+            matrix = std::make_unique<FullMatrix<double>>(
+              fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
+            fe1.get_face_interpolation_matrix(fe2, *matrix, face_no);
           }
       }
 
@@ -354,12 +373,20 @@ namespace DoFTools
         const unsigned int                   subface,
         std::unique_ptr<FullMatrix<double>> &matrix)
       {
+        // TODO: the implementation makes the assumption that all faces have the
+        // same number of dofs
+        AssertDimension(fe1.n_unique_faces(), 1);
+        AssertDimension(fe2.n_unique_faces(), 1);
+        const unsigned int face_no = 0;
+
         if (matrix == nullptr)
           {
-            matrix =
-              std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(),
-                                                   fe1.n_dofs_per_face());
-            fe1.get_subface_interpolation_matrix(fe2, subface, *matrix);
+            matrix = std::make_unique<FullMatrix<double>>(
+              fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
+            fe1.get_subface_interpolation_matrix(fe2,
+                                                 subface,
+                                                 *matrix,
+                                                 face_no);
           }
       }
 
@@ -774,7 +801,7 @@ namespace DoFTools
                 // first of all, make sure that we treat a case which is
                 // possible, i.e. either no dofs on the face at all or no
                 // anisotropic refinement
-                if (cell->get_fe().n_dofs_per_face() == 0)
+                if (cell->get_fe().n_dofs_per_face(face) == 0)
                   continue;
 
                 Assert(cell->face(face)->refinement_case() ==
@@ -836,10 +863,10 @@ namespace DoFTools
                 const FiniteElement<dim> &fe       = cell->get_fe();
                 const unsigned int        fe_index = cell->active_fe_index();
 
-                const unsigned int n_dofs_on_mother = fe.n_dofs_per_face();
+                const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
                 const unsigned int n_dofs_on_children =
                   (5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() +
-                   4 * fe.n_dofs_per_quad());
+                   4 * fe.n_dofs_per_quad(face));
 
                 // TODO[TL]: think about this and the following in case of
                 // anisotropic refinement
@@ -872,7 +899,8 @@ namespace DoFTools
                   for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
                     dofs_on_mother[next_index++] =
                       this_face->line(line)->dof_index(dof, fe_index);
-                for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); ++dof)
+                for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
+                     ++dof)
                   dofs_on_mother[next_index++] =
                     this_face->dof_index(dof, fe_index);
                 AssertDimension(next_index, dofs_on_mother.size());
@@ -937,7 +965,7 @@ namespace DoFTools
                     if (cell->neighbor_child_on_subface(face, child)
                           ->is_artificial())
                       continue;
-                    for (unsigned int dof = 0; dof != fe.n_dofs_per_quad();
+                    for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
                          ++dof)
                       dofs_on_children.push_back(
                         this_face->child(child)->dof_index(dof, fe_index));
@@ -1048,7 +1076,7 @@ namespace DoFTools
                 // first of all, make sure that we treat a case which is
                 // possible, i.e. either no dofs on the face at all or no
                 // anisotropic refinement
-                if (cell->get_fe().n_dofs_per_face() == 0)
+                if (cell->get_fe().n_dofs_per_face(face) == 0)
                   continue;
 
                 Assert(cell->face(face)->refinement_case() ==
@@ -1120,7 +1148,8 @@ namespace DoFTools
                         //
                         // so we are going to constrain the DoFs on the face
                         // children against the DoFs on the face itself
-                        primary_dofs.resize(cell->get_fe().n_dofs_per_face());
+                        primary_dofs.resize(
+                          cell->get_fe().n_dofs_per_face(face));
 
                         cell->face(face)->get_dof_indices(
                           primary_dofs, cell->active_fe_index());
@@ -1166,7 +1195,7 @@ namespace DoFTools
                             // the face DoFs from the cell DoFs.
                             dependent_dofs.resize(
                               subface->get_fe(subface_fe_index)
-                                .n_dofs_per_face());
+                                .n_dofs_per_face(face, c));
                             subface->get_dof_indices(dependent_dofs,
                                                      subface_fe_index);
 
@@ -1275,8 +1304,8 @@ namespace DoFTools
 
                         // first get the interpolation matrix from the mother to
                         // the virtual dofs
-                        Assert(dominating_fe.n_dofs_per_face() <=
-                                 cell->get_fe().n_dofs_per_face(),
+                        Assert(dominating_fe.n_dofs_per_face(face) <=
+                                 cell->get_fe().n_dofs_per_face(face),
                                ExcInternalError());
 
                         ensure_existence_of_face_matrix(
@@ -1318,16 +1347,17 @@ namespace DoFTools
                         // now compute the constraint matrix as the product
                         // between the inverse matrix and the dependent part
                         constraint_matrix.reinit(
-                          cell->get_fe().n_dofs_per_face() -
-                            dominating_fe.n_dofs_per_face(),
-                          dominating_fe.n_dofs_per_face());
+                          cell->get_fe().n_dofs_per_face(face) -
+                            dominating_fe.n_dofs_per_face(face),
+                          dominating_fe.n_dofs_per_face(face));
                         restrict_mother_to_virtual_dependent.mmult(
                           constraint_matrix,
                           restrict_mother_to_virtual_primary_inv);
 
                         // then figure out the global numbers of primary and
                         // dependent dofs and apply constraints
-                        scratch_dofs.resize(cell->get_fe().n_dofs_per_face());
+                        scratch_dofs.resize(
+                          cell->get_fe().n_dofs_per_face(face));
                         cell->face(face)->get_dof_indices(
                           scratch_dofs, cell->active_fe_index());
 
@@ -1335,7 +1365,7 @@ namespace DoFTools
                         primary_dofs.clear();
                         dependent_dofs.clear();
                         for (unsigned int i = 0;
-                             i < cell->get_fe().n_dofs_per_face();
+                             i < cell->get_fe().n_dofs_per_face(face);
                              ++i)
                           if ((*primary_dof_masks[dominating_fe_index]
                                                  [cell
@@ -1346,10 +1376,10 @@ namespace DoFTools
                             dependent_dofs.push_back(scratch_dofs[i]);
 
                         AssertDimension(primary_dofs.size(),
-                                        dominating_fe.n_dofs_per_face());
+                                        dominating_fe.n_dofs_per_face(face));
                         AssertDimension(dependent_dofs.size(),
-                                        cell->get_fe().n_dofs_per_face() -
-                                          dominating_fe.n_dofs_per_face());
+                                        cell->get_fe().n_dofs_per_face(face) -
+                                          dominating_fe.n_dofs_per_face(face));
 
                         filter_constraints(primary_dofs,
                                            dependent_dofs,
@@ -1386,8 +1416,8 @@ namespace DoFTools
 
                             // first get the interpolation matrix from the
                             // subface to the virtual dofs
-                            Assert(dominating_fe.n_dofs_per_face() <=
-                                     subface_fe.n_dofs_per_face(),
+                            Assert(dominating_fe.n_dofs_per_face(face) <=
+                                     subface_fe.n_dofs_per_face(face),
                                    ExcInternalError());
                             ensure_existence_of_subface_matrix(
                               dominating_fe,
@@ -1402,14 +1432,15 @@ namespace DoFTools
                                   [dominating_fe_index][subface_fe_index][sf]);
 
                             constraint_matrix.reinit(
-                              subface_fe.n_dofs_per_face(),
-                              dominating_fe.n_dofs_per_face());
+                              subface_fe.n_dofs_per_face(face),
+                              dominating_fe.n_dofs_per_face(face));
 
                             restrict_subface_to_virtual.mmult(
                               constraint_matrix,
                               restrict_mother_to_virtual_primary_inv);
 
-                            dependent_dofs.resize(subface_fe.n_dofs_per_face());
+                            dependent_dofs.resize(
+                              subface_fe.n_dofs_per_face(face));
                             cell->face(face)->child(sf)->get_dof_indices(
                               dependent_dofs, subface_fe_index);
 
@@ -1471,7 +1502,7 @@ namespace DoFTools
                             // Get DoFs on dominating and dominated side of the
                             // face
                             primary_dofs.resize(
-                              cell->get_fe().n_dofs_per_face());
+                              cell->get_fe().n_dofs_per_face(face));
                             cell->face(face)->get_dof_indices(
                               primary_dofs, cell->active_fe_index());
 
@@ -1482,7 +1513,7 @@ namespace DoFTools
                               break;
 
                             dependent_dofs.resize(
-                              neighbor->get_fe().n_dofs_per_face());
+                              neighbor->get_fe().n_dofs_per_face(face));
                             cell->face(face)->get_dof_indices(
                               dependent_dofs, neighbor->active_fe_index());
 
@@ -1589,8 +1620,8 @@ namespace DoFTools
 
                             // first get the interpolation matrix from main FE
                             // to the virtual dofs
-                            Assert(dominating_fe.n_dofs_per_face() <=
-                                     cell->get_fe().n_dofs_per_face(),
+                            Assert(dominating_fe.n_dofs_per_face(face) <=
+                                     cell->get_fe().n_dofs_per_face(face),
                                    ExcInternalError());
 
                             ensure_existence_of_face_matrix(
@@ -1633,9 +1664,9 @@ namespace DoFTools
                             // now compute the constraint matrix as the product
                             // between the inverse matrix and the dependent part
                             constraint_matrix.reinit(
-                              cell->get_fe().n_dofs_per_face() -
-                                dominating_fe.n_dofs_per_face(),
-                              dominating_fe.n_dofs_per_face());
+                              cell->get_fe().n_dofs_per_face(face) -
+                                dominating_fe.n_dofs_per_face(face),
+                              dominating_fe.n_dofs_per_face(face));
                             restrict_mother_to_virtual_dependent.mmult(
                               constraint_matrix,
                               restrict_mother_to_virtual_primary_inv);
@@ -1643,7 +1674,7 @@ namespace DoFTools
                             // then figure out the global numbers of primary and
                             // dependent dofs and apply constraints
                             scratch_dofs.resize(
-                              cell->get_fe().n_dofs_per_face());
+                              cell->get_fe().n_dofs_per_face(face));
                             cell->face(face)->get_dof_indices(
                               scratch_dofs, cell->active_fe_index());
 
@@ -1651,7 +1682,7 @@ namespace DoFTools
                             primary_dofs.clear();
                             dependent_dofs.clear();
                             for (unsigned int i = 0;
-                                 i < cell->get_fe().n_dofs_per_face();
+                                 i < cell->get_fe().n_dofs_per_face(face);
                                  ++i)
                               if ((*primary_dof_masks[dominating_fe_index]
                                                      [cell->active_fe_index()])
@@ -1661,10 +1692,12 @@ namespace DoFTools
                                 dependent_dofs.push_back(scratch_dofs[i]);
 
                             AssertDimension(primary_dofs.size(),
-                                            dominating_fe.n_dofs_per_face());
-                            AssertDimension(dependent_dofs.size(),
-                                            cell->get_fe().n_dofs_per_face() -
-                                              dominating_fe.n_dofs_per_face());
+                                            dominating_fe.n_dofs_per_face(
+                                              face));
+                            AssertDimension(
+                              dependent_dofs.size(),
+                              cell->get_fe().n_dofs_per_face(face) -
+                                dominating_fe.n_dofs_per_face(face));
 
                             filter_constraints(primary_dofs,
                                                dependent_dofs,
@@ -1674,8 +1707,8 @@ namespace DoFTools
                             // now do the same for another FE this is pretty
                             // much the same we do above to resolve h-refinement
                             // constraints
-                            Assert(dominating_fe.n_dofs_per_face() <=
-                                     neighbor->get_fe().n_dofs_per_face(),
+                            Assert(dominating_fe.n_dofs_per_face(face) <=
+                                     neighbor->get_fe().n_dofs_per_face(face),
                                    ExcInternalError());
 
                             ensure_existence_of_face_matrix(
@@ -1692,15 +1725,15 @@ namespace DoFTools
                                     [neighbor->active_fe_index()]);
 
                             constraint_matrix.reinit(
-                              neighbor->get_fe().n_dofs_per_face(),
-                              dominating_fe.n_dofs_per_face());
+                              neighbor->get_fe().n_dofs_per_face(face),
+                              dominating_fe.n_dofs_per_face(face));
 
                             restrict_secondface_to_virtual.mmult(
                               constraint_matrix,
                               restrict_mother_to_virtual_primary_inv);
 
                             dependent_dofs.resize(
-                              neighbor->get_fe().n_dofs_per_face());
+                              neighbor->get_fe().n_dofs_per_face(face));
                             cell->face(face)->get_dof_indices(
                               dependent_dofs, neighbor->active_fe_index());
 
@@ -1796,6 +1829,14 @@ namespace DoFTools
 
       Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
 
+      // TODO: the implementation makes the assumption that all faces have the
+      // same number of dofs
+      AssertDimension(
+        face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
+      AssertDimension(
+        face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
+      const unsigned int face_no = 0;
+
       // If face_2 does have children, then we need to iterate over these
       // children and set periodic constraints in the inverse direction:
 
@@ -1806,7 +1847,8 @@ namespace DoFTools
                  ExcNotImplemented());
 
           const unsigned int dofs_per_face =
-            face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face();
+            face_1->get_fe(face_1->nth_active_fe_index(0))
+              .n_dofs_per_face(face_no);
           FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
           FullMatrix<double> subface_interpolation(dofs_per_face,
                                                    dofs_per_face);
@@ -1817,7 +1859,10 @@ namespace DoFTools
               // interpolated from face_1 to face_2 by multiplying from the left
               // with the one that interpolates from face_2 to its child
               const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
-              fe.get_subface_interpolation_matrix(fe, c, subface_interpolation);
+              fe.get_subface_interpolation_matrix(fe,
+                                                  c,
+                                                  subface_interpolation,
+                                                  face_no);
               subface_interpolation.mmult(child_transformation, transformation);
 
               set_periodicity_constraints(face_1,
@@ -1852,7 +1897,7 @@ namespace DoFTools
                "zero or equal to the number of components in the finite "
                "element."));
 
-      const unsigned int dofs_per_face = fe.n_dofs_per_face();
+      const unsigned int dofs_per_face = fe.n_dofs_per_face(face_no);
 
       std::vector<types::global_dof_index> dofs_1(dofs_per_face);
       std::vector<types::global_dof_index> dofs_2(dofs_per_face);
@@ -1925,7 +1970,8 @@ namespace DoFTools
           // Obey the component mask
           if ((component_mask.n_selected_components(fe.n_components()) !=
                fe.n_components()) &&
-              !component_mask[fe.face_system_to_component_index(i).first])
+              !component_mask[fe.face_system_to_component_index(i, face_no)
+                                .first])
             continue;
 
           // We have to be careful to treat so called "identity
@@ -2120,9 +2166,14 @@ namespace DoFTools
       const FullMatrix<double> &          matrix,
       const std::vector<unsigned int> &   first_vector_components)
     {
+      // TODO: the implementation makes the assumption that all faces have the
+      // same number of dofs
+      AssertDimension(fe.n_unique_faces(), 1);
+      const unsigned int face_no = 0;
+
       Assert(matrix.m() == matrix.n(), ExcInternalError());
 
-      const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
+      const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
 
       if (matrix.m() == n_dofs_per_face)
         {
@@ -2144,7 +2195,7 @@ namespace DoFTools
       Assert(matrix.m() == spacedim, ExcInternalError())
 
         Quadrature<dim - 1>
-          quadrature(fe.get_unit_face_support_points());
+          quadrature(fe.get_unit_face_support_points(face_no));
 
       // have an array that stores the location of each vector-dof tuple we want
       // to rotate.
@@ -2158,7 +2209,7 @@ namespace DoFTools
           std::vector<unsigned int>::const_iterator comp_it =
             std::find(first_vector_components.begin(),
                       first_vector_components.end(),
-                      fe.face_system_to_component_index(i).first);
+                      fe.face_system_to_component_index(i, face_no).first);
           if (comp_it != first_vector_components.end())
             {
               const unsigned int first_vector_component = *comp_it;
@@ -2176,12 +2227,13 @@ namespace DoFTools
 
               for (unsigned int k = 0; k < n_dofs_per_face; ++k)
                 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
-                    (fe.face_system_to_component_index(k).first >=
+                    (fe.face_system_to_component_index(k, face_no).first >=
                      first_vector_component) &&
-                    (fe.face_system_to_component_index(k).first <
+                    (fe.face_system_to_component_index(k, face_no).first <
                      first_vector_component + spacedim))
                   {
-                    vector_dofs[fe.face_system_to_component_index(k).first -
+                    vector_dofs[fe.face_system_to_component_index(k, face_no)
+                                  .first -
                                 first_vector_component] = k;
                     n_found++;
                     if (n_found == dim)
@@ -2221,6 +2273,14 @@ namespace DoFTools
     const std::vector<unsigned int> &            first_vector_components,
     const number                                 periodicity_factor)
   {
+    // TODO: the implementation makes the assumption that all faces have the
+    // same number of dofs
+    AssertDimension(
+      face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
+    AssertDimension(
+      face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
+    const unsigned int face_no = 0;
+
     static const int dim      = FaceIterator::AccessorType::dimension;
     static const int spacedim = FaceIterator::AccessorType::space_dimension;
 
@@ -2256,7 +2316,8 @@ namespace DoFTools
       {
         Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
         const unsigned int n_dofs_per_face =
-          face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face();
+          face_1->get_fe(face_1->nth_active_fe_index(0))
+            .n_dofs_per_face(face_no);
 
         Assert(matrix.m() == 0 ||
                  (first_vector_components.empty() &&
@@ -2273,7 +2334,8 @@ namespace DoFTools
       {
         Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
         const unsigned int n_dofs_per_face =
-          face_2->get_fe(face_2->nth_active_fe_index(0)).n_dofs_per_face();
+          face_2->get_fe(face_2->nth_active_fe_index(0))
+            .n_dofs_per_face(face_no);
 
         Assert(matrix.m() == 0 ||
                  (first_vector_components.empty() &&
@@ -2372,7 +2434,7 @@ namespace DoFTools
             face_2->get_fe(face_2->nth_active_fe_index(0)) :
             face_1->get_fe(face_1->nth_active_fe_index(0));
 
-        const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
+        const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
 
         // Sometimes we just have nothing to do (for all finite elements, or
         // systems which accidentally don't have any dofs on the boundary).
@@ -3408,7 +3470,7 @@ namespace DoFTools
                    (face->boundary_id() == boundary_id)))
                 {
                   // get indices and physical location on this face
-                  face_dofs.resize(fe.n_dofs_per_face());
+                  face_dofs.resize(fe.n_dofs_per_face(face_no));
                   face->get_dof_indices(face_dofs, cell->active_fe_index());
 
                   // enter those dofs into the list that match the component
index 165a15ca663d8d424e4121455271f725266762ea..c61801583d30b9fff7e2b9b3e02b00e85df4a2e9 100644 (file)
@@ -407,7 +407,8 @@ namespace DoFTools
       for (const unsigned int f : cell->face_indices())
         if (cell->at_boundary(f))
           {
-            const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+            const unsigned int dofs_per_face =
+              cell->get_fe().n_dofs_per_face(f);
             dofs_on_this_face.resize(dofs_per_face);
             cell->face(f)->get_dof_indices(dofs_on_this_face,
                                            cell->active_fe_index());
@@ -504,7 +505,8 @@ namespace DoFTools
         if (boundary_ids.find(cell->face(f)->boundary_id()) !=
             boundary_ids.end())
           {
-            const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+            const unsigned int dofs_per_face =
+              cell->get_fe().n_dofs_per_face(f);
             dofs_on_this_face.resize(dofs_per_face);
             cell->face(f)->get_dof_indices(dofs_on_this_face,
                                            cell->active_fe_index());
index 88598a00cf377362456441602b955bdf89dd49d9..a4d99be7b8e403d21706680c688a4eb5780cb1f6 100644 (file)
@@ -1830,7 +1830,7 @@ namespace GridTools
                                ExcInternalError());
 
                         const unsigned int n_dofs_per_face =
-                          cell->get_fe().n_dofs_per_face();
+                          cell->get_fe().n_dofs_per_face(f, c);
                         local_face_dof_indices.resize(n_dofs_per_face);
 
                         cell->face(f)->child(c)->get_dof_indices(
@@ -1865,7 +1865,7 @@ namespace GridTools
                     unsigned int subface = neighbor_face_no_subface_no.second;
 
                     const unsigned int n_dofs_per_face =
-                      cell->get_fe().n_dofs_per_face();
+                      cell->get_fe().n_dofs_per_face(face_no);
                     local_face_dof_indices.resize(n_dofs_per_face);
 
                     cell->neighbor(f)->face(face_no)->get_dof_indices(
@@ -1884,7 +1884,7 @@ namespace GridTools
                                           // original cell
                           {
                             const unsigned int n_dofs_per_face =
-                              cell->get_fe().n_dofs_per_face();
+                              cell->get_fe().n_dofs_per_face(face_no, c);
                             local_face_dof_indices.resize(n_dofs_per_face);
 
                             Assert(cell->neighbor(f)
index a1f65565175749f103203efd53c8577f8924733e..09d532ae6ff2f64c811f37e767ecc8c6ff05ca1d 100644 (file)
@@ -162,10 +162,14 @@ namespace MGTools
         // n_dofs_per_vertex() is
         // arbitrary, not the other way
         // round.
-        // TODO: This assumes that even in hp context, the dofs per face
-        // coincide!
+        // TODO: This assumes that the dofs per face on all faces coincide!
+        const unsigned int face_no = 0;
+
+        Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim),
+               ExcNotImplemented());
+
         unsigned int increment =
-          fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face();
+          fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face(face_no);
         while (i < fe.get_first_line_index())
           row_lengths[cell_indices[i++]] += increment;
         // From now on, if an object is
@@ -178,23 +182,25 @@ namespace MGTools
         // In all other cases we
         // subtract adjacent faces to be
         // added in the loop below.
-        increment = (dim > 1) ?
-                      fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face() :
-                      fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
-                                               fe.n_dofs_per_face();
-        while (i < fe.get_first_quad_index())
+        increment =
+          (dim > 1) ?
+            fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face(face_no) :
+            fe.n_dofs_per_cell() -
+              GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face(face_no);
+        while (i < fe.get_first_quad_index(face_no))
           row_lengths[cell_indices[i++]] += increment;
 
         // Now quads in 2D and 3D
-        increment = (dim > 2) ?
-                      fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face() :
-                      fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
-                                               fe.n_dofs_per_face();
+        increment =
+          (dim > 2) ?
+            fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face(face_no) :
+            fe.n_dofs_per_cell() -
+              GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face(face_no);
         while (i < fe.get_first_hex_index())
           row_lengths[cell_indices[i++]] += increment;
         // Finally, cells in 3D
-        increment = fe.n_dofs_per_cell() -
-                    GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face();
+        increment = fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
+                                             fe.n_dofs_per_face(face_no);
         while (i < fe.n_dofs_per_cell())
           row_lengths[cell_indices[i++]] += increment;
 
@@ -225,7 +231,8 @@ namespace MGTools
                 for (unsigned int local_dof = 0;
                      local_dof < fe.n_dofs_per_cell();
                      ++local_dof)
-                  row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face();
+                  row_lengths[cell_indices[local_dof]] +=
+                    fe.n_dofs_per_face(face_no);
                 continue;
               }
 
@@ -244,7 +251,7 @@ namespace MGTools
             if (flux_coupling != DoFTools::none)
               {
                 const unsigned int dof_increment =
-                  nfe.n_dofs_per_cell() - nfe.n_dofs_per_face();
+                  nfe.n_dofs_per_cell() - nfe.n_dofs_per_face(face_no);
                 for (unsigned int local_dof = 0;
                      local_dof < fe.n_dofs_per_cell();
                      ++local_dof)
@@ -272,10 +279,12 @@ namespace MGTools
             neighbor->get_mg_dof_indices(neighbor_indices);
             for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell();
                  ++local_dof)
-              row_lengths[cell_indices[local_dof]] += nfe.n_dofs_per_face();
+              row_lengths[cell_indices[local_dof]] +=
+                nfe.n_dofs_per_face(face_no);
             for (unsigned int local_dof = 0; local_dof < nfe.n_dofs_per_cell();
                  ++local_dof)
-              row_lengths[neighbor_indices[local_dof]] += fe.n_dofs_per_face();
+              row_lengths[neighbor_indices[local_dof]] +=
+                fe.n_dofs_per_face(face_no);
           }
       }
     user_flags_triangulation.load_user_flags(old_flags);
@@ -333,6 +342,12 @@ namespace MGTools
         const FiniteElement<dim> &fe       = cell->get_fe();
         const unsigned int        fe_index = cell->active_fe_index();
 
+
+        // TODO: This assumes that the dofs per face on all faces coincide!
+        const unsigned int face_no = 0;
+        Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim),
+               ExcNotImplemented());
+
         Assert(couplings.n_rows() == fe.n_components(),
                ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
         Assert(couplings.n_cols() == fe.n_components(),
@@ -378,8 +393,9 @@ namespace MGTools
                                           fe.first_block_of_base(base) +
                                             mult) != DoFTools::none)
                   {
-                    increment = fe.base_element(base).n_dofs_per_cell() -
-                                dim * fe.base_element(base).n_dofs_per_face();
+                    increment =
+                      fe.base_element(base).n_dofs_per_cell() -
+                      dim * fe.base_element(base).n_dofs_per_face(face_no);
                     row_lengths[cell_indices[i]] += increment;
                   }
             ++i;
@@ -394,7 +410,7 @@ namespace MGTools
         // In all other cases we
         // subtract adjacent faces to be
         // added in the loop below.
-        while (i < fe.get_first_quad_index())
+        while (i < fe.get_first_quad_index(face_no))
           {
             for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
               for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
@@ -407,7 +423,7 @@ namespace MGTools
                       fe.base_element(base).n_dofs_per_cell() -
                       ((dim > 1) ? (dim - 1) :
                                    GeometryInfo<dim>::faces_per_cell) *
-                        fe.base_element(base).n_dofs_per_face();
+                        fe.base_element(base).n_dofs_per_face(face_no);
                     row_lengths[cell_indices[i]] += increment;
                   }
             ++i;
@@ -427,7 +443,7 @@ namespace MGTools
                       fe.base_element(base).n_dofs_per_cell() -
                       ((dim > 2) ? (dim - 2) :
                                    GeometryInfo<dim>::faces_per_cell) *
-                        fe.base_element(base).n_dofs_per_face();
+                        fe.base_element(base).n_dofs_per_face(face_no);
                     row_lengths[cell_indices[i]] += increment;
                   }
             ++i;
@@ -443,9 +459,10 @@ namespace MGTools
                                           fe.first_block_of_base(base) +
                                             mult) != DoFTools::none)
                   {
-                    increment = fe.base_element(base).n_dofs_per_cell() -
-                                GeometryInfo<dim>::faces_per_cell *
-                                  fe.base_element(base).n_dofs_per_face();
+                    increment =
+                      fe.base_element(base).n_dofs_per_cell() -
+                      GeometryInfo<dim>::faces_per_cell *
+                        fe.base_element(base).n_dofs_per_face(face_no);
                     row_lengths[cell_indices[i]] += increment;
                   }
             ++i;
@@ -478,7 +495,8 @@ namespace MGTools
                 for (unsigned int local_dof = 0;
                      local_dof < fe.n_dofs_per_cell();
                      ++local_dof)
-                  row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face();
+                  row_lengths[cell_indices[local_dof]] +=
+                    fe.n_dofs_per_face(face_no);
                 continue;
               }
 
@@ -506,7 +524,7 @@ namespace MGTools
                     {
                       const unsigned int dof_increment =
                         nfe.base_element(base).n_dofs_per_cell() -
-                        nfe.base_element(base).n_dofs_per_face();
+                        nfe.base_element(base).n_dofs_per_face(face_no);
                       row_lengths[cell_indices[local_dof]] += dof_increment;
                     }
 
@@ -548,7 +566,7 @@ namespace MGTools
                         fe.system_to_component_index(local_dof).first,
                         nfe.first_block_of_base(base) + mult) != DoFTools::none)
                     row_lengths[cell_indices[local_dof]] +=
-                      nfe.base_element(base).n_dofs_per_face();
+                      nfe.base_element(base).n_dofs_per_face(face_no);
             for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
               for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                    ++mult)
@@ -559,7 +577,7 @@ namespace MGTools
                         nfe.system_to_component_index(local_dof).first,
                         fe.first_block_of_base(base) + mult) != DoFTools::none)
                     row_lengths[neighbor_indices[local_dof]] +=
-                      fe.base_element(base).n_dofs_per_face();
+                      fe.base_element(base).n_dofs_per_face(face_no);
           }
       }
     user_flags_triangulation.load_user_flags(old_flags);
@@ -1332,7 +1350,6 @@ namespace MGTools
               continue;
             const FiniteElement<dim> &fe    = cell->get_fe();
             const unsigned int        level = cell->level();
-            local_dofs.resize(fe.n_dofs_per_face());
 
             for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
               if (cell->at_boundary(face_no) == true)
@@ -1343,6 +1360,7 @@ namespace MGTools
                   // Face is listed in boundary map
                   if (boundary_ids.find(bi) != boundary_ids.end())
                     {
+                      local_dofs.resize(fe.n_dofs_per_face(face_no));
                       face->get_mg_dof_indices(level, local_dofs);
                       dofs_by_level[level].insert(dofs_by_level[level].end(),
                                                   local_dofs.begin(),
@@ -1405,7 +1423,7 @@ namespace MGTools
 
                     // get indices, physical location and boundary values of
                     // dofs on this face
-                    local_dofs.resize(fe.n_dofs_per_face());
+                    local_dofs.resize(fe.n_dofs_per_face(face_no));
                     face->get_mg_dof_indices(level, local_dofs);
                     if (fe_is_system)
                       {
@@ -1415,7 +1433,8 @@ namespace MGTools
                               numbers::invalid_unsigned_int;
                             if (fe.is_primitive())
                               component =
-                                fe.face_system_to_component_index(i).first;
+                                fe.face_system_to_component_index(i, face_no)
+                                  .first;
                             else
                               {
                                 // Just pick the first of the components
@@ -1473,7 +1492,6 @@ namespace MGTools
     const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
 
     const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-    const unsigned int dofs_per_face = fe.n_dofs_per_face();
 
     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
@@ -1517,7 +1535,8 @@ namespace MGTools
                 // Do refinement face from the coarse side
                 if (neighbor->level() < cell->level())
                   {
-                    for (unsigned int j = 0; j < dofs_per_face; ++j)
+                    for (unsigned int j = 0; j < fe.n_dofs_per_face(face_nr);
+                         ++j)
                       cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
 
                     has_coarser_neighbor = true;

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.