]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor HangingNodes::setup_constraints() 12977/head
authorPeter Munch <peterrmuench@gmail.com>
Thu, 18 Nov 2021 21:57:20 +0000 (22:57 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 23 Nov 2021 10:52:06 +0000 (11:52 +0100)
include/deal.II/matrix_free/cuda_matrix_free.templates.h
include/deal.II/matrix_free/dof_info.h
include/deal.II/matrix_free/dof_info.templates.h
include/deal.II/matrix_free/hanging_nodes_internal.h
include/deal.II/matrix_free/matrix_free.templates.h
source/matrix_free/dof_info.cc
tests/matrix_free/mixed_mesh_hn_algorithm_01.cc [new file with mode: 0644]
tests/matrix_free/mixed_mesh_hn_algorithm_01.output [new file with mode: 0644]

index 4c64f02188259ed4d47f9aa91c94b0882b946197..9efc73c5500467fba9c37ab0e04fd7d87cba5749 100644 (file)
@@ -387,7 +387,7 @@ namespace CUDAWrappers
 
       hanging_nodes.setup_constraints(cell,
                                       partitioner,
-                                      lexicographic_inv,
+                                      {lexicographic_inv},
                                       lexicographic_dof_indices,
                                       cell_id_view);
 
index d4b69bce4c4285e71cd382b4ecf2322f8e9ca160..08178f7eaba6d51b2625c6d09bd8ed12d4116d48 100644 (file)
@@ -222,9 +222,9 @@ namespace internal
       template <int dim>
       bool
       process_hanging_node_constraints(
-        const HangingNodes<dim> &        hanging_nodes,
-        const std::vector<unsigned int> &lexicographic_mapping,
-        const unsigned int               cell_number,
+        const HangingNodes<dim> &                     hanging_nodes,
+        const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+        const unsigned int                            cell_number,
         const TriaIterator<DoFCellAccessor<dim, dim, false>> &cell,
         std::vector<types::global_dof_index> &                dof_indices);
 
index ffb59ed4c3770154e16301dbe2d45e2ed9cd6590..814988d4d79471b8e78d88acff1ae2ff63065167 100644 (file)
@@ -277,9 +277,9 @@ namespace internal
     template <int dim>
     bool
     DoFInfo::process_hanging_node_constraints(
-      const HangingNodes<dim> &        hanging_nodes,
-      const std::vector<unsigned int> &lexicographic_mapping,
-      const unsigned int               cell_number,
+      const HangingNodes<dim> &                     hanging_nodes,
+      const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+      const unsigned int                            cell_number,
       const TriaIterator<DoFCellAccessor<dim, dim, false>> &cell,
       std::vector<types::global_dof_index> &                dof_indices)
     {
index a5bb9c42ebe84e1fe1651053ed2fbf7ed7933a3c..2c3313292b529bb3153871a9a825a91fa5c5ecb8 100644 (file)
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/dofs/dof_accessor.h>
@@ -208,9 +209,38 @@ namespace internal
       setup_constraints(
         const CellIterator &                                      cell,
         const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
-        const std::vector<unsigned int> &     lexicographic_mapping,
-        std::vector<types::global_dof_index> &dof_indices,
-        const ArrayView<ConstraintKinds> &    mask) const;
+        const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+        std::vector<types::global_dof_index> &        dof_indices,
+        const ArrayView<ConstraintKinds> &            mask) const;
+
+      /**
+       * Compute the supported components of all entries of the given
+       * hp::FECollection object.
+       */
+      std::vector<std::vector<bool>>
+      compute_supported_components(
+        const dealii::hp::FECollection<dim> &fe) const;
+
+      /**
+       * Determine the refinement configuration of the given cell.
+       */
+      template <typename CellIterator>
+      ConstraintKinds
+      compute_refinement_configuration(const CellIterator &cell) const;
+
+      /**
+       * Update the DoF indices of a given cell for the given refinement
+       * configuration and for the given components.
+       */
+      template <typename CellIterator>
+      void
+      update_dof_indices(
+        const CellIterator &                                      cell,
+        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
+        const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+        const std::vector<std::vector<bool>> &        component_mask,
+        const ConstraintKinds &                       refinement_configuration,
+        std::vector<types::global_dof_index> &        dof_indices) const;
 
     private:
       /**
@@ -242,6 +272,11 @@ namespace internal
       std::vector<std::vector<
         std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>>
         line_to_cells;
+
+      const dealii::ndarray<unsigned int, 3, 2, 2> local_lines = {
+        {{{{{7, 3}}, {{6, 2}}}},
+         {{{{5, 1}}, {{4, 0}}}},
+         {{{{11, 9}}, {{10, 8}}}}}};
     };
 
 
@@ -345,24 +380,142 @@ namespace internal
 
 
 
+    template <int dim>
+    inline std::vector<std::vector<bool>>
+    HangingNodes<dim>::compute_supported_components(
+      const dealii::hp::FECollection<dim> &fe_collection) const
+    {
+      std::vector<std::vector<bool>> supported_components(
+        fe_collection.size(),
+        std::vector<bool>(fe_collection.n_components(), false));
+
+      for (unsigned int i = 0; i < fe_collection.size(); ++i)
+        {
+          for (unsigned int base_element_index = 0, comp = 0;
+               base_element_index < fe_collection[i].n_base_elements();
+               ++base_element_index)
+            for (unsigned int c = 0;
+                 c < fe_collection[i].element_multiplicity(base_element_index);
+                 ++c, ++comp)
+              if (dim == 1 || dynamic_cast<const FE_Q<dim> *>(
+                                &fe_collection[i].base_element(
+                                  base_element_index)) == nullptr)
+                supported_components[i][comp] = false;
+              else
+                supported_components[i][comp] = true;
+        }
+
+      return supported_components;
+    }
+
+
+
     template <int dim>
     template <typename CellIterator>
-    inline bool
-    HangingNodes<dim>::setup_constraints(
+    inline ConstraintKinds
+    HangingNodes<dim>::compute_refinement_configuration(
+      const CellIterator &cell) const
+    {
+      // TODO: for simplex or mixed meshes: nothing to do
+      if ((dim == 3 && line_to_cells.size() == 0) ||
+          (cell->reference_cell().is_hyper_cube() == false))
+        return ConstraintKinds::unconstrained;
+
+      if (cell->level() == 0)
+        return ConstraintKinds::unconstrained;
+
+      const std::uint16_t subcell =
+        cell->parent()->child_iterator_to_index(cell);
+      const std::uint16_t subcell_x = (subcell >> 0) & 1;
+      const std::uint16_t subcell_y = (subcell >> 1) & 1;
+      const std::uint16_t subcell_z = (subcell >> 2) & 1;
+
+      std::uint16_t face = 0;
+      std::uint16_t edge = 0;
+
+      for (unsigned int direction = 0; direction < dim; ++direction)
+        {
+          const auto side    = (subcell >> direction) & 1U;
+          const auto face_no = direction * 2 + side;
+
+          // ignore if at boundary
+          if (cell->at_boundary(face_no))
+            continue;
+
+          const auto &neighbor = cell->neighbor(face_no);
+
+          // ignore neighbors that are artificial or have the same level or
+          // have children
+          if (neighbor->has_children() || neighbor->is_artificial() ||
+              neighbor->level() == cell->level())
+            continue;
+
+          face |= 1 << direction;
+        }
+
+      if (dim == 3)
+        for (unsigned int direction = 0; direction < dim; ++direction)
+          if (face == 0 || face == (1 << direction))
+            {
+              const unsigned int line_no =
+                direction == 0 ?
+                  (local_lines[0][subcell_y == 0][subcell_z == 0]) :
+                  (direction == 1 ?
+                     (local_lines[1][subcell_x == 0][subcell_z == 0]) :
+                     (local_lines[2][subcell_x == 0][subcell_y == 0]));
+
+              const unsigned int line_index = cell->line(line_no)->index();
+
+              const auto edge_neighbor =
+                std::find_if(line_to_cells[line_index].begin(),
+                             line_to_cells[line_index].end(),
+                             [&cell](const auto &edge_neighbor) {
+                               return edge_neighbor.first->is_artificial() ==
+                                        false &&
+                                      edge_neighbor.first->level() <
+                                        cell->level();
+                             });
+
+              if (edge_neighbor == line_to_cells[line_index].end())
+                continue;
+
+              edge |= 1 << direction;
+            }
+
+      if ((face == 0) && (edge == 0))
+        return ConstraintKinds::unconstrained;
+
+      const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7));
+
+      const auto refinement_configuration = static_cast<ConstraintKinds>(
+        inverted_subcell + (face << 3) + (edge << 6));
+      Assert(check(refinement_configuration, dim), ExcInternalError());
+      return refinement_configuration;
+    }
+
+
+
+    template <int dim>
+    template <typename CellIterator>
+    inline void
+    HangingNodes<dim>::update_dof_indices(
       const CellIterator &                                      cell,
       const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
-      const std::vector<unsigned int> &     lexicographic_mapping,
-      std::vector<types::global_dof_index> &dof_indices,
-      const ArrayView<ConstraintKinds> &    masks) const
+      const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+      const std::vector<std::vector<bool>> &        supported_components,
+      const ConstraintKinds &                       refinement_configuration,
+      std::vector<types::global_dof_index> &        dof_indices) const
     {
-      bool cell_has_hanging_node_constraints = false;
-
-      // for simplex or mixed meshes: nothing to do
-      if (dim == 3 && line_to_cells.size() == 0)
-        return cell_has_hanging_node_constraints;
+      if (std::find(supported_components[cell->active_fe_index()].begin(),
+                    supported_components[cell->active_fe_index()].end(),
+                    true) ==
+          supported_components[cell->active_fe_index()].end())
+        return;
 
       const auto &fe = cell->get_fe();
 
+      AssertDimension(fe.n_unique_faces(), 1);
+
       std::vector<std::vector<unsigned int>>
         component_to_system_index_face_array(fe.n_components());
 
@@ -373,7 +526,6 @@ namespace internal
 
       std::vector<unsigned int> idx_offset = {0};
 
-
       for (unsigned int base_element_index = 0;
            base_element_index < cell->get_fe().n_base_elements();
            ++base_element_index)
@@ -384,386 +536,252 @@ namespace internal
             idx_offset.back() +
             cell->get_fe().base_element(base_element_index).n_dofs_per_cell());
 
-      for (unsigned int base_element_index = 0, comp = 0;
-           base_element_index < cell->get_fe().n_base_elements();
-           ++base_element_index)
-        for (unsigned int c = 0;
-             c < cell->get_fe().element_multiplicity(base_element_index);
-             ++c, ++comp)
+      std::vector<types::global_dof_index> neighbor_dofs_all(idx_offset.back());
+      std::vector<types::global_dof_index> neighbor_dofs_all_temp(
+        idx_offset.back());
+
+      const auto get_face_idx = [](const auto n_dofs_1d,
+                                   const auto face_no,
+                                   const auto i,
+                                   const auto j) -> unsigned int {
+        const auto direction = face_no / 2;
+        const auto side      = face_no % 2;
+        const auto offset    = (side == 1) ? (n_dofs_1d - 1) : 0;
+
+        if (dim == 2)
+          return (direction == 0) ? (n_dofs_1d * i + offset) :
+                                    (n_dofs_1d * offset + i);
+        else if (dim == 3)
+          switch (direction)
+            {
+              case 0:
+                return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset;
+              case 1:
+                return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i;
+              case 2:
+                return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j;
+              default:
+                Assert(false, ExcNotImplemented());
+            }
+
+        Assert(false, ExcNotImplemented());
+
+        return 0;
+      };
+
+      const std::uint16_t kind =
+        static_cast<std::uint16_t>(refinement_configuration);
+      const std::uint16_t subcell   = (kind >> 0) & 7;
+      const std::uint16_t subcell_x = (subcell >> 0) & 1;
+      const std::uint16_t subcell_y = (subcell >> 1) & 1;
+      const std::uint16_t subcell_z = (subcell >> 2) & 1;
+      const std::uint16_t face      = (kind >> 3) & 7;
+      const std::uint16_t edge      = (kind >> 6) & 7;
+
+      for (unsigned int direction = 0; direction < dim; ++direction)
+        if ((face >> direction) & 1U)
           {
-            auto &mask = masks[comp];
-            mask       = ConstraintKinds::unconstrained;
+            const auto side    = ((subcell >> direction) & 1U) == 0;
+            const auto face_no = direction * 2 + side;
+
+            // read DoFs of parent of face, ...
+            cell->neighbor(face_no)
+              ->face(cell->neighbor_face_no(face_no))
+              ->get_dof_indices(neighbor_dofs_all,
+                                cell->neighbor(face_no)->active_fe_index());
+
+            // ... convert the global DoFs to serial ones, and ...
+            if (partitioner)
+              for (auto &index : neighbor_dofs_all)
+                index = partitioner->global_to_local(index);
+
+            for (unsigned int base_element_index = 0, comp = 0;
+                 base_element_index < cell->get_fe().n_base_elements();
+                 ++base_element_index)
+              for (unsigned int c = 0;
+                   c < cell->get_fe().element_multiplicity(base_element_index);
+                   ++c, ++comp)
+                {
+                  if (supported_components[cell->active_fe_index()][comp] ==
+                      false)
+                    continue;
+
+                  const unsigned int n_dofs_1d =
+                    cell->get_fe()
+                      .base_element(base_element_index)
+                      .tensor_degree() +
+                    1;
+                  const unsigned int dofs_per_face =
+                    Utilities::pow(n_dofs_1d, dim - 1);
+                  std::vector<types::global_dof_index> neighbor_dofs(
+                    dofs_per_face);
+                  const auto lex_face_mapping =
+                    FETools::lexicographic_to_hierarchic_numbering<dim - 1>(
+                      n_dofs_1d - 1);
+
+                  // ... extract the DoFs of the current component
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    neighbor_dofs[i] = neighbor_dofs_all
+                      [component_to_system_index_face_array[comp][i]];
+
+                  // fix DoFs depending on orientation, flip, and rotation
+                  if (dim == 2)
+                    {
+                      // TODO: for mixed meshes we need to take care of
+                      // orientation here
+                      Assert(cell->face_orientation(face_no),
+                             ExcNotImplemented());
+                    }
+                  else if (dim == 3)
+                    {
+                      int rotate = 0;                   // TODO
+                      if (cell->face_rotation(face_no)) //
+                        rotate -= 1;                    //
+                      if (cell->face_flip(face_no))     //
+                        rotate -= 2;                    //
+
+                      rotate_face(rotate, n_dofs_1d, neighbor_dofs);
+
+                      if (cell->face_orientation(face_no) == false)
+                        transpose_face(n_dofs_1d - 1, neighbor_dofs);
+                    }
+                  else
+                    {
+                      Assert(false, ExcNotImplemented());
+                    }
+
+                  // update DoF map
+                  for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i)
+                    for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d);
+                         ++j, ++k)
+                      dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) +
+                                  idx_offset[comp]] =
+                        neighbor_dofs[lex_face_mapping[k]];
+                }
+          }
 
-            const auto &fe_base =
-              cell->get_fe().base_element(base_element_index);
+      if (dim == 3)
+        for (unsigned int direction = 0; direction < dim; ++direction)
+          if ((edge >> direction) & 1U)
+            {
+              const unsigned int line_no =
+                direction == 0 ?
+                  (local_lines[0][subcell_y][subcell_z]) :
+                  (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) :
+                                    (local_lines[2][subcell_x][subcell_y]));
+
+              const unsigned int line_index = cell->line(line_no)->index();
+
+              const auto edge_neighbor =
+                std::find_if(line_to_cells[line_index].begin(),
+                             line_to_cells[line_index].end(),
+                             [&cell](const auto &edge_neighbor) {
+                               return edge_neighbor.first->is_artificial() ==
+                                        false &&
+                                      edge_neighbor.first->level() <
+                                        cell->level();
+                             });
+
+              if (edge_neighbor == line_to_cells[line_index].end())
+                continue;
+
+              const auto neighbor_cell       = edge_neighbor->first;
+              const auto local_line_neighbor = edge_neighbor->second;
+
+              DoFCellAccessor<dim, dim, false>(
+                &neighbor_cell->get_triangulation(),
+                neighbor_cell->level(),
+                neighbor_cell->index(),
+                &cell->get_dof_handler())
+                .get_dof_indices(neighbor_dofs_all);
+
+              if (partitioner)
+                for (auto &index : neighbor_dofs_all)
+                  index = partitioner->global_to_local(index);
+
+              for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i)
+                neighbor_dofs_all_temp[i] = neighbor_dofs_all
+                  [lexicographic_mapping[cell->active_fe_index()][i]];
+
+              const bool flipped =
+                cell->line_orientation(line_no) !=
+                neighbor_cell->line_orientation(local_line_neighbor);
+
+              for (unsigned int base_element_index = 0, comp = 0;
+                   base_element_index < cell->get_fe().n_base_elements();
+                   ++base_element_index)
+                for (unsigned int c = 0;
+                     c <
+                     cell->get_fe().element_multiplicity(base_element_index);
+                     ++c, ++comp)
+                  {
+                    if (supported_components[cell->active_fe_index()][comp] ==
+                        false)
+                      continue;
 
-            if (dim == 1 ||
-                dynamic_cast<const FE_Q<dim> *>(&fe_base) == nullptr)
-              continue;
+                    const unsigned int n_dofs_1d =
+                      cell->get_fe()
+                        .base_element(base_element_index)
+                        .tensor_degree() +
+                      1;
+
+                    for (unsigned int i = 0; i < n_dofs_1d; ++i)
+                      dof_indices[line_dof_idx(line_no, i, n_dofs_1d) +
+                                  idx_offset[comp]] = neighbor_dofs_all_temp
+                        [line_dof_idx(local_line_neighbor,
+                                      flipped ? (n_dofs_1d - 1 - i) : i,
+                                      n_dofs_1d) +
+                         idx_offset[comp]];
+                  }
+            }
+    }
 
-            const unsigned int fe_degree = fe_base.tensor_degree();
-            const unsigned int n_dofs_1d = fe_degree + 1;
-            const unsigned int dofs_per_face =
-              Utilities::fixed_power<dim - 1>(n_dofs_1d);
 
-            std::vector<types::global_dof_index> neighbor_dofs_all(
-              idx_offset.back());
-            std::vector<types::global_dof_index> neighbor_dofs_all_temp(
-              idx_offset.back());
 
-            std::vector<types::global_dof_index> neighbor_dofs(dofs_per_face);
+    template <int dim>
+    template <typename CellIterator>
+    inline bool
+    HangingNodes<dim>::setup_constraints(
+      const CellIterator &                                      cell,
+      const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
+      const std::vector<std::vector<unsigned int>> &lexicographic_mapping,
+      std::vector<types::global_dof_index> &        dof_indices,
+      const ArrayView<ConstraintKinds> &            masks) const
+    {
+      // 1) check if finite elements support fast hanging-node algorithm
+      const auto supported_components = compute_supported_components(
+        cell->get_dof_handler().get_fe_collection());
+
+      if ([](const auto &supported_components) {
+            return std::none_of(supported_components.begin(),
+                                supported_components.end(),
+                                [](const auto &a) {
+                                  return *std::max_element(a.begin(), a.end());
+                                });
+          }(supported_components))
+        return false;
 
-            const auto lex_face_mapping =
-              FETools::lexicographic_to_hierarchic_numbering<dim - 1>(
-                fe_degree);
+      // 2) determine the refinement configuration of the cell
+      const auto refinement_configuration =
+        compute_refinement_configuration(cell);
 
-            for (const unsigned int face : GeometryInfo<dim>::face_indices())
-              {
-                if ((!cell->at_boundary(face)) &&
-                    (cell->neighbor(face)->has_children() == false))
-                  {
-                    const auto &neighbor = cell->neighbor(face);
+      if (refinement_configuration == ConstraintKinds::unconstrained)
+        return false;
 
-                    if (neighbor->is_artificial())
-                      continue;
+      // 3) update DoF indices of cell for specified components
+      update_dof_indices(cell,
+                         partitioner,
+                         lexicographic_mapping,
+                         supported_components,
+                         refinement_configuration,
+                         dof_indices);
 
-                    // Neighbor is coarser than us, i.e., face is constrained
-                    if (neighbor->level() < cell->level())
-                      {
-                        const unsigned int neighbor_face =
-                          cell->neighbor_face_no(face);
-
-                        // Find position of face on neighbor
-                        unsigned int subface = 0;
-                        for (;
-                             subface < GeometryInfo<dim>::max_children_per_face;
-                             ++subface)
-                          if (neighbor->neighbor_child_on_subface(neighbor_face,
-                                                                  subface) ==
-                              cell)
-                            break;
-
-                        // Get indices to read
-                        DoFAccessor<dim - 1, dim, dim, false>(
-                          &neighbor->face(neighbor_face)->get_triangulation(),
-                          neighbor->face(neighbor_face)->level(),
-                          neighbor->face(neighbor_face)->index(),
-                          &cell->get_dof_handler())
-                          .get_dof_indices(neighbor_dofs_all);
-
-                        for (unsigned int i = 0; i < dofs_per_face; ++i)
-                          neighbor_dofs[i] = neighbor_dofs_all
-                            [component_to_system_index_face_array[comp][i]];
-
-                        // If the vector is distributed, we need to transform
-                        // the global indices to local ones.
-                        if (partitioner)
-                          for (auto &index : neighbor_dofs)
-                            index = partitioner->global_to_local(index);
-
-                        if (dim == 2)
-                          {
-                            if (face < 2)
-                              {
-                                mask |= ConstraintKinds::face_x;
-                                if (face == 0)
-                                  mask |= ConstraintKinds::subcell_x;
-                                if (subface == 0)
-                                  mask |= ConstraintKinds::subcell_y;
-                              }
-                            else
-                              {
-                                mask |= ConstraintKinds::face_y;
-                                if (face == 2)
-                                  mask |= ConstraintKinds::subcell_y;
-                                if (subface == 0)
-                                  mask |= ConstraintKinds::subcell_x;
-                              }
-
-                            // Reorder neighbor_dofs and copy into faceth face
-                            // of dof_indices
-
-                            // Offset if upper/right face
-                            unsigned int offset =
-                              (face % 2 == 1) ? fe_degree : 0;
-
-                            for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                              {
-                                unsigned int idx = 0;
-                                // If X-line, i.e., if y = 0 or y = fe_degree
-                                if (face > 1)
-                                  idx = n_dofs_1d * offset + i;
-                                // If Y-line, i.e., if x = 0 or x = fe_degree
-                                else
-                                  idx = n_dofs_1d * i + offset;
-
-                                dof_indices[idx + idx_offset[comp]] =
-                                  neighbor_dofs[lex_face_mapping[i]];
-                              }
-                          }
-                        else if (dim == 3)
-                          {
-                            const bool transpose =
-                              !(cell->face_orientation(face));
-
-                            int rotate = 0;
-
-                            if (cell->face_rotation(face))
-                              rotate -= 1;
-                            if (cell->face_flip(face))
-                              rotate -= 2;
-
-                            rotate_face(rotate, n_dofs_1d, neighbor_dofs);
-                            rotate_subface_index(rotate, subface);
-
-                            if (transpose)
-                              {
-                                transpose_face(fe_degree, neighbor_dofs);
-                                transpose_subface_index(subface);
-                              }
-
-                            // YZ-plane
-                            if (face < 2)
-                              {
-                                mask |= ConstraintKinds::face_x;
-                                if (face == 0)
-                                  mask |= ConstraintKinds::subcell_x;
-                                if (subface % 2 == 0)
-                                  mask |= ConstraintKinds::subcell_y;
-                                if (subface / 2 == 0)
-                                  mask |= ConstraintKinds::subcell_z;
-                              }
-                            // XZ-plane
-                            else if (face < 4)
-                              {
-                                mask |= ConstraintKinds::face_y;
-                                if (face == 2)
-                                  mask |= ConstraintKinds::subcell_y;
-                                if (subface % 2 == 0)
-                                  mask |= ConstraintKinds::subcell_z;
-                                if (subface / 2 == 0)
-                                  mask |= ConstraintKinds::subcell_x;
-                              }
-                            // XY-plane
-                            else
-                              {
-                                mask |= ConstraintKinds::face_z;
-                                if (face == 4)
-                                  mask |= ConstraintKinds::subcell_z;
-                                if (subface % 2 == 0)
-                                  mask |= ConstraintKinds::subcell_x;
-                                if (subface / 2 == 0)
-                                  mask |= ConstraintKinds::subcell_y;
-                              }
-
-                            // Offset if upper/right/back face
-                            unsigned int offset =
-                              (face % 2 == 1) ? fe_degree : 0;
-
-                            for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                              {
-                                for (unsigned int j = 0; j < n_dofs_1d; ++j)
-                                  {
-                                    unsigned int idx = 0;
-                                    // If YZ-plane, i.e., if x = 0 or x =
-                                    // fe_degree, and orientation standard
-                                    if (face < 2)
-                                      idx = n_dofs_1d * n_dofs_1d * i +
-                                            n_dofs_1d * j + offset;
-                                    // If XZ-plane, i.e., if y = 0 or y =
-                                    // fe_degree, and orientation standard
-                                    else if (face < 4)
-                                      idx = n_dofs_1d * n_dofs_1d * j +
-                                            n_dofs_1d * offset + i;
-                                    // If XY-plane, i.e., if z = 0 or z =
-                                    // fe_degree, and orientation standard
-                                    else
-                                      idx = n_dofs_1d * n_dofs_1d * offset +
-                                            n_dofs_1d * i + j;
-
-                                    dof_indices[idx + idx_offset[comp]] =
-                                      neighbor_dofs
-                                        [lex_face_mapping[n_dofs_1d * i + j]];
-                                  }
-                              }
-                          }
-                        else
-                          ExcNotImplemented();
-                      }
-                  }
-              }
-
-            // In 3D we can have a situation where only DoFs on an edge are
-            // constrained. Append these here.
-            if (dim == 3)
-              {
-                // For each line on cell, which faces does it belong to, what is
-                // the edge mask, what is the types of the faces it belong to,
-                // and what is the type along the edge.
-                const ConstraintKinds line_to_edge[12][4] = {
-                  {ConstraintKinds::face_x | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_y,
-                   ConstraintKinds::subcell_x | ConstraintKinds::subcell_z,
-                   ConstraintKinds::subcell_y},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_y,
-                   ConstraintKinds::subcell_z,
-                   ConstraintKinds::subcell_y},
-                  {ConstraintKinds::face_y | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_x,
-                   ConstraintKinds::subcell_y | ConstraintKinds::subcell_z,
-                   ConstraintKinds::subcell_x},
-                  {ConstraintKinds::face_y | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_x,
-                   ConstraintKinds::subcell_z,
-                   ConstraintKinds::subcell_x},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_y,
-                   ConstraintKinds::subcell_x,
-                   ConstraintKinds::subcell_y},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_y,
-                   ConstraintKinds::unconstrained,
-                   ConstraintKinds::subcell_y},
-                  {ConstraintKinds::face_y | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_x,
-                   ConstraintKinds::subcell_y,
-                   ConstraintKinds::subcell_x},
-                  {ConstraintKinds::face_y | ConstraintKinds::face_z,
-                   ConstraintKinds::edge_x,
-                   ConstraintKinds::unconstrained,
-                   ConstraintKinds::subcell_x},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_y,
-                   ConstraintKinds::edge_z,
-                   ConstraintKinds::subcell_x | ConstraintKinds::subcell_y,
-                   ConstraintKinds::subcell_z},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_y,
-                   ConstraintKinds::edge_z,
-                   ConstraintKinds::subcell_y,
-                   ConstraintKinds::subcell_z},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_y,
-                   ConstraintKinds::edge_z,
-                   ConstraintKinds::subcell_x,
-                   ConstraintKinds::subcell_z},
-                  {ConstraintKinds::face_x | ConstraintKinds::face_y,
-                   ConstraintKinds::edge_z,
-                   ConstraintKinds::unconstrained,
-                   ConstraintKinds::subcell_z}};
-
-                for (unsigned int local_line = 0;
-                     local_line < GeometryInfo<dim>::lines_per_cell;
-                     ++local_line)
-                  {
-                    // If we don't already have a constraint for as part of a
-                    // face
-                    if ((mask & line_to_edge[local_line][0]) ==
-                        ConstraintKinds::unconstrained)
-                      {
-                        // For each cell which share that edge
-                        const unsigned int line =
-                          cell->line(local_line)->index();
-                        for (const auto &edge_neighbor : line_to_cells[line])
-                          {
-                            // If one of them is coarser than us
-                            const auto neighbor_cell = edge_neighbor.first;
-
-                            if (neighbor_cell->is_artificial())
-                              continue;
-
-                            if (neighbor_cell->level() < cell->level())
-                              {
-                                const unsigned int local_line_neighbor =
-                                  edge_neighbor.second;
-                                mask |= line_to_edge[local_line][1] |
-                                        line_to_edge[local_line][2];
-
-                                bool flipped = false;
-                                if (cell->line(local_line)->vertex_index(0) ==
-                                    neighbor_cell->line(local_line_neighbor)
-                                      ->vertex_index(0))
-                                  {
-                                    // Assuming line directions match axes
-                                    // directions, we have an unflipped edge of
-                                    // first type
-                                    mask |= line_to_edge[local_line][3];
-                                  }
-                                else if (cell->line(local_line)
-                                           ->vertex_index(1) ==
-                                         neighbor_cell
-                                           ->line(local_line_neighbor)
-                                           ->vertex_index(1))
-                                  {
-                                    // We have an unflipped edge of second type
-                                  }
-                                else if (cell->line(local_line)
-                                           ->vertex_index(1) ==
-                                         neighbor_cell
-                                           ->line(local_line_neighbor)
-                                           ->vertex_index(0))
-                                  {
-                                    // We have a flipped edge of second type
-                                    flipped = true;
-                                  }
-                                else if (cell->line(local_line)
-                                           ->vertex_index(0) ==
-                                         neighbor_cell
-                                           ->line(local_line_neighbor)
-                                           ->vertex_index(1))
-                                  {
-                                    // We have a flipped edge of first type
-                                    mask |= line_to_edge[local_line][3];
-                                    flipped = true;
-                                  }
-                                else
-                                  ExcInternalError();
-
-                                // Copy the unconstrained values
-                                DoFCellAccessor<dim, dim, false>(
-                                  &neighbor_cell->get_triangulation(),
-                                  neighbor_cell->level(),
-                                  neighbor_cell->index(),
-                                  &cell->get_dof_handler())
-                                  .get_dof_indices(neighbor_dofs_all);
-                                // If the vector is distributed, we need to
-                                // transform the global indices to local ones.
-                                if (partitioner)
-                                  for (auto &index : neighbor_dofs_all)
-                                    index = partitioner->global_to_local(index);
-
-                                for (unsigned int i = 0;
-                                     i < neighbor_dofs_all_temp.size();
-                                     ++i)
-                                  neighbor_dofs_all_temp[i] =
-                                    neighbor_dofs_all[lexicographic_mapping[i]];
-
-                                for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                                  {
-                                    // Get local dof index along line
-                                    const unsigned int idx =
-                                      line_dof_idx(local_line, i, n_dofs_1d);
-
-                                    dof_indices[idx + idx_offset[comp]] =
-                                      neighbor_dofs_all_temp
-                                        [line_dof_idx(local_line_neighbor,
-                                                      flipped ? fe_degree - i :
-                                                                i,
-                                                      n_dofs_1d) +
-                                         idx_offset[comp]];
-                                  }
-
-                                // Stop looping over edge neighbors
-                                break;
-                              }
-                          }
-                      }
-                  }
-              }
-            Assert(check(mask, dim), ExcInternalError());
+      // 4)  TODO: copy refinement configuration to all components
+      for (unsigned int c = 0; c < supported_components[0].size(); ++c)
+        if (supported_components[cell->active_fe_index()][c])
+          masks[c] = refinement_configuration;
 
-            cell_has_hanging_node_constraints |=
-              mask != ConstraintKinds::unconstrained;
-          }
-      return cell_has_hanging_node_constraints;
+      return true;
     }
 
 
index 3a7dc80fc15d9e52d234262f6a2905a33ccfb9b1..7b5c1911ebd1e65547837e6b54e64a89d3fe2e8a 100644 (file)
@@ -1023,7 +1023,7 @@ namespace internal
     MatrixFreeFunctions::FaceSetup<dim> &               face_setup,
     MatrixFreeFunctions::ConstraintValues<double> &     constraint_values,
     const bool use_vector_data_exchanger_full,
-    const bool use_fast_hanging_node_algorithm)
+    const bool use_fast_hanging_node_algorithm_in)
   {
     if (do_face_integrals)
       face_setup.initialize(dof_handler[0]->get_triangulation(),
@@ -1049,6 +1049,34 @@ namespace internal
 
     bool cell_categorization_enabled = !cell_vectorization_category.empty();
 
+    bool use_fast_hanging_node_algorithm = use_fast_hanging_node_algorithm_in;
+
+    if (use_fast_hanging_node_algorithm)
+      {
+        const auto &reference_cells =
+          dof_handler[0]->get_triangulation().get_reference_cells();
+        use_fast_hanging_node_algorithm =
+          std::all_of(reference_cells.begin(),
+                      reference_cells.end(),
+                      [](const auto &r) {
+                        return r.is_hyper_cube() || r.is_simplex();
+                      });
+      }
+
+    if (use_fast_hanging_node_algorithm)
+      for (unsigned int no = 0; no < n_dof_handlers; ++no)
+        {
+          const dealii::hp::FECollection<dim> &fes =
+            dof_handler[no]->get_fe_collection();
+
+          use_fast_hanging_node_algorithm &=
+            std::all_of(fes.begin(), fes.end(), [&fes](const auto &fe) {
+              return fes[0].compare_for_domination(fe) ==
+                     FiniteElementDomination::Domination::
+                       either_element_can_dominate;
+            });
+        }
+
     for (unsigned int no = 0; no < n_dof_handlers; ++no)
       {
         const dealii::hp::FECollection<dim> &fes =
@@ -1171,9 +1199,8 @@ namespace internal
         hanging_nodes = std::make_unique<
           dealii::internal::MatrixFreeFunctions::HangingNodes<dim>>(tria);
         for (unsigned int no = 0; no < n_dof_handlers; ++no)
-          if (dof_handler[no]->get_fe_collection().size() == 1)
-            dof_info[no].hanging_node_constraint_masks.resize(
-              n_active_cells * dof_handler[no]->get_fe().n_components());
+          dof_info[no].hanging_node_constraint_masks.resize(
+            n_active_cells * dof_handler[no]->get_fe().n_components());
       }
 
     for (unsigned int counter = 0; counter < n_active_cells; ++counter)
@@ -1212,15 +1239,14 @@ namespace internal
 
                 bool cell_has_hanging_node_constraints = false;
 
-                if (dim > 1 && use_fast_hanging_node_algorithm &&
-                    dofh->get_fe_collection().size() == 1)
+                if (dim > 1 && use_fast_hanging_node_algorithm)
                   {
                     local_dof_indices_resolved = local_dof_indices;
 
                     cell_has_hanging_node_constraints =
                       dof_info[no].process_hanging_node_constraints(
                         *hanging_nodes,
-                        lexicographic[no][0],
+                        lexicographic[no],
                         counter,
                         cell_it,
                         local_dof_indices_resolved);
index 8e0eaf5a7383ac1eb60422b633164f54e7cff844..6bd7669f1149892bda8ffa202e7a65ad765d2946 100644 (file)
@@ -1519,21 +1519,21 @@ namespace internal
     template bool
     DoFInfo::process_hanging_node_constraints<1>(
       const HangingNodes<1> &                           hanging_nodes,
-      const std::vector<unsigned int> &                 lexicographic_mapping,
+      const std::vector<std::vector<unsigned int>> &    lexicographic_mapping,
       const unsigned int                                cell_number,
       const TriaIterator<DoFCellAccessor<1, 1, false>> &cell,
       std::vector<types::global_dof_index> &            dof_indices);
     template bool
     DoFInfo::process_hanging_node_constraints<2>(
       const HangingNodes<2> &                           hanging_nodes,
-      const std::vector<unsigned int> &                 lexicographic_mapping,
+      const std::vector<std::vector<unsigned int>> &    lexicographic_mapping,
       const unsigned int                                cell_number,
       const TriaIterator<DoFCellAccessor<2, 2, false>> &cell,
       std::vector<types::global_dof_index> &            dof_indices);
     template bool
     DoFInfo::process_hanging_node_constraints<3>(
       const HangingNodes<3> &                           hanging_nodes,
-      const std::vector<unsigned int> &                 lexicographic_mapping,
+      const std::vector<std::vector<unsigned int>> &    lexicographic_mapping,
       const unsigned int                                cell_number,
       const TriaIterator<DoFCellAccessor<3, 3, false>> &cell,
       std::vector<types::global_dof_index> &            dof_indices);
diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc b/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc
new file mode 100644 (file)
index 0000000..6190362
--- /dev/null
@@ -0,0 +1,165 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test that the fast matrix-free hanging-node algorithm is also working on
+// adaptively refined 2D mixed meshes.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+
+
+int
+main()
+{
+  initlog();
+
+  const unsigned int dim    = 2;
+  const unsigned int degree = 1;
+
+  Triangulation<dim> tria_0, tria_1, tria_2, tria_3, tria;
+
+  GridGenerator::subdivided_hyper_rectangle(tria_0,
+                                            {1, 1},
+                                            {0.0, 0.0},
+                                            {0.5, 0.5});
+  GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_1,
+                                                           {1, 1},
+                                                           {0.5, 0.0},
+                                                           {1.0, 0.5});
+  GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_2,
+                                                           {1, 1},
+                                                           {0.0, 0.5},
+                                                           {0.5, 1.0});
+  GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_3,
+                                                           {1, 1},
+                                                           {0.5, 0.5},
+                                                           {1.0, 1.0});
+
+  GridGenerator::merge_triangulations({&tria_0, &tria_1, &tria_2, &tria_3},
+                                      tria);
+
+  auto cell = tria.begin();
+
+  cell->set_refine_flag();
+  cell++;
+  cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  if (false)
+    {
+      GridOut       grid_out;
+      std::ofstream out("mesh.vtk");
+      grid_out.write_vtk(tria, out);
+    }
+
+  DoFHandler<dim> dof_handler(tria);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->reference_cell() == ReferenceCells::Triangle)
+        cell->set_active_fe_index(0);
+      else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
+        cell->set_active_fe_index(1);
+      else
+        Assert(false, ExcNotImplemented());
+    }
+
+  const hp::MappingCollection<2> mapping(MappingFE<2>(FE_SimplexP<2>(1)),
+                                         MappingFE<2>(FE_Q<2>(1)));
+  const hp::FECollection<2>      fe(FE_SimplexP<2>{degree}, FE_Q<2>{degree});
+  const hp::QCollection<2> quadrature_formula(QGaussSimplex<2>(degree + 1),
+                                              QGauss<2>(degree + 1));
+
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+  const auto print = [](const auto &label, const auto &matrix_free) {
+    deallog << label << std::endl;
+
+    for (unsigned int c = 0;
+         c < matrix_free.get_dof_info(0).row_starts.size() - 1;
+         ++c)
+      {
+        deallog
+          << std::setw(3)
+          << (matrix_free.get_dof_info(0)
+                    .hanging_node_constraint_masks.size() == 0 ?
+                0 :
+                static_cast<unsigned int>(
+                  matrix_free.get_dof_info(0).hanging_node_constraint_masks[c]))
+          << " : ";
+
+        for (unsigned int i = matrix_free.get_dof_info(0).row_starts[c].first;
+             i < matrix_free.get_dof_info(0).row_starts[c + 1].first;
+             ++i)
+          deallog << std::setw(3) << matrix_free.get_dof_info(0).dof_indices[i]
+                  << " ";
+        deallog << std::endl;
+      }
+    deallog << std::endl;
+  };
+
+  {
+    typename MatrixFree<dim, double, VectorizedArray<double, 1>>::AdditionalData
+      additional_data;
+    additional_data.mapping_update_flags = update_gradients | update_values;
+
+    MatrixFree<dim, double, VectorizedArray<double, 1>> matrix_free;
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quadrature_formula, additional_data);
+
+    print("use_fast_hanging_node_algorithm = true", matrix_free);
+  }
+
+  {
+    typename MatrixFree<dim, double, VectorizedArray<double, 1>>::AdditionalData
+      additional_data;
+    additional_data.mapping_update_flags = update_gradients | update_values;
+    additional_data.use_fast_hanging_node_algorithm = false;
+
+    MatrixFree<dim, double, VectorizedArray<double, 1>> matrix_free;
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quadrature_formula, additional_data);
+
+    for (const auto &i :
+         matrix_free.get_dof_info(0).hanging_node_constraint_masks)
+      deallog << static_cast<unsigned int>(i) << std::endl;
+    deallog << std::endl;
+
+    print("use_fast_hanging_node_algorithm = false", matrix_free);
+  }
+}
diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.output b/tests/matrix_free/mixed_mesh_hn_algorithm_01.output
new file mode 100644 (file)
index 0000000..384befa
--- /dev/null
@@ -0,0 +1,32 @@
+
+DEAL::use_fast_hanging_node_algorithm = true
+DEAL::0   : 12  13  14  
+DEAL::0   : 13  2   2   1   
+DEAL::0   : 14  2   1   1   
+DEAL::0   : 13  2   1   14  
+DEAL::0   : 0   1   2   
+DEAL::0   : 3   1   4   
+DEAL::0   : 5   4   1   
+DEAL::0   : 1   0   5   
+DEAL::0   : 6   5   0   
+DEAL::0   : 7   8   9   10  
+DEAL::0   : 8   12  10  14  
+DEAL::17  : 9   10  3   1   
+DEAL::16  : 10  14  3   1   
+DEAL::
+DEAL::
+DEAL::use_fast_hanging_node_algorithm = false
+DEAL::0   : 12  13  14  
+DEAL::0   : 13  2   2   1   
+DEAL::0   : 14  2   1   1   
+DEAL::0   : 13  2   1   14  
+DEAL::0   : 0   1   2   
+DEAL::0   : 3   1   4   
+DEAL::0   : 5   4   1   
+DEAL::0   : 1   0   5   
+DEAL::0   : 6   5   0   
+DEAL::0   : 7   8   9   10  
+DEAL::0   : 8   12  10  14  
+DEAL::0   : 9   10  3   3   1   
+DEAL::0   : 10  14  3   1   1   
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.