]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a Poisson test on pyramid grid 11317/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 11 Dec 2020 09:43:50 +0000 (10:43 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 11 Dec 2020 21:21:06 +0000 (22:21 +0100)
tests/simplex/poisson_01.cc
tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output
tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output
tests/simplex/simplex_grids.h

index e1d36fe5a5e9c9164126aec37835312fb39b9dc5..44e9f0580049d4f5cded51d3bd6e99f47c53a889 100644 (file)
@@ -49,6 +49,7 @@
 #include <deal.II/lac/trilinos_sparsity_pattern.h>
 
 #include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
 
 #include <deal.II/simplex/fe_lib.h>
 #include <deal.II/simplex/grid_generator.h>
@@ -250,6 +251,8 @@ test(const Triangulation<dim, spacedim> &tria,
   SolverCG<VectorType> 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<dim, spacedim> &tria,
 
   deallog << std::endl;
 
+
   if (false)
     {
+      Vector<double> difference(tria.n_active_cells());
+
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        Functions::ZeroFunction<dim>(),
+                                        difference,
+                                        quad,
+                                        VectorTools::NormType::L2_norm);
+
+      deallog << VectorTools::compute_global_error(
+                   tria, difference, VectorTools::NormType::L2_norm)
+              << std::endl;
       DataOut<dim> data_out;
 
       data_out.attach_dof_handler(dof_handler);
@@ -493,6 +510,93 @@ test_wedge(const MPI_Comm &comm, const Parameters<dim> &params)
   test(*tria, fe, quad, face_quad, mapping, params.p2[0], false);
 }
 
+template <int dim, int spacedim = dim>
+void
+test_pyramid(const MPI_Comm &comm, const Parameters<dim> &params)
+{
+  const unsigned int tria_type = 2;
+
+  // 1) Create triangulation...
+  Triangulation<dim, spacedim> *tria;
+
+  // a) serial triangulation
+  Triangulation<dim, spacedim> tr_1;
+
+  // b) shared triangulation (with artificial cells)
+  parallel::shared::Triangulation<dim> tr_2(
+    MPI_COMM_WORLD,
+    ::Triangulation<dim>::none,
+    true,
+    parallel::shared::Triangulation<dim>::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<dim> 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<dim, spacedim> 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<dim> fe(params.degree);
+
+  Simplex::QGaussPyramid<dim> quad(params.degree + 1);
+
+  Quadrature<dim - 1> face_quad; // not needed
+
+  Simplex::FE_PyramidP<dim> fe_mapping(1);
+  MappingFE<dim>            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<unsigned int>{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);
+    }
   }
 }
index c89f970280e237f9721d4fde41ba87c359a50f3d..f239fef1ac4483ae415f7cc8238837ee4b5a4e4e 100644 (file)
@@ -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::
index 5631eb08ac5888e6832b54dc4b52d78bed5bbee4..01f8dc29c4cf0541a69069e856f7c42f0fccc385 100644 (file)
@@ -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::
index a3bc82af8912674cc351cd2ec742858b92c60567..7c916003108193aedb7f07c36300ee929458191c 100644 (file)
@@ -129,6 +129,127 @@ namespace dealii
 
 
 
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle_with_pyramids(
+      Triangulation<dim, spacedim> &   tria,
+      const std::vector<unsigned int> &repetitions,
+      const Point<dim> &               p1,
+      const Point<dim> &               p2,
+      const bool                       colorize = false)
+    {
+      AssertDimension(dim, spacedim);
+
+      AssertThrow(colorize == false, ExcNotImplemented());
+
+      std::vector<Point<spacedim>> vertices;
+      std::vector<CellData<dim>>   cells;
+
+      if (dim == 3)
+        {
+          // determine cell sizes
+          const Point<dim> 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<spacedim>(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<spacedim>(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<unsigned int, 9> 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<dim> tri;
+                    tri.vertices = {
+                      quad[0], quad[2], quad[4], quad[6], quad[8]};
+                    cells.push_back(tri);
+                  }
+                  // TRI cell 1
+                  {
+                    CellData<dim> tri;
+                    tri.vertices = {
+                      quad[2], quad[3], quad[6], quad[7], quad[8]};
+                    cells.push_back(tri);
+                  }
+                  // TRI cell 2
+                  {
+                    CellData<dim> tri;
+                    tri.vertices = {
+                      quad[6], quad[7], quad[4], quad[5], quad[8]};
+                    cells.push_back(tri);
+                  }
+                  // TRI cell 3
+                  {
+                    CellData<dim> tri;
+                    tri.vertices = {
+                      quad[1], quad[5], quad[3], quad[7], quad[8]};
+                    cells.push_back(tri);
+                  }
+                  // TRI cell 4
+                  {
+                    CellData<dim> tri;
+                    tri.vertices = {
+                      quad[0], quad[1], quad[2], quad[3], quad[8]};
+                    cells.push_back(tri);
+                  }
+                  // TRI cell 5
+                  {
+                    CellData<dim> 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 <int dim, int spacedim>
     void
     subdivided_hyper_rectangle_with_simplices_mix(

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.