]> https://gitweb.dealii.org/ - dealii.git/commitdiff
simplex: verify hanging nodes on triangular and hybrid meshes. 11676/head
authorMarc Fehling <mafehling.git@gmail.com>
Mon, 1 Feb 2021 01:15:18 +0000 (18:15 -0700)
committerMarc Fehling <marc.fehling@gmx.net>
Thu, 4 Feb 2021 22:12:13 +0000 (15:12 -0700)
15 files changed:
source/grid/grid_tools.cc
tests/simplex/hanging_nodes.h [new file with mode: 0644]
tests/simplex/hanging_nodes_01.cc
tests/simplex/hanging_nodes_01.with_simplex_support=on.output
tests/simplex/hanging_nodes_02.cc
tests/simplex/hanging_nodes_02.with_simplex_support=on.output
tests/simplex/hanging_nodes_03.cc [new file with mode: 0644]
tests/simplex/hanging_nodes_03.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_01.cc [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_02.cc [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_03.cc [new file with mode: 0644]
tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/simplex_grids.h

index 31bf8691d2fd50b60b37446e60f19f7291595a6d..843c490b52881e6c8643e6cb0b2eea3462224807 100644 (file)
@@ -552,7 +552,7 @@ namespace GridTools
     for (const auto &cell : tria.cell_iterators_on_level(0))
       {
         // Save cell data
-        CellData<dim> cell_data;
+        CellData<dim> cell_data(cell->n_vertices());
         for (const unsigned int cell_vertex_n : cell->vertex_indices())
           {
             Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
diff --git a/tests/simplex/hanging_nodes.h b/tests/simplex/hanging_nodes.h
new file mode 100644 (file)
index 0000000..309ac2e
--- /dev/null
@@ -0,0 +1,238 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// miscellaneous helper funcions for hanging_nodes_* tests
+
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+
+// ----- setup -----
+
+/**
+ * On a Triangulation @p tria, refine those cells which are descendants of the
+ * same coarse cell by a specified amount of times.
+ *
+ * The index of an entry in @p n_refinements corresponds to the globally unique
+ * coarse cell id, while the entry itself describes the number of refinements.
+ */
+template <int dim>
+void
+refine(const std::vector<unsigned int> &n_refinements, Triangulation<dim> &tria)
+{
+  AssertDimension(n_refinements.size(), tria.n_cells(0));
+
+  const unsigned int max_level =
+    *std::max_element(n_refinements.begin(), n_refinements.end());
+
+  while (tria.n_levels() <= max_level)
+    {
+      for (const auto &cell : tria.active_cell_iterators())
+        {
+          const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id();
+
+          if (static_cast<unsigned int>(cell->level()) <
+              n_refinements[coarse_cell_id])
+            cell->set_refine_flag();
+        }
+
+      tria.execute_coarsening_and_refinement();
+    }
+
+#ifdef DEBUG
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id();
+
+      AssertDimension(cell->level(), n_refinements[coarse_cell_id]);
+    }
+#endif
+}
+
+
+/**
+ * On a DoFHandler @p dofh, assign the same active FE index to all cells which
+ * are descendants of the same coarse cell.
+ *
+ * The index of an entry in @p fe_indices corresponds to the globally unique
+ * coarse cell id, while the entry itself is the active FE index.
+ */
+template <int dim>
+void
+set_active_fe_indices(const std::vector<unsigned int> &fe_indices,
+                      DoFHandler<dim> &                dofh)
+{
+  AssertDimension(fe_indices.size(), dofh.get_triangulation().n_cells(0));
+
+  for (const auto &cell : dofh.active_cell_iterators())
+    {
+      const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id();
+
+      cell->set_active_fe_index(fe_indices[coarse_cell_id]);
+    }
+}
+
+
+
+// ----- diagnostics -----
+
+/**
+ * Print the indices of degrees of freedom for each face on each cell in
+ * deallog.
+ *
+ * If a face is divided into subfaces on a neighboring cell which is refined,
+ * then also print the dof indices on the subfaces from the perspective of the
+ * neighboring cell.
+ */
+template <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+  std::vector<types::global_dof_index> dof_indices;
+
+  for (const auto &cell : dofh.active_cell_iterators())
+    {
+      const auto &       fe       = cell->get_fe();
+      const unsigned int fe_index = cell->active_fe_index();
+
+      for (unsigned int f = 0; f < cell->n_faces(); ++f)
+        {
+          const auto &face = cell->face(f);
+          Assert(face->fe_index_is_active(fe_index), ExcInternalError());
+
+          const unsigned int n_dofs =
+            internal::DoFAccessorImplementation::Implementation::n_dof_indices(
+              *face, fe_index);
+          dof_indices.resize(n_dofs);
+          face->get_dof_indices(dof_indices, fe_index);
+
+          deallog << "cell:" << cell->active_cell_index() << " face:" << f
+                  << " dofs:";
+          for (const auto &i : dof_indices)
+            deallog << i << " ";
+          deallog << std::endl;
+
+          // also print dofs on subfaces if neighboring cell is refined.
+          // in this case, only one fe should be active on the subface.
+          for (unsigned int sf = 0; sf < face->n_children(); ++sf)
+            {
+              const auto &       subface = face->child(sf);
+              const unsigned int subface_fe_index =
+                subface->nth_active_fe_index(0);
+              Assert(subface->n_active_fe_indices() == 1, ExcInternalError());
+
+              const unsigned int n_dofs = internal::DoFAccessorImplementation::
+                Implementation::n_dof_indices(*subface, subface_fe_index);
+              dof_indices.resize(n_dofs);
+              subface->get_dof_indices(dof_indices, subface_fe_index);
+
+              deallog << "    subface:" << sf << " dofs:";
+              for (const auto &i : dof_indices)
+                deallog << i << " ";
+              deallog << std::endl;
+            }
+        }
+    }
+}
+
+
+/**
+ * Print the index and physical coordinate of each degree of freedom in deallog.
+ */
+template <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+  hp::MappingCollection<dim> mapping;
+  for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
+    mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
+
+  std::vector<Point<dim>> points(dofh.n_dofs());
+  DoFTools::map_dofs_to_support_points(mapping, dofh, points);
+
+  for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
+    deallog << "dof:" << i << " point:" << points[i] << std::endl;
+}
+
+
+
+// ----- tests -----
+
+/**
+ * Verify hanging node constraints on locally refined meshes.
+ *
+ * Creates a Triangulation based on @p grid_generator, refines the coarse cells
+ * @p n_refinements times, and assigns @p fe_indices to all descendants of the
+ * coarse cells based on the provided @p fe_collection.
+ *
+ * Makes hanging node constraints and prints them in deallog.
+ */
+template <int dim>
+void
+test(const std::vector<unsigned int> &                n_refinements,
+     const std::vector<unsigned int> &                fe_indices,
+     const hp::FECollection<dim> &                    fe_collection,
+     const std::function<void(Triangulation<dim> &)> &grid_generator)
+{
+  // setup grid
+  Triangulation<dim> tria;
+  grid_generator(tria);
+
+  AssertDimension(n_refinements.size(), tria.n_cells());
+  AssertDimension(fe_indices.size(), tria.n_cells());
+
+  refine(n_refinements, tria);
+
+#if false
+  GridOut grid_out;
+  grid_out.write_vtk(tria, deallog.get_file_stream());
+#endif
+
+  DoFHandler<dim> dofh(tria);
+  set_active_fe_indices(fe_indices, dofh);
+  dofh.distribute_dofs(fe_collection);
+  deallog << "ndofs: " << dofh.n_dofs() << std::endl;
+
+#if false
+  print_dof_points(dofh);
+  print_dof_indices_on_faces(dofh);
+#endif
+
+  // hanging node constraints
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dofh, constraints);
+  constraints.close();
+  constraints.print(deallog.get_file_stream());
+
+  deallog << "OK" << std::endl;
+}
index 20762147c64157cdd1cdf387292e82108315d643..ac96c5eb37bd8c6ba78e226dcd5efc8e80b2cdc1 100644 (file)
@@ -31,8 +31,7 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.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/grid/tria.h>
 
 
 #include "../tests.h"
 
-
-// ----- diagnostics -----
-
-template <int dim>
-void
-print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
-{
-  std::vector<types::global_dof_index> dof_indices;
-
-  for (const auto &cell : dofh.active_cell_iterators())
-    for (unsigned int f = 0; f < cell->n_faces(); ++f)
-      {
-        const auto &face = cell->face(f);
-
-        if (face->has_children())
-          {
-            for (unsigned int sf = 0; sf < face->n_children(); ++sf)
-              {
-                const auto &subface = face->child(sf);
-                Assert(subface->n_active_fe_indices() == 1, ExcInternalError());
-                const unsigned int subface_fe_index =
-                  subface->nth_active_fe_index(0);
-                const auto &subface_fe = dofh.get_fe(subface_fe_index);
-
-                dof_indices.resize(subface_fe.n_dofs_per_face(f));
-                subface->get_dof_indices(dof_indices, subface_fe_index);
-
-                deallog << "cell:" << cell->active_cell_index() << " face:" << f
-                        << " subface:" << sf << " dofs:";
-                for (const auto &i : dof_indices)
-                  deallog << i << " ";
-                deallog << std::endl;
-              }
-          }
-        else
-          {
-            Assert(face->n_active_fe_indices() == 1, ExcInternalError());
-            const unsigned int face_fe_index = face->nth_active_fe_index(0);
-            const auto &       face_fe       = dofh.get_fe(face_fe_index);
-
-            dof_indices.resize(face_fe.n_dofs_per_face(f));
-            face->get_dof_indices(dof_indices, face_fe_index);
-
-            deallog << "cell:" << cell->active_cell_index() << " face:" << f
-                    << " dofs:";
-            for (const auto &i : dof_indices)
-              deallog << i << " ";
-            deallog << std::endl;
-          }
-      }
-}
-
-
-template <int dim>
-void
-print_dof_points(const DoFHandler<dim> &dofh)
-{
-  std::vector<Point<dim>> points(dofh.n_dofs());
-  DoFTools::map_dofs_to_support_points(MappingFE<dim>(dofh.get_fe()),
-                                       dofh,
-                                       points);
-
-  for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
-    deallog << "dof:" << i << " point:" << points[i] << std::endl;
-}
-
-
-
-// ----- test -----
-
-template <int dim>
-void
-test()
-{
-  // setup grid
-  Triangulation<dim> tria;
-  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
-
-  tria.begin_active()->set_refine_flag();
-  tria.execute_coarsening_and_refinement();
-
-#if false
-  GridOut grid_out;
-  grid_out.write_vtk(tria, deallog.get_file_stream());
-#endif
-
-  DoFHandler<dim> dofh(tria);
-  dofh.distribute_dofs(Simplex::FE_P<dim>(1));
-  deallog << "ndofs: " << dofh.n_dofs() << std::endl;
-
-#if false
-  print_dof_points(dofh);
-  print_dof_indices_on_faces(dofh);
-#endif
-
-  // hanging node constraints
-  AffineConstraints<double> constraints;
-  DoFTools::make_hanging_node_constraints(dofh, constraints);
-  constraints.print(deallog.get_file_stream());
-
-  deallog << "OK" << std::endl;
-}
+#include "hanging_nodes.h"
 
 
 int
@@ -153,6 +51,22 @@ main()
   initlog();
 
   deallog.push("2d");
-  test<2>();
+  {
+    const unsigned int dim = 2;
+
+    const auto subdivided_hyper_cube_with_simplices =
+      [](Triangulation<dim> &tria) {
+        GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+      };
+
+    test<dim>({1, 0},
+              {0, 0},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(1)),
+              subdivided_hyper_cube_with_simplices);
+    test<dim>({1, 0},
+              {0, 0},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(2)),
+              subdivided_hyper_cube_with_simplices);
+  }
   deallog.pop();
 }
index 8034a9e25676a3dfd446579df6af0d1bfc979325..eb6b0d13acc87eea3161bfafcf0feca985014957 100644 (file)
@@ -1,5 +1,14 @@
 
 DEAL:2d::ndofs: 7
-    6 2:  0.500000
     6 1:  0.500000
+    6 2:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 19
+    12 4:  1.00000
+    14 1:  -0.125000
+    14 2:  0.375000
+    14 4:  0.750000
+    17 1:  0.375000
+    17 2:  -0.125000
+    17 4:  0.750000
 DEAL:2d::OK
index 024848424ed16f1f971be7270f9ee82ac8a80d58..838092d0f8127580753e455825bdc310534d9294 100644 (file)
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.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/grid/tria.h>
 
 #include <deal.II/hp/fe_collection.h>
-#include <deal.II/hp/mapping_collection.h>
 
 #include <deal.II/lac/affine_constraints.h>
 
 
 #include "../tests.h"
 
-
-// ----- diagnostics -----
-
-template <int dim>
-void
-print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
-{
-  std::vector<types::global_dof_index> dof_indices;
-
-  for (const auto &cell : dofh.active_cell_iterators())
-    for (unsigned int f = 0; f < cell->n_faces(); ++f)
-      {
-        const auto &face = cell->face(f);
-
-        Assert(!face->has_children(), ExcInternalError());
-
-        const unsigned int fe_index = cell->active_fe_index();
-        const auto &       fe       = cell->get_fe();
-
-        dof_indices.resize(fe.n_dofs_per_face(f));
-        face->get_dof_indices(dof_indices, fe_index);
-
-        deallog << "cell:" << cell->active_cell_index() << " face:" << f
-                << " dofs:";
-        for (const auto &i : dof_indices)
-          deallog << i << " ";
-        deallog << std::endl;
-      }
-}
-
-
-template <int dim>
-void
-print_dof_points(const DoFHandler<dim> &dofh)
-{
-  hp::MappingCollection<dim> mapping;
-  for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
-    mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
-
-  std::vector<Point<dim>> points(dofh.n_dofs());
-  DoFTools::map_dofs_to_support_points(mapping, dofh, points);
-
-  for (unsigned int i = 0; i < dofh.n_dofs(); ++i)
-    deallog << "dof:" << i << " point:" << points[i] << std::endl;
-}
-
-
-// ----- test -----
-
-template <int dim>
-void
-test(const hp::FECollection<dim> &fes)
-{
-  // setup grid
-  Triangulation<dim> tria;
-  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
-
-#if false
-  GridOut grid_out;
-  grid_out.write_vtk(tria, deallog.get_file_stream());
-#endif
-
-  DoFHandler<dim> dofh(tria);
-  dofh.begin_active()->set_active_fe_index(1);
-
-  dofh.distribute_dofs(fes);
-  deallog << "ndofs: " << dofh.n_dofs() << std::endl;
-
-#if false
-  print_dof_points(dofh);
-  print_dof_indices_on_faces(dofh);
-#endif
-
-  // hanging node constraints
-  AffineConstraints<double> constraints;
-  DoFTools::make_hanging_node_constraints(dofh, constraints);
-  constraints.print(deallog.get_file_stream());
-
-  deallog << "OK" << std::endl;
-}
+#include "hanging_nodes.h"
 
 
 int
@@ -133,7 +53,24 @@ main()
   initlog();
 
   deallog.push("2d");
-  test<2>(hp::FECollection<2>(Simplex::FE_P<2>(1), Simplex::FE_P<2>(2)));
-  test<2>(hp::FECollection<2>(Simplex::FE_P<2>(2), Simplex::FE_P<2>(1)));
+  {
+    const unsigned int dim = 2;
+
+    const auto subdivided_hyper_cube_with_simplices =
+      [](Triangulation<dim> &tria) {
+        GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+      };
+
+    test<dim>({0, 0},
+              {0, 1},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(2),
+                                    Simplex::FE_P<dim>(1)),
+              subdivided_hyper_cube_with_simplices);
+    test<dim>({0, 0},
+              {0, 1},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(1),
+                                    Simplex::FE_P<dim>(2)),
+              subdivided_hyper_cube_with_simplices);
+  }
   deallog.pop();
 }
index 0bc8c9bb95bfacfadc9f932e99c53528ef8346d7..bcddbbf9013654581f07732d14a2a03029f7cb69 100644 (file)
@@ -1,7 +1,7 @@
 
 DEAL:2d::ndofs: 7
-    2 6:  0.500000
     2 5:  0.500000
+    2 6:  0.500000
 DEAL:2d::OK
 DEAL:2d::ndofs: 7
     5 1:  0.500000
diff --git a/tests/simplex/hanging_nodes_03.cc b/tests/simplex/hanging_nodes_03.cc
new file mode 100644 (file)
index 0000000..c0e802b
--- /dev/null
@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// verify hanging node constraints on locally hp-refined simplex mesh
+//
+// dofs will be enumerated as follows
+//  scenario 1:
+//   9---1---0
+//   |\      |
+//   |  \    |
+//   6--2,8  3
+//   |\  |\  |
+//   |  \|  \|
+//   4---5---7
+
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+#include "hanging_nodes.h"
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    const unsigned int dim = 2;
+
+    const auto subdivided_hyper_cube_with_simplices =
+      [](Triangulation<dim> &tria) {
+        GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+      };
+
+    test<dim>({1, 0},
+              {0, 1},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(1),
+                                    Simplex::FE_P<dim>(2)),
+              subdivided_hyper_cube_with_simplices);
+    test<dim>({1, 0},
+              {0, 1},
+              hp::FECollection<dim>(Simplex::FE_P<dim>(2),
+                                    Simplex::FE_P<dim>(1)),
+              subdivided_hyper_cube_with_simplices);
+  }
+  deallog.pop();
+}
diff --git a/tests/simplex/hanging_nodes_03.with_simplex_support=on.output b/tests/simplex/hanging_nodes_03.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..589460a
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL:2d::ndofs: 10
+    2 7:  0.500000
+    2 9:  0.500000
+    8 7:  0.500000
+    8 9:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 16
+    9 1:  0.500000
+    9 2:  0.500000
+    11 1:  0.250000
+    11 2:  0.750000
+    14 1:  0.750000
+    14 2:  0.250000
+DEAL:2d::OK
diff --git a/tests/simplex/hanging_nodes_hybrid_01.cc b/tests/simplex/hanging_nodes_hybrid_01.cc
new file mode 100644 (file)
index 0000000..6a98030
--- /dev/null
@@ -0,0 +1,83 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// verify hanging node constraints on locally h-refined hybrid meshes
+/*
+ * dofs will be enumerated as follows with d = 1
+ *  scenario 1:        scenario 2:        scenario 3:
+ *   2-------3          2-------3          7---8---9
+ *   |       |\         |       |\         |   |   |\
+ *   |       | \        |       | 6        |   |   | \
+ *   |       |  \       |       |/|\       |   |   |  \
+ *   |       |   4      |       4 | 7      3---4---6   0
+ *   |       |  /       |       |\|/       |   |   |  /
+ *   |       | /        |       | 5        |   |   | /
+ *   |       |/         |       |/         |   |   |/
+ *   0-------1          0-------1          1---2---5
+ */
+
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+#include "hanging_nodes.h"
+#include "simplex_grids.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    const unsigned int dim = 2;
+
+    const auto cube_and_pyramid = [](Triangulation<dim> &tria) {
+      GridGenerator::cube_and_pyramid(tria, 1);
+    };
+
+    for (unsigned int d = 1; d <= 2; ++d)
+      {
+        deallog << "degree: " << d << std::endl;
+        test<dim>({0, 1},
+                  {0, 1},
+                  hp::FECollection<dim>(FE_Q<dim>(d), Simplex::FE_P<dim>(d)),
+                  cube_and_pyramid);
+        test<dim>({1, 0},
+                  {0, 1},
+                  hp::FECollection<dim>(FE_Q<dim>(d), Simplex::FE_P<dim>(d)),
+                  cube_and_pyramid);
+        test<dim>({0, 1},
+                  {1, 0},
+                  hp::FECollection<dim>(Simplex::FE_P<dim>(d), FE_Q<dim>(d)),
+                  cube_and_pyramid);
+        test<dim>({1, 0},
+                  {1, 0},
+                  hp::FECollection<dim>(Simplex::FE_P<dim>(d), FE_Q<dim>(d)),
+                  cube_and_pyramid);
+      }
+  }
+  deallog.pop();
+}
diff --git a/tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..3ff5c89
--- /dev/null
@@ -0,0 +1,55 @@
+
+DEAL:2d::degree: 1
+DEAL:2d::ndofs: 8
+    4 1:  0.500000
+    4 3:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 10
+    6 5:  0.500000
+    6 9:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 8
+    3 2:  0.500000
+    3 5:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 10
+    7 0:  0.500000
+    7 1:  0.500000
+DEAL:2d::OK
+DEAL:2d::degree: 2
+DEAL:2d::ndofs: 22
+    9 5:  1.00000
+    11 1:  0.375000
+    11 3:  -0.125000
+    11 5:  0.750000
+    15 1:  -0.125000
+    15 3:  0.375000
+    15 5:  0.750000
+DEAL:2d::OK
+DEAL:2d::ndofs: 29
+    14 1:  1.00000
+    15 1:  0.750000
+    15 13:  0.375000
+    15 25:  -0.125000
+    26 1:  0.750000
+    26 13:  -0.125000
+    26 25:  0.375000
+DEAL:2d::OK
+DEAL:2d::ndofs: 22
+    8 3:  1.00000
+    10 3:  0.750000
+    10 7:  0.375000
+    10 13:  -0.125000
+    15 3:  0.750000
+    15 7:  -0.125000
+    15 13:  0.375000
+DEAL:2d::OK
+DEAL:2d::ndofs: 29
+    15 3:  1.00000
+    16 0:  0.375000
+    16 1:  -0.125000
+    16 3:  0.750000
+    26 0:  -0.125000
+    26 1:  0.375000
+    26 3:  0.750000
+DEAL:2d::OK
diff --git a/tests/simplex/hanging_nodes_hybrid_02.cc b/tests/simplex/hanging_nodes_hybrid_02.cc
new file mode 100644 (file)
index 0000000..cf7dc74
--- /dev/null
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// verify hanging node constraints on locally p-refined hybrid meshes
+/*
+ * dofs will be enumerated as follows for degrees = {1, 2}
+ *  scenario 1:        scenario 2:
+ *   2-------3          1---5---8
+ *   |       |\         |       |\
+ *   |       | 6        |       | \
+ *   |       |  \       |       |  \
+ *   |       5   4      2   6   3   9
+ *   |       |  /       |       |  /
+ *   |       | 7        |       | /
+ *   |       |/         |       |/
+ *   0-------1          0---4---7
+ */
+
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+#include "hanging_nodes.h"
+#include "simplex_grids.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    const unsigned int dim = 2;
+
+    const auto cube_and_pyramid = [](Triangulation<dim> &tria) {
+      GridGenerator::cube_and_pyramid(tria, 1);
+    };
+
+    for (unsigned int q = 1; q <= 4; ++q)
+      for (unsigned int p = 1; p <= 2; ++p)
+        {
+          if (q == p)
+            continue;
+
+          deallog << "q_degree: " << q << ", p_degree: " << p << std::endl;
+          test<dim>({0, 0},
+                    {0, 1},
+                    hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+                    cube_and_pyramid);
+          test<dim>({0, 0},
+                    {1, 0},
+                    hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+                    cube_and_pyramid);
+        }
+  }
+  deallog.pop();
+}
diff --git a/tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..1d51692
--- /dev/null
@@ -0,0 +1,83 @@
+
+DEAL:2d::q_degree: 1, p_degree: 2
+DEAL:2d::ndofs: 8
+    5 1:  0.500000
+    5 3:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 8
+    5 1:  0.500000
+    5 3:  0.500000
+DEAL:2d::OK
+DEAL:2d::q_degree: 2, p_degree: 1
+DEAL:2d::ndofs: 10
+    3 7:  0.500000
+    3 8:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 10
+    3 7:  0.500000
+    3 8:  0.500000
+DEAL:2d::OK
+DEAL:2d::q_degree: 3, p_degree: 1
+DEAL:2d::ndofs: 17
+    4 14:  0.723607
+    4 15:  0.276393
+    5 14:  0.276393
+    5 15:  0.723607
+DEAL:2d::OK
+DEAL:2d::ndofs: 17
+    4 14:  0.723607
+    4 15:  0.276393
+    5 14:  0.276393
+    5 15:  0.723607
+DEAL:2d::OK
+DEAL:2d::q_degree: 3, p_degree: 2
+DEAL:2d::ndofs: 20
+    4 14:  0.323607
+    4 15:  -0.123607
+    4 17:  0.800000
+    5 14:  -0.123607
+    5 15:  0.323607
+    5 17:  0.800000
+DEAL:2d::OK
+DEAL:2d::ndofs: 20
+    4 14:  0.323607
+    4 15:  -0.123607
+    4 17:  0.800000
+    5 14:  -0.123607
+    5 15:  0.323607
+    5 17:  0.800000
+DEAL:2d::OK
+DEAL:2d::q_degree: 4, p_degree: 1
+DEAL:2d::ndofs: 26
+    5 23:  0.827327
+    5 24:  0.172673
+    6 23:  0.500000
+    6 24:  0.500000
+    7 23:  0.172673
+    7 24:  0.827327
+DEAL:2d::OK
+DEAL:2d::ndofs: 26
+    5 23:  0.827327
+    5 24:  0.172673
+    6 23:  0.500000
+    6 24:  0.500000
+    7 23:  0.172673
+    7 24:  0.827327
+DEAL:2d::OK
+DEAL:2d::q_degree: 4, p_degree: 2
+DEAL:2d::ndofs: 28
+    5 22:  0.541613
+    5 23:  -0.113041
+    5 25:  0.571429
+    6 22:  -0.113041
+    6 23:  0.541613
+    6 25:  0.571429
+DEAL:2d::OK
+DEAL:2d::ndofs: 28
+    5 22:  0.541613
+    5 23:  -0.113041
+    5 25:  0.571429
+    6 22:  -0.113041
+    6 23:  0.541613
+    6 25:  0.571429
+DEAL:2d::OK
diff --git a/tests/simplex/hanging_nodes_hybrid_03.cc b/tests/simplex/hanging_nodes_hybrid_03.cc
new file mode 100644 (file)
index 0000000..63440a7
--- /dev/null
@@ -0,0 +1,74 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// verify hanging node constraints on locally hp-refined hybrid meshes
+
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+#include "hanging_nodes.h"
+#include "simplex_grids.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    const unsigned int dim = 2;
+
+    const auto cube_and_pyramid = [](Triangulation<dim> &tria) {
+      GridGenerator::cube_and_pyramid(tria, 1);
+    };
+
+    for (unsigned int q = 1; q <= 4; ++q)
+      for (unsigned int p = 1; p <= 2; ++p)
+        {
+          if (q == p)
+            continue;
+
+          deallog << "q_degree: " << q << ", p_degree: " << p << std::endl;
+          test<dim>({0, 1},
+                    {0, 1},
+                    hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+                    cube_and_pyramid);
+          test<dim>({1, 0},
+                    {0, 1},
+                    hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+                    cube_and_pyramid);
+          test<dim>({0, 1},
+                    {1, 0},
+                    hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+                    cube_and_pyramid);
+          test<dim>({1, 0},
+                    {1, 0},
+                    hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+                    cube_and_pyramid);
+        }
+  }
+  deallog.pop();
+}
diff --git a/tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..cfc2558
--- /dev/null
@@ -0,0 +1,289 @@
+
+DEAL:2d::q_degree: 1, p_degree: 2
+DEAL:2d::ndofs: 17
+    4 1:  0.500000
+    4 3:  0.500000
+    6 1:  0.750000
+    6 3:  0.250000
+    10 1:  0.250000
+    10 3:  0.750000
+DEAL:2d::OK
+DEAL:2d::ndofs: 13
+    1 8:  0.500000
+    1 12:  0.500000
+    9 8:  0.500000
+    9 12:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 17
+    4 1:  0.500000
+    4 3:  0.500000
+    6 1:  0.750000
+    6 3:  0.250000
+    10 1:  0.250000
+    10 3:  0.750000
+DEAL:2d::OK
+DEAL:2d::ndofs: 13
+    1 8:  0.500000
+    1 12:  0.500000
+    9 8:  0.500000
+    9 12:  0.500000
+DEAL:2d::OK
+DEAL:2d::q_degree: 2, p_degree: 1
+DEAL:2d::ndofs: 13
+    3 7:  0.500000
+    3 10:  0.500000
+    8 7:  0.500000
+    8 10:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 26
+    12 0:  0.500000
+    12 1:  0.500000
+    13 0:  0.750000
+    13 1:  0.250000
+    23 0:  0.250000
+    23 1:  0.750000
+DEAL:2d::OK
+DEAL:2d::ndofs: 13
+    3 7:  0.500000
+    3 10:  0.500000
+    8 7:  0.500000
+    8 10:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 26
+    12 0:  0.500000
+    12 1:  0.500000
+    13 0:  0.750000
+    13 1:  0.250000
+    23 0:  0.250000
+    23 1:  0.750000
+DEAL:2d::OK
+DEAL:2d::q_degree: 3, p_degree: 1
+DEAL:2d::ndofs: 20
+    4 14:  0.723607
+    4 17:  0.276393
+    5 14:  0.276393
+    5 17:  0.723607
+    15 14:  0.500000
+    15 17:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 50
+    19 0:  0.500000
+    19 1:  0.500000
+    20 0:  0.861803
+    20 1:  0.138197
+    21 0:  0.638197
+    21 1:  0.361803
+    42 0:  0.361803
+    42 1:  0.638197
+    43 0:  0.138197
+    43 1:  0.861803
+DEAL:2d::OK
+DEAL:2d::ndofs: 20
+    4 14:  0.723607
+    4 17:  0.276393
+    5 14:  0.276393
+    5 17:  0.723607
+    15 14:  0.500000
+    15 17:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 50
+    19 0:  0.500000
+    19 1:  0.500000
+    20 0:  0.861803
+    20 1:  0.138197
+    21 0:  0.638197
+    21 1:  0.361803
+    42 0:  0.361803
+    42 1:  0.638197
+    43 0:  0.138197
+    43 1:  0.861803
+DEAL:2d::OK
+DEAL:2d::q_degree: 3, p_degree: 2
+DEAL:2d::ndofs: 29
+    5 4:  1.00000
+    5 14:  -0.447214
+    5 20:  0.447214
+    15 4:  1.25000
+    15 14:  -0.404508
+    15 20:  0.154508
+    17 4:  0.937500
+    17 14:  0.0716186
+    17 20:  -0.00911863
+    22 4:  0.937500
+    22 14:  -0.428381
+    22 20:  0.490881
+DEAL:2d::OK
+DEAL:2d::ndofs: 53
+    22 3:  1.00000
+    23 0:  0.623607
+    23 1:  -0.100000
+    23 3:  0.476393
+    24 0:  0.176393
+    24 1:  -0.100000
+    24 3:  0.923607
+    45 0:  -0.100000
+    45 1:  0.176393
+    45 3:  0.923607
+    46 0:  -0.100000
+    46 1:  0.623607
+    46 3:  0.476393
+DEAL:2d::OK
+DEAL:2d::ndofs: 29
+    5 4:  1.00000
+    5 14:  -0.447214
+    5 20:  0.447214
+    15 4:  1.25000
+    15 14:  -0.404508
+    15 20:  0.154508
+    17 4:  0.937500
+    17 14:  0.0716186
+    17 20:  -0.00911863
+    22 4:  0.937500
+    22 14:  -0.428381
+    22 20:  0.490881
+DEAL:2d::OK
+DEAL:2d::ndofs: 53
+    22 3:  1.00000
+    23 0:  0.623607
+    23 1:  -0.100000
+    23 3:  0.476393
+    24 0:  0.176393
+    24 1:  -0.100000
+    24 3:  0.923607
+    45 0:  -0.100000
+    45 1:  0.176393
+    45 3:  0.923607
+    46 0:  -0.100000
+    46 1:  0.623607
+    46 3:  0.476393
+DEAL:2d::OK
+DEAL:2d::q_degree: 4, p_degree: 1
+DEAL:2d::ndofs: 29
+    5 23:  0.827327
+    5 26:  0.172673
+    6 23:  0.500000
+    6 26:  0.500000
+    7 23:  0.172673
+    7 26:  0.827327
+    24 23:  0.500000
+    24 26:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 82
+    28 0:  0.500000
+    28 1:  0.500000
+    29 0:  0.913663
+    29 1:  0.0863366
+    30 0:  0.750000
+    30 1:  0.250000
+    31 0:  0.586337
+    31 1:  0.413663
+    67 0:  0.413663
+    67 1:  0.586337
+    68 0:  0.250000
+    68 1:  0.750000
+    69 0:  0.0863366
+    69 1:  0.913663
+DEAL:2d::OK
+DEAL:2d::ndofs: 29
+    5 23:  0.827327
+    5 26:  0.172673
+    6 23:  0.500000
+    6 26:  0.500000
+    7 23:  0.172673
+    7 26:  0.827327
+    24 23:  0.500000
+    24 26:  0.500000
+DEAL:2d::OK
+DEAL:2d::ndofs: 82
+    28 0:  0.500000
+    28 1:  0.500000
+    29 0:  0.913663
+    29 1:  0.0863366
+    30 0:  0.750000
+    30 1:  0.250000
+    31 0:  0.586337
+    31 1:  0.413663
+    67 0:  0.413663
+    67 1:  0.586337
+    68 0:  0.250000
+    68 1:  0.750000
+    69 0:  0.0863366
+    69 1:  0.913663
+DEAL:2d::OK
+DEAL:2d::q_degree: 4, p_degree: 2
+DEAL:2d::ndofs: 38
+    6 5:  1.75000
+    6 23:  -0.947822
+    6 29:  0.197822
+    7 5:  1.00000
+    7 23:  -0.654654
+    7 29:  0.654654
+    24 5:  1.75000
+    24 23:  -0.947822
+    24 29:  0.197822
+    26 5:  1.31250
+    26 23:  -0.335866
+    26 29:  0.0233665
+    31 5:  1.31250
+    31 23:  -0.835866
+    31 29:  0.523366
+DEAL:2d::OK
+DEAL:2d::ndofs: 85
+    31 3:  1.00000
+    32 0:  0.755898
+    32 1:  -0.0714286
+    32 3:  0.315530
+    33 0:  0.375000
+    33 1:  -0.125000
+    33 3:  0.750000
+    34 0:  0.101245
+    34 1:  -0.0714286
+    34 3:  0.970184
+    70 0:  -0.0714286
+    70 1:  0.101245
+    70 3:  0.970184
+    71 0:  -0.125000
+    71 1:  0.375000
+    71 3:  0.750000
+    72 0:  -0.0714286
+    72 1:  0.755898
+    72 3:  0.315530
+DEAL:2d::OK
+DEAL:2d::ndofs: 38
+    6 5:  1.75000
+    6 23:  -0.947822
+    6 29:  0.197822
+    7 5:  1.00000
+    7 23:  -0.654654
+    7 29:  0.654654
+    24 5:  1.75000
+    24 23:  -0.947822
+    24 29:  0.197822
+    26 5:  1.31250
+    26 23:  -0.335866
+    26 29:  0.0233665
+    31 5:  1.31250
+    31 23:  -0.835866
+    31 29:  0.523366
+DEAL:2d::OK
+DEAL:2d::ndofs: 85
+    31 3:  1.00000
+    32 0:  0.755898
+    32 1:  -0.0714286
+    32 3:  0.315530
+    33 0:  0.375000
+    33 1:  -0.125000
+    33 3:  0.750000
+    34 0:  0.101245
+    34 1:  -0.0714286
+    34 3:  0.970184
+    70 0:  -0.0714286
+    70 1:  0.101245
+    70 3:  0.970184
+    71 0:  -0.125000
+    71 1:  0.375000
+    71 3:  0.750000
+    72 0:  -0.0714286
+    72 1:  0.755898
+    72 3:  0.315530
+DEAL:2d::OK
index 8d0e42dc1a14071da9b57c869a155d4303735cdd..48e08c6697213eff43940f9dbed55843a3c9d821 100644 (file)
@@ -376,5 +376,60 @@ namespace dealii
           AssertThrow(false, ExcNotImplemented())
         }
     }
+
+
+
+    /**
+     * Adds a simplex cell to face @p face_no of a hyper_cube cell.
+     */
+    template <int dim, int spacedim>
+    void
+    cube_and_pyramid(Triangulation<dim, spacedim> &tria,
+                     const unsigned int            face_no = 1)
+    {
+      Assert(face_no % 2 == 1,
+             ExcMessage("Only works for odd face numbers. "
+                        "GridReordering::reorder_cells() is not prepared for "
+                        "simplices yet (uses GeometryInfo)."));
+
+      // create cube
+      Triangulation<dim, spacedim> tria_cube;
+      GridGenerator::hyper_cube(tria_cube);
+      const auto cube = tria_cube.begin_active();
+      AssertIndexRange(face_no, cube->n_faces());
+
+      // extract vertices from specified face, store their midpoint
+      std::vector<Point<spacedim>> vertices;
+      Point<spacedim>              midpoint;
+      const auto                   shared_face = cube->face(face_no);
+      for (unsigned int v = 0; v < shared_face->n_vertices(); ++v)
+        {
+          const auto &vertex = shared_face->vertex(v);
+          vertices.push_back(vertex);
+          midpoint += vertex;
+        }
+      midpoint /= vertices.size();
+
+      // add simplex cell in coordinate direction
+      const unsigned int coordinate = face_no / 2;
+#ifdef DEBUG
+      // all vertices should be in a plane
+      for (const auto &vertex : vertices)
+        Assert(midpoint(coordinate) == vertex(coordinate), ExcInternalError());
+#endif
+
+      // add another vertex as tip of triangle/pyramid
+      Point<spacedim> tip = midpoint;
+      tip(coordinate) += (face_no % 2 == 1) ? 1. : -1.;
+      vertices.push_back(tip);
+
+      CellData<dim> simplex(vertices.size());
+      std::iota(simplex.vertices.begin(), simplex.vertices.end(), 0);
+
+      Triangulation<dim, spacedim> tria_simplex;
+      tria_simplex.create_triangulation(vertices, {simplex}, SubCellData());
+
+      GridGenerator::merge_triangulations(tria_cube, tria_simplex, tria);
+    }
   } // namespace GridGenerator
 } // namespace dealii

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.