From d5164d971b422a90ccaf3c6537dae308796224e2 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 2 Dec 2020 22:42:37 +0100 Subject: [PATCH] Add Poisson/Helmholtz test: MatrixFree + wedge mesh --- tests/simplex/matrix_free_01.cc | 77 ++++++++++++++----- ...rix_free_01.with_simplex_support=on.output | 20 +++-- tests/simplex/simplex_grids.h | 27 +++++++ 3 files changed, 97 insertions(+), 27 deletions(-) diff --git a/tests/simplex/matrix_free_01.cc b/tests/simplex/matrix_free_01.cc index 94e0f4031d..8f0df19794 100644 --- 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 void -test(const unsigned int degree, const bool do_helmholtz) +test(const unsigned int v, const unsigned int degree, const bool do_helmholtz) { Triangulation tria; - GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 16 : 8); - - Simplex::FE_P fe(degree); - Simplex::QGauss quad(degree + 1); - MappingFE mapping(Simplex::FE_P(1)); + std::shared_ptr> fe; + std::shared_ptr> quad; + std::shared_ptr> fe_mapping; + + if (v == 0) + { + GridGenerator::subdivided_hyper_cube_with_simplices(tria, + dim == 2 ? 16 : 8); + fe = std::make_shared>(degree); + quad = std::make_shared>(degree + 1); + fe_mapping = std::make_shared>(1); + } + else if (v == 1) + { + GridGenerator::subdivided_hyper_cube_with_wedges(tria, dim == 2 ? 16 : 8); + fe = std::make_shared>(degree); + quad = std::make_shared>(degree + 1); + fe_mapping = std::make_shared>(1); + } + else + Assert(false, ExcNotImplemented()); + + MappingFE mapping(*fe_mapping); DoFHandler dof_handler(tria); - dof_handler.distribute_dofs(fe); + dof_handler.distribute_dofs(*fe); AffineConstraints constraints; +#if false VectorTools::interpolate_boundary_values( mapping, dof_handler, 0, Functions::ZeroFunction(), 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 matrix_free; matrix_free.reinit( - mapping, dof_handler, constraints, quad, additional_data); + mapping, dof_handler, constraints, *quad, additional_data); PoissonOperator 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 fe_values(mapping, fe, quad, flags); + FEValues fe_values(mapping, *fe, *quad, flags); FullMatrix cell_matrix; Vector 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 ebff9c041c..5bcc1f1895 100644 --- 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 babbb358e6..a3bc82af89 100644 --- 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 + void + subdivided_hyper_cube_with_wedges(Triangulation &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 void subdivided_hyper_rectangle_with_simplices_mix( -- 2.39.5