From 766ae37c2086cc9d55b9e986ef5a40dd0943be3c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 11 Dec 2020 10:43:50 +0100 Subject: [PATCH] Add a Poisson test on pyramid grid --- tests/simplex/poisson_01.cc | 118 +++++++++++++++++ ...ilinos=true.with_simplex_support=on.output | 6 + ...ilinos=true.with_simplex_support=on.output | 6 + tests/simplex/simplex_grids.h | 121 ++++++++++++++++++ 4 files changed, 251 insertions(+) diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc index e1d36fe5a5..44e9f05800 100644 --- a/tests/simplex/poisson_01.cc +++ b/tests/simplex/poisson_01.cc @@ -49,6 +49,7 @@ #include #include +#include #include #include @@ -250,6 +251,8 @@ test(const Triangulation &tria, SolverCG solver(solver_control); solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + // deallog << solution.l2_norm() << " " << system_rhs.l2_norm() << std::endl; + deallog << " with " << solver_control.last_step() << " CG iterations needed to obtain convergence" << std::endl; @@ -263,8 +266,22 @@ test(const Triangulation &tria, deallog << std::endl; + if (false) { + Vector difference(tria.n_active_cells()); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Functions::ZeroFunction(), + difference, + quad, + VectorTools::NormType::L2_norm); + + deallog << VectorTools::compute_global_error( + tria, difference, VectorTools::NormType::L2_norm) + << std::endl; DataOut data_out; data_out.attach_dof_handler(dof_handler); @@ -493,6 +510,93 @@ test_wedge(const MPI_Comm &comm, const Parameters ¶ms) test(*tria, fe, quad, face_quad, mapping, params.p2[0], false); } +template +void +test_pyramid(const MPI_Comm &comm, const Parameters ¶ms) +{ + const unsigned int tria_type = 2; + + // 1) Create triangulation... + Triangulation *tria; + + // a) serial triangulation + Triangulation tr_1; + + // b) shared triangulation (with artificial cells) + parallel::shared::Triangulation tr_2( + MPI_COMM_WORLD, + ::Triangulation::none, + true, + parallel::shared::Triangulation::partition_custom_signal); + + tr_2.signals.create.connect([&]() { + GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm), + tr_2); + }); + + // c) distributed triangulation + parallel::fullydistributed::Triangulation tr_3(comm); + + + // ... choose the right triangulation + if (tria_type == 0 || tria_type == 2) + tria = &tr_1; + else if (tria_type == 1) + tria = &tr_2; + + // ... create triangulation + if (params.use_grid_generator) + { + // ...via Simplex::GridGenerator + GridGenerator::subdivided_hyper_rectangle_with_pyramids( + *tria, params.repetitions, params.p1, params.p2, false); + } + else + { + // ...via GridIn + GridIn grid_in; + grid_in.attach_triangulation(*tria); + std::ifstream input_file(params.file_name_in); + grid_in.read_ucd(input_file); + // std::ifstream input_file("test_tet_geometry.unv"); + // grid_in.read_unv(input_file); + } + + // ... partition serial triangulation and create distributed triangulation + if (tria_type == 0 || tria_type == 2) + { + GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm), + tr_1); + + auto construction_data = TriangulationDescription::Utilities:: + create_description_from_triangulation(tr_1, comm); + + tr_3.create_triangulation(construction_data); + + tria = &tr_3; + } + + // 2) Output generated triangulation via GridOut + GridOut grid_out; + std::ofstream out(params.file_name_out + "." + + std::to_string(Utilities::MPI::this_mpi_process(comm)) + + ".vtk"); + grid_out.write_vtk(*tria, out); + + // 3) Select components + Simplex::FE_PyramidP fe(params.degree); + + Simplex::QGaussPyramid quad(params.degree + 1); + + Quadrature face_quad; // not needed + + Simplex::FE_PyramidP fe_mapping(1); + MappingFE mapping(fe_mapping); + + // 4) Perform test (independent of mesh type) + test(*tria, fe, quad, face_quad, mapping, params.p2[0], false); +} + int main(int argc, char **argv) { @@ -563,5 +667,19 @@ main(int argc, char **argv) params.p2 = Point<3>(3.2, 1, 1); test_wedge(comm, params); } + + // test PYRAMID + { + deallog << "Solve problem on PYRAMID mesh:" << std::endl; + + params.file_name_out = "mesh-pyramid"; + params.repetitions = std::vector{10, 10, 10}; + params.p1 = Point<3>(3.3, 0, 0); + params.p2 = Point<3>(4.3, 1, 1); + + params.degree = 1; + + test_pyramid(comm, params); + } } } diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output index c89f970280..f239fef1ac 100644 --- a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output +++ b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output @@ -29,3 +29,9 @@ DEAL:cg::Starting value 0.0132550 DEAL:cg::Convergence step 186 value 8.72381e-13 DEAL:: with 186 CG iterations needed to obtain convergence DEAL:: +DEAL::Solve problem on PYRAMID mesh: +DEAL:: on parallel::fullydistributed::Triangulation +DEAL:cg::Starting value 0.0230669 +DEAL:cg::Convergence step 79 value 7.70260e-13 +DEAL:: with 79 CG iterations needed to obtain convergence +DEAL:: diff --git a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output index 5631eb08ac..01f8dc29c4 100644 --- a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output +++ b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output @@ -29,3 +29,9 @@ DEAL:cg::Starting value 0.0132550 DEAL:cg::Convergence step 186 value 8.72414e-13 DEAL:: with 186 CG iterations needed to obtain convergence DEAL:: +DEAL::Solve problem on PYRAMID mesh: +DEAL:: on parallel::fullydistributed::Triangulation +DEAL:cg::Starting value 0.0230669 +DEAL:cg::Convergence step 79 value 7.71392e-13 +DEAL:: with 79 CG iterations needed to obtain convergence +DEAL:: diff --git a/tests/simplex/simplex_grids.h b/tests/simplex/simplex_grids.h index a3bc82af89..7c91600310 100644 --- a/tests/simplex/simplex_grids.h +++ b/tests/simplex/simplex_grids.h @@ -129,6 +129,127 @@ namespace dealii + template + void + subdivided_hyper_rectangle_with_pyramids( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false) + { + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 3) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[2]); + + // create vertices + for (unsigned int k = 0; k <= repetitions[2]; ++k) + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back(Point(p1[0] + dx[0] * i, + p1[1] + dx[1] * j, + p1[2] + dx[2] * k)); + for (unsigned int k = 0; k < repetitions[2]; ++k) + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + vertices.push_back(Point(p1[0] + dx[0] * (i + 0.5), + p1[1] + dx[1] * (j + 0.5), + p1[2] + dx[2] * (k + 0.5))); + + // create cells + for (unsigned int k = 0; k < repetitions[2]; ++k) + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference HEX cell + std::array quad{ + {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (repetitions[2] + 1) * (repetitions[1] + 1) * + (repetitions[0] + 1) + + (repetitions[1] * repetitions[0] * k + + repetitions[0] * j + i)}}; + + + // TRI cell 0 + { + CellData tri; + tri.vertices = { + quad[0], quad[2], quad[4], quad[6], quad[8]}; + cells.push_back(tri); + } + // TRI cell 1 + { + CellData tri; + tri.vertices = { + quad[2], quad[3], quad[6], quad[7], quad[8]}; + cells.push_back(tri); + } + // TRI cell 2 + { + CellData tri; + tri.vertices = { + quad[6], quad[7], quad[4], quad[5], quad[8]}; + cells.push_back(tri); + } + // TRI cell 3 + { + CellData tri; + tri.vertices = { + quad[1], quad[5], quad[3], quad[7], quad[8]}; + cells.push_back(tri); + } + // TRI cell 4 + { + CellData tri; + tri.vertices = { + quad[0], quad[1], quad[2], quad[3], quad[8]}; + cells.push_back(tri); + } + // TRI cell 5 + { + CellData tri; + tri.vertices = { + quad[0], quad[4], quad[1], quad[5], quad[8]}; + cells.push_back(tri); + } + } + } + else + { + AssertThrow(colorize == false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + } + + + template void subdivided_hyper_rectangle_with_simplices_mix( -- 2.39.5