MatrixFree: support FE_PyramidP

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 403e25aceb3c391d8b3def19180584b89abc19de..9a0bcf9b71eddad1328a137eea4fa22100c48a84 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -93,7 +93,8 @@ namespace internal
 
       if (dynamic_cast<const Simplex::FE_P<dim, dim> *>(&fe) != nullptr ||
           dynamic_cast<const Simplex::FE_DGP<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const Simplex::FE_WedgeP<dim, dim> *>(&fe) != nullptr)
+          dynamic_cast<const Simplex::FE_WedgeP<dim, dim> *>(&fe) != nullptr ||
+          dynamic_cast<const Simplex::FE_PyramidP<dim, dim> *>(&fe) != nullptr)
         {
           scalar_lexicographic.resize(fe.n_dofs_per_cell());
           for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i)
@@ -206,7 +207,8 @@ namespace internal
               // indicative of their support
               if (dynamic_cast<const Simplex::FE_P<dim> *>(fe_poly_ptr) ||
                   dynamic_cast<const Simplex::FE_DGP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const Simplex::FE_WedgeP<dim> *>(fe_poly_ptr))
+                  dynamic_cast<const Simplex::FE_WedgeP<dim> *>(fe_poly_ptr) ||
+                  dynamic_cast<const Simplex::FE_PyramidP<dim> *>(fe_poly_ptr))
                 return true;
 #endif
 
@@ -247,6 +249,8 @@ namespace internal
           dynamic_cast<const Simplex::FE_DGP<dim> *>(
             &fe_in.base_element(base_element_number)) ||
           dynamic_cast<const Simplex::FE_WedgeP<dim> *>(
+            &fe_in.base_element(base_element_number)) ||
+          dynamic_cast<const Simplex::FE_PyramidP<dim> *>(
             &fe_in.base_element(base_element_number)))
         {
           // specialization for arbitrary finite elements and quadrature rules

diff --git a/tests/simplex/matrix_free_01.cc b/tests/simplex/matrix_free_01.cc
index 8f0df19794a5fc3fe4ee62ceb6865ffdb6c0b7f8..daeeb53115599067e14ec2322f22df85a0e28bb4 100644 (file)
--- a/tests/simplex/matrix_free_01.cc
+++ b/tests/simplex/matrix_free_01.cc
@@ -148,6 +148,14 @@ test(const unsigned int v, const unsigned int degree, const bool do_helmholtz)
       quad       = std::make_shared<Simplex::QGaussWedge<dim>>(degree + 1);
       fe_mapping = std::make_shared<Simplex::FE_WedgeP<dim>>(1);
     }
+  else if (v == 2)
+    {
+      GridGenerator::subdivided_hyper_cube_with_pyramids(tria,
+                                                         dim == 2 ? 16 : 8);
+      fe         = std::make_shared<Simplex::FE_PyramidP<dim>>(degree);
+      quad       = std::make_shared<Simplex::QGaussPyramid<dim>>(degree + 1);
+      fe_mapping = std::make_shared<Simplex::FE_PyramidP<dim>>(1);
+    }
   else
     Assert(false, ExcNotImplemented());
 
@@ -304,16 +312,18 @@ main(int argc, char **argv)
 
   Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
 
-  for (unsigned int i = 0; i < 2; ++i)
+  for (unsigned int i = 0; i <= 2; ++i)
     {
       if (i == 0)
         deallog.push("SIMPLEX");
       else if (i == 1)
         deallog.push("WEDGE  ");
+      else if (i == 2)
+        deallog.push("PYRAMID");
       else
         Assert(false, ExcNotImplemented());
 
-      if (i == 0)
+      if (i == 0) // 2D makes only sense for simplex
         {
           test<2>(i, /*degree=*/1, /*do_helmholtz*/ false);
           test<2>(i, /*degree=*/1, /*do_helmholtz*/ true);
@@ -323,8 +333,13 @@ main(int argc, char **argv)
 
       test<3>(i, /*degree=*/1, /*do_helmholtz*/ false);
       test<3>(i, /*degree=*/1, /*do_helmholtz*/ true);
-      test<3>(i, /*degree=*/2, /*do_helmholtz*/ false);
-      test<3>(i, /*degree=*/2, /*do_helmholtz*/ true);
+
+      if (i !=
+          2) // for pyramids no quadratic elements have been implemented yet
+        {
+          test<3>(i, /*degree=*/2, /*do_helmholtz*/ false);
+          test<3>(i, /*degree=*/2, /*do_helmholtz*/ true);
+        }
 
       deallog.pop();
     }

diff --git a/tests/simplex/matrix_free_01.with_simplex_support=on.output b/tests/simplex/matrix_free_01.with_simplex_support=on.output
index 5bcc1f1895a71c9f80fad8a8df0db0f0f68cf207..83ee05c055cde84d0cea4d86dedecccc50d07d2f 100644 (file)
--- a/tests/simplex/matrix_free_01.with_simplex_support=on.output
+++ b/tests/simplex/matrix_free_01.with_simplex_support=on.output
@@ -11,3 +11,5 @@ DEAL:WEDGE  ::dim=3 degree=1 Type=Possion   : Convergence step 7 value 0.0002242
 DEAL:WEDGE  ::dim=3 degree=1 Type=Helmholtz : Convergence step 7 value 0.000211401.
 DEAL:WEDGE  ::dim=3 degree=2 Type=Possion   : Convergence step 20 value 0.000133408.
 DEAL:WEDGE  ::dim=3 degree=2 Type=Helmholtz : Convergence step 20 value 0.000125844.
+DEAL:PYRAMID::dim=3 degree=1 Type=Possion   : Convergence step 10 value 0.000197898.
+DEAL:PYRAMID::dim=3 degree=1 Type=Helmholtz : Convergence step 10 value 0.000186534.

diff --git a/tests/simplex/simplex_grids.h b/tests/simplex/simplex_grids.h
index 7c916003108193aedb7f07c36300ee929458191c..8d0e42dc1a14071da9b57c869a155d4303735cdd 100644 (file)
--- a/tests/simplex/simplex_grids.h
+++ b/tests/simplex/simplex_grids.h
@@ -250,6 +250,31 @@ namespace dealii
 
 
 
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_cube_with_pyramids(Triangulation<dim, spacedim> &tria,
+                                        const unsigned int repetitions,
+                                        const double       p1       = 0.0,
+                                        const double       p2       = 1.0,
+                                        const bool         colorize = false)
+    {
+      if (dim == 3)
+        {
+          subdivided_hyper_rectangle_with_pyramids(
+            tria,
+            {{repetitions, repetitions, repetitions}},
+            {p1, p1, p1},
+            {p2, p2, p2},
+            colorize);
+        }
+      else
+        {
+          AssertThrow(false, ExcNotImplemented())
+        }
+    }
+
+
+
     template <int dim, int spacedim>
     void
     subdivided_hyper_rectangle_with_simplices_mix(