Enable periodicity for 1D

diff --git a/doc/news/changes/minor/20190727PeterMunch b/doc/news/changes/minor/20190727PeterMunch
new file mode 100644 (file)
index 0000000..7aa7d1c
--- /dev/null
+++ b/doc/news/changes/minor/20190727PeterMunch
@@ -0,0 +1,3 @@
+New: Enable periodicity for 1D triangulations (only for dealii::DoFHandler).
+<br>
+(Peter Munch, 2019/07/27)

diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h
index f8a4b7a6905c37b1bee1cabe807274158102aa22..be8285f8baad2fe1fd5257d9c6e6470f70b81961 100644 (file)
--- a/include/deal.II/dofs/dof_accessor.templates.h
+++ b/include/deal.II/dofs/dof_accessor.templates.h
@@ -35,6 +35,7 @@
 #include <deal.II/lac/read_write_vector.h>
 
 #include <limits>
+#include <type_traits>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
@@ -2565,8 +2566,11 @@ inline unsigned int
 DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>::
   n_active_fe_indices() const
 {
-  Assert(false, ExcNotImplemented());
-  return numbers::invalid_unsigned_int;
+  Assert((std::is_same<DoFHandlerType<1, spacedim>,
+                       dealii::DoFHandler<1, spacedim>>::value == true),
+         ExcNotImplemented());
+
+  return 1;
 }
 
 
@@ -2578,8 +2582,11 @@ inline unsigned int
 DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>::
   nth_active_fe_index(const unsigned int /*n*/) const
 {
-  Assert(false, ExcNotImplemented());
-  return numbers::invalid_unsigned_int;
+  Assert((std::is_same<DoFHandlerType<1, spacedim>,
+                       dealii::DoFHandler<1, spacedim>>::value == true),
+         ExcNotImplemented());
+
+  return 0;
 }
 
 

diff --git a/source/dofs/dof_tools_constraints.inst.in b/source/dofs/dof_tools_constraints.inst.in
index 685a98d02225014d593821905753ec34a36b0bdc..931ef21f2663074317e738c911b3bf8b94f0d904 100644 (file)
--- a/source/dofs/dof_tools_constraints.inst.in
+++ b/source/dofs/dof_tools_constraints.inst.in
@@ -49,7 +49,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS)
       AffineConstraints<std::complex<double>> &);
 #endif
 
-#if deal_II_dimension != 1
     template void DoFTools::make_periodicity_constraints(
       const DH<deal_II_dimension>::face_iterator &,
       const DH<deal_II_dimension>::face_iterator &,
@@ -83,7 +82,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS)
       const int,
       AffineConstraints<double> &,
       const ComponentMask &);
-#endif
   }
 
 for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS;

diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 31615cb9a66dc66c4005c0404454b939c113aac6..17efe915c0b3322a888f6abaf630cfad85e15603 100644 (file)
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -950,7 +950,12 @@ namespace
     Assert(face_1->at_boundary() && face_2->at_boundary(),
            ExcMessage("Periodic faces must be on the boundary"));
 
-    Assert(std::abs(cell_1->level() - cell_2->level()) < 2, ExcInternalError());
+    // Check if the requirement that each edge can only have at most one hanging
+    // node, and as a consequence neighboring cells can differ by at most
+    // one refinement level is enforced. In 1d, there are no hanging nodes and
+    // so neighboring cells can differ by more than one refinement level.
+    Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
+           ExcInternalError());
 
     // insert periodic face pair for both cells
     using CellFace =
@@ -968,127 +973,24 @@ namespace
     Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
     periodic_face_map.insert(periodic_faces);
 
-    // A lookup table on how to go through the child cells depending on the
-    // orientation:
-    // see Documentation of GeometryInfo for details
-
-    static const int lookup_table_2d[2][2] =
-      //               flip:
+    if (dim == 1)
       {
-        {0, 1}, // false
-        {1, 0}  // true
-      };
-
-    static const int lookup_table_3d[2][2][2][4] =
-      //                           orientation flip  rotation
-      {{{
-          {0, 2, 1, 3}, // false       false false
-          {2, 3, 0, 1}  // false       false true
-        },
-        {
-          {3, 1, 2, 0}, // false       true  false
-          {1, 0, 3, 2}  // false       true  true
-        }},
-       {{
-          {0, 1, 2, 3}, // true        false false
-          {1, 3, 0, 2}  // true        false true
-        },
-        {
-          {3, 2, 1, 0}, // true        true  false
-          {2, 0, 3, 1}  // true        true  true
-        }}};
-
-    if (cell_1->has_children())
-      {
-        if (cell_2->has_children())
+        if (cell_1->has_children())
           {
-            // In the case that both faces have children, we loop over all
-            // children and apply update_periodic_face_map_recursively
-            // recursively:
-
-            Assert(face_1->n_children() ==
-                       GeometryInfo<dim>::max_children_per_face &&
-                     face_2->n_children() ==
-                       GeometryInfo<dim>::max_children_per_face,
-                   ExcNotImplemented());
-
-            for (unsigned int i = 0;
-                 i < GeometryInfo<dim>::max_children_per_face;
-                 ++i)
+            if (cell_2->has_children())
               {
-                // Lookup the index for the second face
-                unsigned int j = 0;
-                switch (dim)
-                  {
-                    case 2:
-                      j = lookup_table_2d[face_flip][i];
-                      break;
-                    case 3:
-                      j = lookup_table_3d[face_orientation][face_flip]
-                                         [face_rotation][i];
-                      break;
-                    default:
-                      AssertThrow(false, ExcNotImplemented());
-                  }
-
-                // find subcell ids that belong to the subface indices
-                unsigned int child_cell_1 =
-                  GeometryInfo<dim>::child_cell_on_face(
-                    cell_1->refinement_case(),
-                    n_face_1,
-                    i,
-                    cell_1->face_orientation(n_face_1),
-                    cell_1->face_flip(n_face_1),
-                    cell_1->face_rotation(n_face_1),
-                    face_1->refinement_case());
-                unsigned int child_cell_2 =
-                  GeometryInfo<dim>::child_cell_on_face(
-                    cell_2->refinement_case(),
-                    n_face_2,
-                    j,
-                    cell_2->face_orientation(n_face_2),
-                    cell_2->face_flip(n_face_2),
-                    cell_2->face_rotation(n_face_2),
-                    face_2->refinement_case());
-
-                Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
-                         face_1->child(i),
-                       ExcInternalError());
-                Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
-                         face_2->child(j),
-                       ExcInternalError());
-
-                // precondition: subcell has the same orientation as cell (so
-                // that the face numbers coincide) recursive call
                 update_periodic_face_map_recursively<dim, spacedim>(
-                  cell_1->child(child_cell_1),
-                  cell_2->child(child_cell_2),
+                  cell_1->child(n_face_1),
+                  cell_2->child(n_face_2),
                   n_face_1,
                   n_face_2,
                   orientation,
                   periodic_face_map);
               }
-          }
-        else // only face_1 has children
-          {
-            for (unsigned int i = 0;
-                 i < GeometryInfo<dim>::max_children_per_face;
-                 ++i)
+            else // only face_1 has children
               {
-                // find subcell ids that belong to the subface indices
-                unsigned int child_cell_1 =
-                  GeometryInfo<dim>::child_cell_on_face(
-                    cell_1->refinement_case(),
-                    n_face_1,
-                    i,
-                    cell_1->face_orientation(n_face_1),
-                    cell_1->face_flip(n_face_1),
-                    cell_1->face_rotation(n_face_1),
-                    face_1->refinement_case());
-
-                // recursive call
                 update_periodic_face_map_recursively<dim, spacedim>(
-                  cell_1->child(child_cell_1),
+                  cell_1->child(n_face_1),
                   cell_2,
                   n_face_1,
                   n_face_2,
@@ -1097,6 +999,138 @@ namespace
               }
           }
       }
+    else // dim == 2 || dim == 3
+      {
+        // A lookup table on how to go through the child cells depending on the
+        // orientation:
+        // see Documentation of GeometryInfo for details
+
+        static const int lookup_table_2d[2][2] =
+          //               flip:
+          {
+            {0, 1}, // false
+            {1, 0}  // true
+          };
+
+        static const int lookup_table_3d[2][2][2][4] =
+          //                           orientation flip  rotation
+          {{{
+              {0, 2, 1, 3}, // false       false false
+              {2, 3, 0, 1}  // false       false true
+            },
+            {
+              {3, 1, 2, 0}, // false       true  false
+              {1, 0, 3, 2}  // false       true  true
+            }},
+           {{
+              {0, 1, 2, 3}, // true        false false
+              {1, 3, 0, 2}  // true        false true
+            },
+            {
+              {3, 2, 1, 0}, // true        true  false
+              {2, 0, 3, 1}  // true        true  true
+            }}};
+
+        if (cell_1->has_children())
+          {
+            if (cell_2->has_children())
+              {
+                // In the case that both faces have children, we loop over all
+                // children and apply update_periodic_face_map_recursively
+                // recursively:
+
+                Assert(face_1->n_children() ==
+                           GeometryInfo<dim>::max_children_per_face &&
+                         face_2->n_children() ==
+                           GeometryInfo<dim>::max_children_per_face,
+                       ExcNotImplemented());
+
+                for (unsigned int i = 0;
+                     i < GeometryInfo<dim>::max_children_per_face;
+                     ++i)
+                  {
+                    // Lookup the index for the second face
+                    unsigned int j = 0;
+                    switch (dim)
+                      {
+                        case 2:
+                          j = lookup_table_2d[face_flip][i];
+                          break;
+                        case 3:
+                          j = lookup_table_3d[face_orientation][face_flip]
+                                             [face_rotation][i];
+                          break;
+                        default:
+                          AssertThrow(false, ExcNotImplemented());
+                      }
+
+                    // find subcell ids that belong to the subface indices
+                    unsigned int child_cell_1 =
+                      GeometryInfo<dim>::child_cell_on_face(
+                        cell_1->refinement_case(),
+                        n_face_1,
+                        i,
+                        cell_1->face_orientation(n_face_1),
+                        cell_1->face_flip(n_face_1),
+                        cell_1->face_rotation(n_face_1),
+                        face_1->refinement_case());
+                    unsigned int child_cell_2 =
+                      GeometryInfo<dim>::child_cell_on_face(
+                        cell_2->refinement_case(),
+                        n_face_2,
+                        j,
+                        cell_2->face_orientation(n_face_2),
+                        cell_2->face_flip(n_face_2),
+                        cell_2->face_rotation(n_face_2),
+                        face_2->refinement_case());
+
+                    Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
+                             face_1->child(i),
+                           ExcInternalError());
+                    Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
+                             face_2->child(j),
+                           ExcInternalError());
+
+                    // precondition: subcell has the same orientation as cell
+                    // (so that the face numbers coincide) recursive call
+                    update_periodic_face_map_recursively<dim, spacedim>(
+                      cell_1->child(child_cell_1),
+                      cell_2->child(child_cell_2),
+                      n_face_1,
+                      n_face_2,
+                      orientation,
+                      periodic_face_map);
+                  }
+              }
+            else // only face_1 has children
+              {
+                for (unsigned int i = 0;
+                     i < GeometryInfo<dim>::max_children_per_face;
+                     ++i)
+                  {
+                    // find subcell ids that belong to the subface indices
+                    unsigned int child_cell_1 =
+                      GeometryInfo<dim>::child_cell_on_face(
+                        cell_1->refinement_case(),
+                        n_face_1,
+                        i,
+                        cell_1->face_orientation(n_face_1),
+                        cell_1->face_flip(n_face_1),
+                        cell_1->face_rotation(n_face_1),
+                        face_1->refinement_case());
+
+                    // recursive call
+                    update_periodic_face_map_recursively<dim, spacedim>(
+                      cell_1->child(child_cell_1),
+                      cell_2,
+                      n_face_1,
+                      n_face_2,
+                      orientation,
+                      periodic_face_map);
+                  }
+              }
+          }
+      }
   }
 
 

diff --git a/tests/bits/periodicity_01.cc b/tests/bits/periodicity_01.cc
index 5af7e592b72c5d421e88ea2000e50d4e12ee69f0..e381160fc046c88169e1065086155c80a4692fb7 100644 (file)
--- a/tests/bits/periodicity_01.cc
+++ b/tests/bits/periodicity_01.cc
@@ -43,11 +43,12 @@ test()
   Triangulation<dim>        triangulation;
   std::vector<unsigned int> repetitions(dim, 1);
   repetitions[0] = 2;
-  GridGenerator::subdivided_hyper_rectangle(triangulation,
-                                            repetitions,
-                                            Point<dim>(),
-                                            (dim == 2 ? Point<dim>(2, 1) :
-                                                        Point<dim>(2, 1, 1)));
+  GridGenerator::subdivided_hyper_rectangle(
+    triangulation,
+    repetitions,
+    Point<dim>(),
+    (dim == 1 ? Point<dim>(2) :
+                (dim == 2 ? Point<dim>(2, 1) : Point<dim>(2, 1, 1))));
 
   FE_Q<dim>       fe(1);
   DoFHandler<dim> dof_handler(triangulation);
@@ -57,6 +58,8 @@ test()
   DoFTools::make_periodicity_constraints(dof_handler.begin(0)->face(0),
                                          (++dof_handler.begin(0))->face(1),
                                          cm);
+
+  deallog << "dim " << std::to_string(dim) << ":" << std::endl;
   cm.print(deallog.get_file_stream());
 }
 
@@ -66,7 +69,7 @@ int
 main()
 {
   initlog();
-
+  test<1>();
   test<2>();
   test<3>();
   return 0;

diff --git a/tests/bits/periodicity_01.output b/tests/bits/periodicity_01.output
index eeb0da966cdb81abab861aea77f28aad4a8e2f48..97d6535c965d9f844d2947eee20273ccdf035a7a 100644 (file)
--- a/tests/bits/periodicity_01.output
+++ b/tests/bits/periodicity_01.output
@@ -1,6 +1,10 @@
 
+DEAL::dim 1:
+    2 0:  1.00000
+DEAL::dim 2:
     4 0:  1.00000
     5 2:  1.00000
+DEAL::dim 3:
     8 0:  1.00000
     9 2:  1.00000
     10 4:  1.00000

diff --git a/tests/grid/periodicity_01.cc b/tests/grid/periodicity_01.cc
index 588e79f458f1d2897dfd5a8be4e7a70ce945d57e..cc5edd7229acd82363d01a5ec9961b754486e733 100644 (file)
--- a/tests/grid/periodicity_01.cc
+++ b/tests/grid/periodicity_01.cc
@@ -67,7 +67,7 @@ test()
   unsigned int neighbor_count = 0;
 
   for (const auto &cell : tria.active_cell_iterators())
-    for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+    for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
       if (cell->has_periodic_neighbor(f))
         ++neighbor_count;
 
@@ -85,6 +85,7 @@ int
 main()
 {
   initlog();
+  test<1>();
   test<2>();
   test<3>();
   return 0;

diff --git a/tests/grid/periodicity_01.output b/tests/grid/periodicity_01.output
index cf04d7bf8763023f9105cf42fc4ee48a95919e5b..aee6245e1e5446301f77cb51c551267aa79f47a7 100644 (file)
--- a/tests/grid/periodicity_01.output
+++ b/tests/grid/periodicity_01.output
@@ -1,3 +1,4 @@
 
+DEAL::Found 2 faces with periodic neighbor
 DEAL::Found 48 faces with periodic neighbor
 DEAL::Found 1280 faces with periodic neighbor