]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a poisson test on wedge grid 11269/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 2 Dec 2020 19:42:40 +0000 (20:42 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 2 Dec 2020 19:42:40 +0000 (20:42 +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 e020592b7df493bbc234e80e73565f74ad3f3f83..e1d36fe5a5e9c9164126aec37835312fb39b9dc5 100644 (file)
@@ -56,6 +56,8 @@
 
 #include "../tests.h"
 
+#include "simplex_grids.h"
+
 using namespace dealii;
 
 template <int dim>
@@ -95,7 +97,8 @@ test(const Triangulation<dim, spacedim> &tria,
      const Quadrature<dim> &             quad,
      const Quadrature<dim - 1> &         face_quad,
      const Mapping<dim, spacedim> &      mapping,
-     const double                        r_boundary)
+     const double                        r_boundary,
+     const bool                          do_use_fe_face_values = true)
 {
   std::string label =
     (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
@@ -188,8 +191,9 @@ test(const Triangulation<dim, spacedim> &tria,
 
   std::shared_ptr<FEFaceValues<dim, spacedim>> fe_face_values;
 
-  fe_face_values.reset(
-    new FEFaceValues<dim, spacedim>(mapping, fe, face_quad, flag));
+  if (do_use_fe_face_values)
+    fe_face_values.reset(
+      new FEFaceValues<dim, spacedim>(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<dim, spacedim> &tria,
     hex_mesh &= (cell->n_vertices() == GeometryInfo<dim>::vertices_per_cell);
 
   deallog << std::endl;
+
+  if (false)
+    {
+      DataOut<dim> 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 <int dim, int spacedim = dim>
@@ -389,6 +406,93 @@ test_hex(const MPI_Comm &comm, const Parameters<dim> &params)
   test(tria, fe, quad, quad_face, mapping, params.p2[0]);
 }
 
+template <int dim, int spacedim = dim>
+void
+test_wedge(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_wedges(
+        *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_WedgeP<dim> fe(params.degree);
+
+  Simplex::QGaussWedge<dim> quad(params.degree + 1);
+
+  Quadrature<dim - 1> face_quad; // not needed
+
+  Simplex::FE_WedgeP<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)
 {
@@ -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);
+    }
   }
 }
index a45012d350738df6f787053c837477e36c307dea..c89f970280e237f9721d4fde41ba87c359a50f3d 100644 (file)
@@ -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::
index 0e41dea463479b4e4b4e133c6154b8d41f56bc79..5631eb08ac5888e6832b54dc4b52d78bed5bbee4 100644 (file)
@@ -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::
index 0164b2d6748c7bfe01939bfa4fc376dab86e661a..babbb358e641a1a30d9cff01f5ea45d487cc5c01 100644 (file)
@@ -20,6 +20,88 @@ namespace dealii
 {
   namespace GridGenerator
   {
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle_with_wedges(
+      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));
+
+          // 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, 8> 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<dim> tri;
+                    tri.vertices = {
+                      quad[0], quad[1], quad[2], quad[4], quad[5], quad[6]};
+                    cells.push_back(tri);
+                  }
+
+                  // TRI cell 1
+                  {
+                    CellData<dim> 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 <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.