GridTools::compute_active_cell_halo_layer() for periodic meshes

diff --git a/doc/news/changes/minor/20240212Munch b/doc/news/changes/minor/20240212Munch
new file mode 100644 (file)
index 0000000..4f8f066
--- /dev/null
+++ b/doc/news/changes/minor/20240212Munch
@@ -0,0 +1,4 @@
+Fixed: The function GridTools::compute_active_cell_halo_layer() now
+also supports perodic meshes.
+<br>
+(Peter Munch, 2024/02/12)

diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc
index aeeacbbc260090e6a54cc738e579e5855048c3cf..802af7c112c7f4f1955b57a2925e2750bc00d546 100644 (file)
--- a/source/dofs/dof_tools_constraints.cc
+++ b/source/dofs/dof_tools_constraints.cc
@@ -3240,6 +3240,22 @@ namespace DoFTools
             return;
           }
 
+      // In the case of shared Trinagulation with artificial cells all
+      // cells have valid DoF indices, i.e., the check above does not work.
+      if (const auto tria = dynamic_cast<
+            const parallel::shared::Triangulation<dim, spacedim> *>(
+            &face_1->get_triangulation()))
+        if (tria->with_artificial_cells() &&
+            (affine_constraints.get_local_lines().size() != 0))
+          for (unsigned int i = 0; i < dofs_per_face; ++i)
+            if ((affine_constraints.get_local_lines().is_element(dofs_1[i]) ==
+                 false) ||
+                (affine_constraints.get_local_lines().is_element(dofs_2[i]) ==
+                 false))
+              {
+                return;
+              }
+
       // Well, this is a hack:
       //
       // There is no

diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
index 46f47c93af12c6ac44274c446ef7ccbac5db3e7e..76b825ef4f65d25cb338cff14f925b81058e5bf7 100644 (file)
--- a/source/grid/grid_tools_dof_handlers.cc
+++ b/source/grid/grid_tools_dof_handlers.cc
@@ -824,13 +824,24 @@ namespace GridTools
     std::vector<bool> locally_active_vertices_on_subdomain(
       mesh.get_triangulation().n_vertices(), false);
 
+    std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
+    std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
+    GridTools::collect_coinciding_vertices(mesh.get_triangulation(),
+                                           coinciding_vertex_groups,
+                                           vertex_to_coinciding_vertex_group);
+
     // Find the cells for which the predicate is true
     // These are the cells around which we wish to construct
     // the halo layer
     for (const auto &cell : mesh.active_cell_iterators())
       if (predicate(cell)) // True predicate --> Part of subdomain
         for (const auto v : cell->vertex_indices())
-          locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;
+          {
+            locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;
+            for (const auto vv : coinciding_vertex_groups
+                   [vertex_to_coinciding_vertex_group[cell->vertex_index(v)]])
+              locally_active_vertices_on_subdomain[vv] = true;
+          }
 
     // Find the cells that do not conform to the predicate
     // but share a vertex with the selected subdomain

diff --git a/tests/sharedtria/pbc_01.cc b/tests/sharedtria/pbc_01.cc
new file mode 100644 (file)
index 0000000..cf6958e
--- /dev/null
+++ b/tests/sharedtria/pbc_01.cc
@@ -0,0 +1,75 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2024 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 shared triangulations (with artificial cells) for periodic
+// boundary conditions.
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                            Triangulation<dim>::none,
+                                            true);
+
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+
+  std::vector<
+    GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+    face_pairs;
+  for (unsigned int d = 0; d < dim; ++d)
+    GridTools::collect_periodic_faces(tria, 2 * d, 2 * d + 1, d, face_pairs);
+  tria.add_periodicity(face_pairs);
+
+  tria.refine_global(6 - dim);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+  AffineConstraints<double> constraints;
+
+  constraints.reinit(DoFTools::extract_locally_relevant_dofs(dof_handler));
+  for (unsigned int d = 0; d < dim; ++d)
+    DoFTools::make_periodicity_constraints(
+      dof_handler, 2 * d, 2 * d + 1, d, constraints);
+  constraints.close();
+
+  deallog << "OK" << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MPILogInitAll all;
+
+  test<1>();
+  test<2>();
+  test<3>();
+}

diff --git a/tests/sharedtria/pbc_01.mpirun=2.output b/tests/sharedtria/pbc_01.mpirun=2.output
new file mode 100644 (file)
index 0000000..5693ef9
--- /dev/null
+++ b/tests/sharedtria/pbc_01.mpirun=2.output
@@ -0,0 +1,9 @@
+
+DEAL:0::OK
+DEAL:0::OK
+DEAL:0::OK
+
+DEAL:1::OK
+DEAL:1::OK
+DEAL:1::OK
+