Extend parallel program to simplices (#138)

fdrmrc · Sep 29, 2024 · da242ce · da242ce
1 parent 274b677
commit da242ce
Showing 1 changed file with 45 additions and 10 deletions.
diff --git a/examples/diffusion_reaction.cc b/examples/diffusion_reaction.cc
@@ -17,13 +17,17 @@
 // where Omega = [0,1]^3
 //
 
+#define HEX
+
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/distributed/fully_distributed_tria.h>
 #include <deal.II/distributed/tria.h>
 
+#include <deal.II/fe/mapping_fe.h>
+
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 
@@ -307,9 +311,16 @@ class DiffusionReactionProblem
   void
   output_results() const;
 
-  const MPI_Comm                                 comm;
-  const unsigned int                             n_ranks;
-  FE_DGQ<dim>                                    fe_dg;
+  const MPI_Comm comm;
+#ifdef HEX
+  MappingQ1<dim> mapping;
+  FE_DGQ<dim>    fe_dg;
+#else
+  MappingFE<dim>     mapping;
+  FE_SimplexDGP<dim> fe_dg;
+#endif
+  const unsigned int n_ranks;
+
   const unsigned int                             n_local_agglomerates;
   const double                                   reaction_coefficient;
   parallel::fullydistributed::Triangulation<dim> tria_pft;
@@ -340,8 +351,13 @@ DiffusionReactionProblem<dim>::DiffusionReactionProblem(
   const double       reaction_coefficient,
   const MPI_Comm     communicator)
   : comm(communicator)
-  , n_ranks(Utilities::MPI::n_mpi_processes(comm))
+#ifdef HEX
+  , mapping()
+#else
+  , mapping(FE_SimplexP<dim>{1})
+#endif
   , fe_dg(degree)
+  , n_ranks(Utilities::MPI::n_mpi_processes(comm))
   , n_local_agglomerates(n_local_agglomerates)
   , reaction_coefficient(reaction_coefficient)
   , tria_pft(comm)
@@ -359,8 +375,15 @@ DiffusionReactionProblem<dim>::make_fine_grid(
   const unsigned int n_global_refinements)
 {
   Triangulation<dim> tria;
+#ifdef HEX
   GridGenerator::hyper_cube(tria);
   tria.refine_global(n_global_refinements);
+#else
+  Triangulation<dim> tria_hex;
+  GridGenerator::hyper_cube(tria_hex, 0., 1.);
+  tria_hex.refine_global(n_global_refinements - 1);
+  GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria);
+#endif
 
   // Partition serial triangulation:
   GridTools::partition_triangulation(n_ranks, tria);
@@ -414,13 +437,21 @@ DiffusionReactionProblem<dim>::assemble_system()
 {
   constraints.close();
 
-  const unsigned int quadrature_degree      = 2 * fe_dg.get_degree() + 1;
-  const unsigned int face_quadrature_degree = 2 * fe_dg.get_degree() + 1;
+  const unsigned int quadrature_degree      = fe_dg.get_degree() + 1;
+  const unsigned int face_quadrature_degree = fe_dg.get_degree() + 1;
+#ifdef HEX
   ah->initialize_fe_values(QGauss<dim>(quadrature_degree),
-                           update_values | update_gradients |
-                             update_JxW_values | update_quadrature_points,
-                           QGauss<dim - 1>(face_quadrature_degree),
-                           update_JxW_values);
+                           update_gradients | update_JxW_values |
+                             update_quadrature_points | update_JxW_values |
+                             update_values,
+                           QGauss<dim - 1>(face_quadrature_degree));
+#else
+  ah->initialize_fe_values(QGaussSimplex<dim>(quadrature_degree),
+                           update_gradients | update_JxW_values |
+                             update_quadrature_points | update_JxW_values |
+                             update_values,
+                           QGaussSimplex<dim - 1>(face_quadrature_degree));
+#endif
 
   ah->distribute_agglomerated_dofs(fe_dg);
 
@@ -834,7 +865,11 @@ main(int argc, char *argv[])
   // number of agglomerates in each local subdomain
   const unsigned int n_local_agglomerates = 10;
   const double       reaction_coefficient = .5;
+#ifdef HEX
   for (unsigned int degree : {1, 2, 3, 4})
+#else
+  for (unsigned int degree : {1, 2, 3}) // degree > 4 not implemented
+#endif
     {
       DiffusionReactionProblem<dim> problem(n_local_agglomerates,
                                             degree,