Add Poisson/Helmholtz test: MatrixFree + wedge mesh

diff --git a/tests/simplex/matrix_free_01.cc b/tests/simplex/matrix_free_01.cc
index 94e0f4031d7bb67d2e99ddf81f787cbafee6d325..8f0df19794a5fc3fe4ee62ceb6865ffdb6c0b7f8 100644 (file)
--- a/tests/simplex/matrix_free_01.cc
+++ b/tests/simplex/matrix_free_01.cc
@@ -45,6 +45,8 @@
 
 #include "../tests.h"
 
+#include "./simplex_grids.h"
+
 using namespace dealii;
 
 
@@ -123,22 +125,44 @@ private:
 
 template <int dim>
 void
-test(const unsigned int degree, const bool do_helmholtz)
+test(const unsigned int v, const unsigned int degree, const bool do_helmholtz)
 {
   Triangulation<dim> tria;
 
-  GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 16 : 8);
-
-  Simplex::FE_P<dim>   fe(degree);
-  Simplex::QGauss<dim> quad(degree + 1);
-  MappingFE<dim>       mapping(Simplex::FE_P<dim>(1));
+  std::shared_ptr<FiniteElement<dim>> fe;
+  std::shared_ptr<Quadrature<dim>>    quad;
+  std::shared_ptr<FiniteElement<dim>> fe_mapping;
+
+  if (v == 0)
+    {
+      GridGenerator::subdivided_hyper_cube_with_simplices(tria,
+                                                          dim == 2 ? 16 : 8);
+      fe         = std::make_shared<Simplex::FE_P<dim>>(degree);
+      quad       = std::make_shared<Simplex::QGauss<dim>>(degree + 1);
+      fe_mapping = std::make_shared<Simplex::FE_P<dim>>(1);
+    }
+  else if (v == 1)
+    {
+      GridGenerator::subdivided_hyper_cube_with_wedges(tria, dim == 2 ? 16 : 8);
+      fe         = std::make_shared<Simplex::FE_WedgeP<dim>>(degree);
+      quad       = std::make_shared<Simplex::QGaussWedge<dim>>(degree + 1);
+      fe_mapping = std::make_shared<Simplex::FE_WedgeP<dim>>(1);
+    }
+  else
+    Assert(false, ExcNotImplemented());
+
+  MappingFE<dim> mapping(*fe_mapping);
 
   DoFHandler<dim> dof_handler(tria);
-  dof_handler.distribute_dofs(fe);
+  dof_handler.distribute_dofs(*fe);
 
   AffineConstraints<double> constraints;
+#if false
   VectorTools::interpolate_boundary_values(
     mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+#else
+  DoFTools::make_zero_boundary_constraints(dof_handler, 0, constraints);
+#endif
   constraints.close();
 
   const auto solve_and_postprocess =
@@ -173,7 +197,7 @@ test(const unsigned int degree, const bool do_helmholtz)
 
     MatrixFree<dim, double> matrix_free;
     matrix_free.reinit(
-      mapping, dof_handler, constraints, quad, additional_data);
+      mapping, dof_handler, constraints, *quad, additional_data);
 
     PoissonOperator<dim> poisson_operator(matrix_free, do_helmholtz);
 
@@ -203,7 +227,7 @@ test(const unsigned int degree, const bool do_helmholtz)
 
     const auto flags = update_values | update_gradients | update_JxW_values;
 
-    FEValues<dim> fe_values(mapping, fe, quad, flags);
+    FEValues<dim> fe_values(mapping, *fe, *quad, flags);
 
     FullMatrix<double>                   cell_matrix;
     Vector<double>                       cell_rhs;
@@ -276,17 +300,32 @@ main(int argc, char **argv)
 {
   initlog();
 
-  deallog.depth_file(1);
+  deallog.depth_file(2);
 
   Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
 
-  test<2>(/*degree=*/1, /*do_helmholtz*/ false);
-  test<2>(/*degree=*/1, /*do_helmholtz*/ true);
-  test<2>(/*degree=*/2, /*do_helmholtz*/ false);
-  test<2>(/*degree=*/2, /*do_helmholtz*/ true);
-
-  test<3>(/*degree=*/1, /*do_helmholtz*/ false);
-  test<3>(/*degree=*/1, /*do_helmholtz*/ true);
-  test<3>(/*degree=*/2, /*do_helmholtz*/ false);
-  test<3>(/*degree=*/2, /*do_helmholtz*/ true);
+  for (unsigned int i = 0; i < 2; ++i)
+    {
+      if (i == 0)
+        deallog.push("SIMPLEX");
+      else if (i == 1)
+        deallog.push("WEDGE  ");
+      else
+        Assert(false, ExcNotImplemented());
+
+      if (i == 0)
+        {
+          test<2>(i, /*degree=*/1, /*do_helmholtz*/ false);
+          test<2>(i, /*degree=*/1, /*do_helmholtz*/ true);
+          test<2>(i, /*degree=*/2, /*do_helmholtz*/ false);
+          test<2>(i, /*degree=*/2, /*do_helmholtz*/ true);
+        }
+
+      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);
+
+      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 ebff9c041c7b0d62c7838415273ea7df499770ab..5bcc1f1895a71c9f80fad8a8df0db0f0f68cf207 100644 (file)
--- a/tests/simplex/matrix_free_01.with_simplex_support=on.output
+++ b/tests/simplex/matrix_free_01.with_simplex_support=on.output
@@ -1,9 +1,13 @@
 
-DEAL::dim=2 degree=1 Type=Possion   : Convergence step 15 value 0.000430255.
-DEAL::dim=2 degree=1 Type=Helmholtz : Convergence step 14 value 0.000558792.
-DEAL::dim=2 degree=2 Type=Possion   : Convergence step 36 value 0.000319106.
-DEAL::dim=2 degree=2 Type=Helmholtz : Convergence step 36 value 0.000294440.
-DEAL::dim=3 degree=1 Type=Possion   : Convergence step 7 value 0.000281643.
-DEAL::dim=3 degree=1 Type=Helmholtz : Convergence step 7 value 0.000268188.
-DEAL::dim=3 degree=2 Type=Possion   : Convergence step 20 value 0.000146315.
-DEAL::dim=3 degree=2 Type=Helmholtz : Convergence step 20 value 0.000137731.
+DEAL:SIMPLEX::dim=2 degree=1 Type=Possion   : Convergence step 15 value 0.000430255.
+DEAL:SIMPLEX::dim=2 degree=1 Type=Helmholtz : Convergence step 14 value 0.000558792.
+DEAL:SIMPLEX::dim=2 degree=2 Type=Possion   : Convergence step 36 value 0.000319106.
+DEAL:SIMPLEX::dim=2 degree=2 Type=Helmholtz : Convergence step 36 value 0.000294440.
+DEAL:SIMPLEX::dim=3 degree=1 Type=Possion   : Convergence step 7 value 0.000281643.
+DEAL:SIMPLEX::dim=3 degree=1 Type=Helmholtz : Convergence step 7 value 0.000268188.
+DEAL:SIMPLEX::dim=3 degree=2 Type=Possion   : Convergence step 20 value 0.000146315.
+DEAL:SIMPLEX::dim=3 degree=2 Type=Helmholtz : Convergence step 20 value 0.000137731.
+DEAL:WEDGE  ::dim=3 degree=1 Type=Possion   : Convergence step 7 value 0.000224296.
+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.

diff --git a/tests/simplex/simplex_grids.h b/tests/simplex/simplex_grids.h
index babbb358e641a1a30d9cff01f5ea45d487cc5c01..a3bc82af8912674cc351cd2ec742858b92c60567 100644 (file)
--- a/tests/simplex/simplex_grids.h
+++ b/tests/simplex/simplex_grids.h
@@ -102,6 +102,33 @@ namespace dealii
       tria.create_triangulation(vertices, cells, SubCellData());
     }
 
+
+
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_cube_with_wedges(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_wedges(
+            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(