From cc15bcfe74e4795b94d09eafff4a5ac436fde636 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 2 Dec 2020 20:42:40 +0100 Subject: [PATCH] Add a poisson test on wedge grid --- tests/simplex/poisson_01.cc | 120 +++++++++++++++++- ...ilinos=true.with_simplex_support=on.output | 10 +- ...ilinos=true.with_simplex_support=on.output | 10 +- tests/simplex/simplex_grids.h | 82 ++++++++++++ 4 files changed, 215 insertions(+), 7 deletions(-) diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc index e020592b7d..e1d36fe5a5 100644 --- a/tests/simplex/poisson_01.cc +++ b/tests/simplex/poisson_01.cc @@ -56,6 +56,8 @@ #include "../tests.h" +#include "simplex_grids.h" + using namespace dealii; template @@ -95,7 +97,8 @@ test(const Triangulation &tria, const Quadrature & quad, const Quadrature & face_quad, const Mapping & mapping, - const double r_boundary) + const double r_boundary, + const bool do_use_fe_face_values = true) { std::string label = (dynamic_cast *>( @@ -188,8 +191,9 @@ test(const Triangulation &tria, std::shared_ptr> fe_face_values; - fe_face_values.reset( - new FEFaceValues(mapping, fe, face_quad, flag)); + if (do_use_fe_face_values) + fe_face_values.reset( + new FEFaceValues(mapping, fe, face_quad, flag)); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quad.size(); @@ -258,6 +262,19 @@ test(const Triangulation &tria, hex_mesh &= (cell->n_vertices() == GeometryInfo::vertices_per_cell); deallog << std::endl; + + if (false) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + + data_out.build_patches(mapping); + + std::ofstream output("result.vtk"); + data_out.write_vtk(output); + } } template @@ -389,6 +406,93 @@ test_hex(const MPI_Comm &comm, const Parameters ¶ms) test(tria, fe, quad, quad_face, mapping, params.p2[0]); } +template +void +test_wedge(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_wedges( + *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_WedgeP fe(params.degree); + + Simplex::QGaussWedge quad(params.degree + 1); + + Quadrature face_quad; // not needed + + Simplex::FE_WedgeP 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) { @@ -449,5 +553,15 @@ main(int argc, char **argv) params.p2 = Point<3>(2.1, 1, 1); test_hex(comm, params); } + + // test WEDGE + { + deallog << "Solve problem on WEDGE mesh:" << std::endl; + + params.file_name_out = "mesh-wedge"; + params.p1 = Point<3>(2.2, 0, 0); + params.p2 = Point<3>(3.2, 1, 1); + test_wedge(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 a45012d350..c89f970280 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 @@ -2,7 +2,7 @@ DEAL::Solve problem on TRI mesh: DEAL:: on parallel::fullydistributed::Triangulation DEAL:cg::Starting value 0.245798 -DEAL:cg::Convergence step 114 value 6.95996e-13 +DEAL:cg::Convergence step 114 value 6.62438e-13 DEAL:: with 114 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on QUAD mesh: @@ -14,7 +14,7 @@ DEAL:: DEAL::Solve problem on TET mesh: DEAL:: on parallel::fullydistributed::Triangulation DEAL:cg::Starting value 0.0607616 -DEAL:cg::Convergence step 156 value 9.07635e-13 +DEAL:cg::Convergence step 156 value 9.07675e-13 DEAL:: with 156 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on HEX mesh: @@ -23,3 +23,9 @@ DEAL:cg::Starting value 0.0573704 DEAL:cg::Convergence step 134 value 8.23917e-13 DEAL:: with 134 CG iterations needed to obtain convergence DEAL:: +DEAL::Solve problem on WEDGE mesh: +DEAL:: on parallel::fullydistributed::Triangulation +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:: 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 0e41dea463..5631eb08ac 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 @@ -2,7 +2,7 @@ DEAL::Solve problem on TRI mesh: DEAL:: on parallel::fullydistributed::Triangulation DEAL:cg::Starting value 0.245798 -DEAL:cg::Convergence step 114 value 6.45793e-13 +DEAL:cg::Convergence step 114 value 6.64515e-13 DEAL:: with 114 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on QUAD mesh: @@ -14,7 +14,7 @@ DEAL:: DEAL::Solve problem on TET mesh: DEAL:: on parallel::fullydistributed::Triangulation DEAL:cg::Starting value 0.0607616 -DEAL:cg::Convergence step 156 value 9.07640e-13 +DEAL:cg::Convergence step 156 value 9.07649e-13 DEAL:: with 156 CG iterations needed to obtain convergence DEAL:: DEAL::Solve problem on HEX mesh: @@ -23,3 +23,9 @@ DEAL:cg::Starting value 0.0573704 DEAL:cg::Convergence step 134 value 8.23917e-13 DEAL:: with 134 CG iterations needed to obtain convergence DEAL:: +DEAL::Solve problem on WEDGE mesh: +DEAL:: on parallel::fullydistributed::Triangulation +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:: diff --git a/tests/simplex/simplex_grids.h b/tests/simplex/simplex_grids.h index 0164b2d674..babbb358e6 100644 --- a/tests/simplex/simplex_grids.h +++ b/tests/simplex/simplex_grids.h @@ -20,6 +20,88 @@ namespace dealii { namespace GridGenerator { + template + void + subdivided_hyper_rectangle_with_wedges( + 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)); + + // 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}}; + + + // TRI cell 0 + { + CellData tri; + tri.vertices = { + quad[0], quad[1], quad[2], quad[4], quad[5], quad[6]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = { + quad[3], quad[2], quad[1], quad[7], quad[6], quad[5]}; + 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