]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add simplex/mixed mesh examples to simplex module page 12160/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 9 May 2021 10:02:04 +0000 (12:02 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 17 May 2021 19:01:39 +0000 (21:01 +0200)
doc/doxygen/headers/simplex.h
doc/doxygen/options.dox.in
examples/CMakeLists.txt
examples/doxygen/step_3_mixed.cc [new file with mode: 0644]
examples/doxygen/step_3_simplex.cc [new file with mode: 0644]

index b6a7be3b73eb68fb3dd097b6b9e27e66849584b4..6cbcd744afb0a798714f5db6b7437420bf13b3f8 100644 (file)
  *
  * @brief This module describes the experimental simplex support in deal.II.
  *
+ * Simplex and mixed meshes in deal.II are still experimental, i.e., work
+ * in progress. Large parts of the library have been ported to be able to
+ * operate on such kind of meshes. However, there are still many functions
+ * that need to be generalized. You can get a good overview of the ported
+ * functionalities by taking a look at the tests in the folder
+ * "tests/simplex". In the following, we provide two very basic examples
+ * to get started and provide some implementation details.
+ *
+ * @section simplex_reference_example_simplex Example: simplex mesh
+ *
+ * The following code shows how to work with simplex meshes:
+ *
+ * @include step_3_simplex.cc
+ *
+ * @section simplex_reference_example_mixed Example: mixed mesh
+ *
+ * The following code shows how to work with mixed meshes:
+ *
+ * @include step_3_mixed.cc
+ *
  * @section simplex_reference_cells Reference cells
  *
  * In 2D, we provide triangles and quadrilaterals with the following possible
index 56a50b7c6a1cc1ad1d6413e3cf947e0affe60340..af3962aaa987411862e32279b5deb11d4734f2b4 100644 (file)
@@ -71,7 +71,8 @@ WARN_IF_DOC_ERROR      = YES
 INPUT                  =
 RECURSIVE              = YES
 EXCLUDE_PATTERNS       = *.templates.h
-EXAMPLE_PATH           = @CMAKE_BINARY_DIR@/doc/doxygen/tutorial
+EXAMPLE_PATH           = @CMAKE_BINARY_DIR@/doc/doxygen/tutorial \
+                         @CMAKE_SOURCE_DIR@/examples/doxygen
 EXAMPLE_RECURSIVE      = NO
 IMAGE_PATH             =
 INPUT_FILTER           = ${CMAKE_SOURCE_DIR}/doc/doxygen/scripts/filter
index 824abea40d5ee24aabd77c21b3f0491eb9aa6607..0a85323f315d14b03a98d68bc46cce9bc86b92a9 100644 (file)
@@ -124,6 +124,24 @@ IF(DEAL_II_COMPONENT_EXAMPLES)
       ENDIF()
 
     ENDFOREACH()
+
+    # the same as above but for the examples folder
+    FILE(GLOB _steps
+      ${CMAKE_CURRENT_SOURCE_DIR}/doxygen/*.cc)
+    FOREACH(_step ${_steps})
+      GET_FILENAME_COMPONENT(_name ${_step} NAME_WE)
+
+        FOREACH(_build ${DEAL_II_BUILD_TYPES})
+          STRING(TOLOWER ${_build} _build_lowercase)
+          ADD_EXECUTABLE(${_name}.${_build_lowercase} ${_step})
+          DEAL_II_INSOURCE_SETUP_TARGET(${_name}.${_build_lowercase} ${_build})
+
+          SET_TARGET_PROPERTIES(${_name}.${_build_lowercase}
+            PROPERTIES
+            RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/${DEAL_II_EXECUTABLE_RELDIR}"
+            )
+        ENDFOREACH()
+    ENDFOREACH()
   ENDIF()
 
   MESSAGE(STATUS "Setting up examples - Done")
diff --git a/examples/doxygen/step_3_mixed.cc b/examples/doxygen/step_3_mixed.cc
new file mode 100644 (file)
index 0000000..0d61428
--- /dev/null
@@ -0,0 +1,460 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 1999 - 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * <br>
+ *
+ * <i>
+ * This program was contributed by Peter Munch. This work and the required
+ * generalizations of the internal data structures of deal.II form part of the
+ * project "Virtual Materials Design" funded by the Helmholtz Association of
+ * German Research Centres.
+ * </i>
+ *
+ *
+ * <a name="Intro"></a>
+ * <h1>Introduction</h1>
+ *
+ * <h3>Motivation</h3>
+ *
+ * The motivation for using simplex meshes (as done in step-3simplex) is
+ * straightforward: many freely available mesh-generation tools are very good in
+ * creating good-quality meshes in such a format, while they struggle with
+ * hex-only meshes. Hex-only meshes, on the other hand, are characterized with
+ * better numerical properties (e.g., less degrees of freedom for the same
+ * degree of accuracy and possibly better performance since the tensor-product
+ * structure can be exploited) and are, as a consequence, the natural choice for
+ * rather simple geometries and for meshes described by a coarse mesh with a few
+ * cells (like a hyperball) and obtained in their final form through iterative
+ * local refinement.
+ *
+ * Mixed meshes try to combine the best of both worlds by partitioning the
+ * geometry in parts that can be easily meshed by hypercube cells
+ * (quadrilaterals in 2D, hexahedrons in 3D) and in parts that can not be meshed
+ * easily, requiring simplices (triangles in 2D, tetrahedrons in 3D). Since one
+ * assumes that the region requiring simplices is rather small compared to the
+ * rest of the domain where more efficient and accurate methods can be applied,
+ * one can expect that the overall efficiency is hardly impacted by such an
+ * approach.
+ *
+ * One should note that in 3D, one also needs a transition region between
+ * hypercube and simplex regions. Here, one can use wedges/prisms and/or
+ * pyramids.
+ *
+ *
+ * <h3>Working with mixed meshes</h3>
+ *
+ * <i>
+ * In the following, we concentrate, for the sake of simplicity, on 2D meshes:
+ * they can only contain triangles and quadrilaterals. However, as detailed in
+ * the outlook, an extension of the presented approach to the 3D case is
+ * straightforward.
+ * </i>
+ *
+ * The complexity of working with mixed meshes in 2D results from the fact
+ * that it contains of two
+ * types of geometrical objects: quadrilaterals and triangles. How to deal with
+ * quadrilaterals, we have discussed in step-3: we selected an appropriate
+ * finite element, quadrature rule and mapping object, e.g., FE_Q, QGauss, and
+ * MappingFE (initialized with FE_Q). For simplex meshes, we selected in
+ * step-3simplex FE_SimplexP, QGaussSimplex, and MappingFE (intialized with
+ * FE_SimplexP).
+ *
+ * For mixed meshes, we need multiple finite elements, quadrature rules, and
+ * mapping objects (one set for triangles and one set for quadrilaterals) in the
+ * same program. To ease the work with the multitude of objects (in particular
+ * in 3D, we need at least four of each), you can collect the objects and group
+ * them together in hp::FECollection, hp::QCollection, and
+ * hp::MappingCollection.
+ *
+ * Just like in the context of finite elements, quadrature rules, and mapping
+ * objects, we need multiple FEValues objects: the collection of FEValues is
+ * called hp::FEValues. It returns the FEValues object needed for the current
+ * cell via the method hp::FEValues::get_present_fe_values().
+ *
+ * For hp::FEValues, to be able to select the right finite element/quadrature
+ * rule/ mapping object set, it queries the active_fe_index of the given cell
+ * during hp::FEValues::reinit(). The indices have to be set - as shown below -
+ * before calling DoFHandler::distribute_dofs() by the user.
+ *
+ * <i>
+ * The namespace name of hp::FECollection, hp::QCollection,
+ * hp::MappingCollection, and hp::FEValues indicates that these classes have not
+ * been written for mixed meshes in the first place, but for problems where each
+ * (hypercube) cell could have a different type of finite element assigned - in
+ * the simplest case, all cells have the same element type but different
+ * polynomial degrees p (the reason for the letter "p" in "hp"). An extension of
+ * this infrastructure to work not only on different element types but also on
+ * different geometrical objects was a natural choice. For further details on
+ * hp-methods, see step-27.
+ * </i>
+ *
+ * <h3>Mesh generation</h3>
+ *
+ * Just like in step-3simplex, we read an externally generated mesh. For this
+ * tutorial, we have created the mesh (square with width and height of one;
+ * quadrilaterals on the left half and triangles on the right half) with Gmsh
+ * with the following journal file "box_2D_mixed.geo":
+ *
+ * @code
+ * Rectangle(1) = {0.0, 0, 0, 0.5, 1, 0};
+ * Rectangle(2) = {0.5, 0, 0, 0.5, 1, 0};
+ * Recombine Surface{1};
+ * Physical Surface("All") = {1, 2};
+ * Mesh 2;
+ * Coherence Mesh;
+ * Save "box_2D_mixed.msh";
+ * @endcode
+ *
+ * The journal file can be processed by Gmsh generating the actual mesh with the
+ * ending ".msh":
+ *
+ * @code
+ * gmsh box_2D_mixed.geo
+ * @endcode
+ *
+ * We have included in the tutorial folder both the journal file and the mesh
+ * file in the event that one does not have access to Gmsh.
+ *
+ */
+
+
+// @sect3{Include files}
+
+// Include files, as used in step-3:
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+// Include files, as added in step-3simplex:
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_in.h>
+
+// Include files that we need in this tutorial to be able to deal with
+// collections of finite elements, quadrature rules, mapping objects, and
+// FEValues.
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+using namespace dealii;
+
+// @sect3{The <code>Step3</code> class}
+//
+// This is the main class of the tutorial. Since it is very similar to the
+// version from step-3 and step-3simplex, we will only point out and explain
+// the relevant differences that allow to perform simulations on mixed meshes.
+
+class Step3
+{
+public:
+  Step3();
+
+  void run();
+
+private:
+  void make_grid();
+  void setup_system();
+  void assemble_system();
+  void solve();
+  void output_results() const;
+
+  Triangulation<2> triangulation;
+
+  // As already explained, we are not working with mapping objects, finite
+  // elements, and quadrature rules directly but with collections of them.
+  const hp::MappingCollection<2> mapping;
+  const hp::FECollection<2>      fe;
+  const hp::QCollection<2>       quadrature_formula;
+
+  DoFHandler<2> dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  Vector<double> solution;
+  Vector<double> system_rhs;
+};
+
+
+// @sect4{Step3::Step3}
+//
+// In the constructor of the Step3 class, we fill the collections. Here, we
+// position the objects related to triangles in the first place (index 0) and
+// the ones related to quadrilaterals in the second place (index 1).
+Step3::Step3()
+  : mapping(MappingFE<2>(FE_SimplexP<2>(1)), MappingFE<2>(FE_Q<2>(1)))
+  , fe(FE_SimplexP<2>(2), FE_Q<2>(2))
+  , quadrature_formula(QGaussSimplex<2>(3), QGauss<2>(3))
+  , dof_handler(triangulation)
+{}
+
+
+// @sect4{Step3::make_grid}
+//
+// Read the external mesh file "box_2D_mixed.msh" as in step-3simplex.
+void Step3::make_grid()
+{
+  GridIn<2>(triangulation).read("box_2D_mixed.msh");
+
+  std::cout << "Number of active cells: " << triangulation.n_active_cells()
+            << std::endl;
+}
+
+
+// @sect4{Step3::setup_system}
+//
+// In contrast to step-3 and step-3simplex, we need here a preprocessing step
+// that assigns an active_fe_index to each cell consistently according to the
+// indices in the collections and the cell type.
+void Step3::setup_system()
+{
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->reference_cell() == ReferenceCells::Triangle)
+        cell->set_active_fe_index(0);
+      else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
+        cell->set_active_fe_index(1);
+      else
+        Assert(false, ExcNotImplemented());
+    }
+
+  dof_handler.distribute_dofs(fe);
+  std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+            << std::endl;
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+  sparsity_pattern.copy_from(dsp);
+
+  system_matrix.reinit(sparsity_pattern);
+
+  solution.reinit(dof_handler.n_dofs());
+  system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+// @sect4{Step3::assemble_system}
+//
+// The following function looks similar to the version in step-3 and
+// step-3simplex with the following two differences:
+//  - We do not work with FEValues directly but with the collection class
+//    hp::FEValues. It gives us - after it has been initialized with the current
+//    cell - a reference to the right FEValues object (constructed
+//    with the correct mapping object, finite element, and quadrature rule),
+//    which can be used as usual to compute the cell integrals.
+//  - The cell-local stiffness matrix and the right-hand-side vector have
+//    different sizes depending on the cell type (6 DoFs vs. 9 DoFs) so that
+//    they might need to be resized for each cell.
+//
+// Apart from these two changes, the code has not changed. Not even, the
+// cell integrals have been changed depending on whether one operates on
+// hypercube, simplex, or mixed meshes.
+void Step3::assemble_system()
+{
+  hp::FEValues<2> hp_fe_values(mapping,
+                               fe,
+                               quadrature_formula,
+                               update_values | update_gradients |
+                                 update_JxW_values);
+
+  FullMatrix<double>                   cell_matrix;
+  Vector<double>                       cell_rhs;
+  std::vector<types::global_dof_index> local_dof_indices;
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      hp_fe_values.reinit(cell);
+
+      const auto &fe_values = hp_fe_values.get_present_fe_values();
+
+      const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+      local_dof_indices.resize(dofs_per_cell);
+
+      cell_matrix = 0;
+      cell_rhs    = 0;
+
+      for (const unsigned int q_index : fe_values.quadrature_point_indices())
+        {
+          for (const unsigned int i : fe_values.dof_indices())
+            for (const unsigned int j : fe_values.dof_indices())
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+                 fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+                 fe_values.JxW(q_index));           // dx
+
+          for (const unsigned int i : fe_values.dof_indices())
+            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+                            1. *                                // f(x_q)
+                            fe_values.JxW(q_index));            // dx
+        }
+      cell->get_dof_indices(local_dof_indices);
+
+      for (const unsigned int i : fe_values.dof_indices())
+        for (const unsigned int j : fe_values.dof_indices())
+          system_matrix.add(local_dof_indices[i],
+                            local_dof_indices[j],
+                            cell_matrix(i, j));
+
+      for (const unsigned int i : fe_values.dof_indices())
+        system_rhs(local_dof_indices[i]) += cell_rhs(i);
+    }
+
+
+  std::map<types::global_dof_index, double> boundary_values;
+  VectorTools::interpolate_boundary_values(
+    mapping, dof_handler, 0, Functions::ZeroFunction<2>(), boundary_values);
+  MatrixTools::apply_boundary_values(boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+}
+
+
+// @sect4{Step3::solve}
+//
+// Nothing has changed here.
+void Step3::solve()
+{
+  SolverControl            solver_control(1000, 1e-12);
+  SolverCG<Vector<double>> solver(solver_control);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+}
+
+
+// @sect4{Step3::output_results}
+//
+// Nothing has changed here.
+void Step3::output_results() const
+{
+  DataOut<2> data_out;
+
+  DataOutBase::VtkFlags flags;
+  flags.write_higher_order_cells = true;
+  data_out.set_flags(flags);
+
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution, "solution");
+  data_out.build_patches(mapping, 2);
+  std::ofstream output("solution.vtk");
+  data_out.write_vtk(output);
+}
+
+
+// @sect4{Step3::run}
+//
+// Nothing has changed here.
+void Step3::run()
+{
+  make_grid();
+  setup_system();
+  assemble_system();
+  solve();
+  output_results();
+}
+
+
+// @sect3{The <code>main</code> function}
+//
+// Nothing has changed here.
+int main()
+{
+  deallog.depth_console(2);
+
+  Step3 laplace_problem;
+  laplace_problem.run();
+
+  return 0;
+}
+
+/**
+ * <h1>Results</h1>
+ *
+ * The following figures show the mesh and the result obtained by executing this
+ * program:
+ *
+ * <table align="center" class="doxtable" style="width:65%">
+ *   <tr>
+ *     <td>
+ *         @image html step_3_mixed_0.png
+ *     </td>
+ *     <td>
+ *         @image html step_3_mixed_1.png
+ *     </td>
+ *   </tr>
+ * </table>
+ *
+ * Not surprisingly, the result looks as expected.
+ *
+ *
+ * <h3>Possibilities for extensions</h3>
+ *
+ * In this tutorial, we presented how to use the deal.II simplex infrastructure
+ * to solve a simple Poisson problem on a mixed mesh in 2D. In this scope, we
+ * could only present a small section of the capabilities. In the following, we
+ * point out further capabilities briefly.
+ *
+ *
+ * <h4>Pure hypercube and simplex meshes</h4>
+ *
+ * In this tutorial, we worked on a mesh consisting both of quadrilaterals and
+ * triangles. However, the program and the underlying concepts also work if the
+ * mesh only contains either quadrilaterals or triangles. Interested users can
+ * try this out: we have provided appropriate journal files and meshes for such
+ * cases.
+ *
+ *
+ * <h4>3D meshes</h4>
+ *
+ * In 3D, meshes might also consist of wedges/prisms and pyramids. Therefore,
+ * the above introduced collections might consist of four components.
+ *
+ * For wedge/prism and pyramid cell types, following finite-element and
+ * quadrature-rule classes are available:
+ *  - wedge: FE_WedgeP, FE_WedgeDGP, QGaussWedge, MappingFE
+ *  - pyramid: FE_PyramidP, FE_PyramidDGP, QGaussPyramid, MappingFE
+ *
+ * <h4>Parallelization, face integrals, discontinuous Galerkin methods, and
+ * matrix-free operator evaluation</h4>
+ *
+ * Regarding these aspects, the same comments are valid that are described in
+ * step-3simplex for pure simplex meshes.
+ *
+ */
diff --git a/examples/doxygen/step_3_simplex.cc b/examples/doxygen/step_3_simplex.cc
new file mode 100644 (file)
index 0000000..c695fe6
--- /dev/null
@@ -0,0 +1,432 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 1999 - 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * <br>
+ *
+ * <i>
+ * This program was contributed by Peter Munch. This work and the required
+ * generalizations of the internal data structures of deal.II form part of the
+ * project "Virtual Materials Design" funded by the Helmholtz Association of
+ * German Research Centres.
+ * </i>
+ *
+ *
+ * <a name="Intro"></a>
+ * <h1>Introduction</h1>
+ *
+ * <h3>Motivation</h3>
+ *
+ * Many freely available mesh-generation tools produce meshes that consist of
+ * simplices (triangles in 2D; tetrahedra in 3D). The reason for this is that
+ * generating such kind of meshes for complex geometries is simpler than the
+ * generation of hex-only meshes. This tutorial shows how to work on such kind
+ * of meshes with the experimental simplex features in deal.II. For this
+ * purpose, we solve the Poisson problem from step-3 in 2D with a mesh only
+ * consisting of triangles.
+ *
+ *
+ * <h3>Working on simplex meshes</h3>
+ *
+ * To be able to work on simplex meshes, one has to select appropriate finite
+ * elements, quadrature rules, and mapping objects. In step-3, we used FE_Q,
+ * QGauss, and (implicitly by not specifying a mapping) MappingQ1. The
+ * equivalent classes for the first two classes in the context of simplices are
+ * FE_SimplexP and QGaussSimplex, which we will utilize here. For mapping
+ * purposes, we use the class MappingFE, which implements an isoparametric
+ * mapping. We initialize it with an FE_SimplexP object so that it can be
+ * applied on simplex meshes.
+ *
+ *
+ * <h3>Mesh generation</h3>
+ *
+ * In contrast to step-3, we do not use a function from the GridGenerator
+ * namespace, but rather read an externally generated mesh. For this tutorial,
+ * we have created the mesh (square with width and height of one) with Gmsh with
+ * the following journal file "box_2D_tri.geo":
+ *
+ * @code
+ * Rectangle(1) = {0, 0, 0, 1, 1, 0};
+ * Mesh 2;
+ * Save "box_2D_tri.msh";
+ * @endcode
+ *
+ * The journal file can be processed by Gmsh generating the actual mesh with the
+ * ending ".geo":
+ *
+ * @code
+ * gmsh box_2D_tri.geo
+ * @endcode
+ *
+ * We have included in the tutorial folder both the journal file and the mesh
+ * file in the event that one does not have access to Gmsh.
+ *
+ * The mesh can be simply read by deal.II with methods provided by the GridIn
+ * class, as shown below.
+ *
+ */
+
+
+// @sect3{Include files}
+
+// Include files, as used in step-3:
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+// Include files that contain appropriate quadrature rules, finite elements,
+// and mapping objects for simplex meshes.
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+// The following file contains the class GridIn, which allows us to read
+// external meshes.
+#include <deal.II/grid/grid_in.h>
+
+using namespace dealii;
+
+// @sect3{The <code>Step3</code> class}
+//
+// This is the main class of the tutorial. Since it is very similar to the
+// version from step-3, we will only point out and explain the relevant
+// differences that allow to perform simulations on simplex meshes.
+
+class Step3
+{
+public:
+  Step3();
+
+  void run();
+
+private:
+  void make_grid();
+  void setup_system();
+  void assemble_system();
+  void solve();
+  void output_results() const;
+
+  Triangulation<2> triangulation;
+
+  // Here, we select a mapping object, a finite element, and a quadrature rule
+  // that are compatible with simplex meshes.
+  const MappingFE<2>     mapping;
+  const FE_SimplexP<2>   fe;
+  const QGaussSimplex<2> quadrature_formula;
+
+  DoFHandler<2> dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  Vector<double> solution;
+  Vector<double> system_rhs;
+};
+
+
+// @sect4{Step3::Step3}
+//
+// In the constructor, we set the polynomial degree of the finite element and
+// the number of quadrature points. Furthermore, we initialize the MappingFE
+// object with a (linear) FE_SimplexP object so that it can work on simplex
+// meshes.
+Step3::Step3()
+  : mapping(FE_SimplexP<2>(1))
+  , fe(2)
+  , quadrature_formula(3)
+  , dof_handler(triangulation)
+{}
+
+
+// @sect4{Step3::make_grid}
+//
+// Read the external mesh file "box_2D_tri.msh" as in step-3.
+void Step3::make_grid()
+{
+  GridIn<2>(triangulation).read("box_2D_tri.msh");
+
+  std::cout << "Number of active cells: " << triangulation.n_active_cells()
+            << std::endl;
+}
+
+
+// @sect4{Step3::setup_system}
+//
+// From here on, nothing has changed. Not even, the
+// cell integrals have been changed depending on whether one operates on
+// hypercube or simplex meshes. This is astonishing and is possible due to the
+// design of the following two classes:
+//  - DoFHandler: this class stores degrees of freedom in a flexible way and
+//    allows simple access to them depending on the element type independent of
+//    the cell type.
+//  - FEValues: this class hides the details of finite element, quadrature rule,
+//    and mapping (even if the implementations are inherently different)
+//    behind a unified interface.
+void Step3::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+  std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+            << std::endl;
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+  sparsity_pattern.copy_from(dsp);
+
+  system_matrix.reinit(sparsity_pattern);
+
+  solution.reinit(dof_handler.n_dofs());
+  system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+// @sect4{Step3::assemble_system}
+//
+// Nothing has changed here.
+void Step3::assemble_system()
+{
+  FEValues<2> fe_values(mapping,
+                        fe,
+                        quadrature_formula,
+                        update_values | update_gradients | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+
+      cell_matrix = 0;
+      cell_rhs    = 0;
+
+      for (const unsigned int q_index : fe_values.quadrature_point_indices())
+        {
+          for (const unsigned int i : fe_values.dof_indices())
+            for (const unsigned int j : fe_values.dof_indices())
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+                 fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+                 fe_values.JxW(q_index));           // dx
+
+          for (const unsigned int i : fe_values.dof_indices())
+            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+                            1. *                                // f(x_q)
+                            fe_values.JxW(q_index));            // dx
+        }
+      cell->get_dof_indices(local_dof_indices);
+
+      for (const unsigned int i : fe_values.dof_indices())
+        for (const unsigned int j : fe_values.dof_indices())
+          system_matrix.add(local_dof_indices[i],
+                            local_dof_indices[j],
+                            cell_matrix(i, j));
+
+      for (const unsigned int i : fe_values.dof_indices())
+        system_rhs(local_dof_indices[i]) += cell_rhs(i);
+    }
+
+
+  std::map<types::global_dof_index, double> boundary_values;
+  VectorTools::interpolate_boundary_values(
+    mapping, dof_handler, 0, Functions::ZeroFunction<2>(), boundary_values);
+  MatrixTools::apply_boundary_values(boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+}
+
+
+// @sect4{Step3::solve}
+//
+// Nothing has changed here.
+void Step3::solve()
+{
+  SolverControl            solver_control(1000, 1e-12);
+  SolverCG<Vector<double>> solver(solver_control);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+}
+
+
+// @sect4{Step3::output_results}
+//
+// Nothing has changed here.
+void Step3::output_results() const
+{
+  DataOut<2> data_out;
+
+  DataOutBase::VtkFlags flags;
+  flags.write_higher_order_cells = true;
+  data_out.set_flags(flags);
+
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution, "solution");
+  data_out.build_patches(mapping, 2);
+  std::ofstream output("solution.vtk");
+  data_out.write_vtk(output);
+}
+
+
+// @sect4{Step3::run}
+//
+// Nothing has changed here.
+void Step3::run()
+{
+  make_grid();
+  setup_system();
+  assemble_system();
+  solve();
+  output_results();
+}
+
+
+// @sect3{The <code>main</code> function}
+//
+// Nothing has changed here.
+int main()
+{
+  deallog.depth_console(2);
+
+  Step3 laplace_problem;
+  laplace_problem.run();
+
+  return 0;
+}
+
+/*
+ * <h1>Results</h1>
+ *
+ * The following figures show the mesh and the result obtained by executing this
+ * program:
+ *
+ * <table align="center" class="doxtable" style="width:65%">
+ *   <tr>
+ *     <td>
+ *         @image html step_3_simplex_0.png
+ *     </td>
+ *     <td>
+ *         @image html step_3_simplex_1.png
+ *     </td>
+ *   </tr>
+ * </table>
+ *
+ * Not surprisingly, the result looks as expected.
+ *
+ *
+ * <h3>Possibilities for extensions</h3>
+ *
+ * In this tutorial, we presented how to use the deal.II simplex infrastructure
+ * to solve a simple Poisson problem on a simplex mesh in 2D. In this scope, we
+ * could only present a small section of the capabilities. In the following, we
+ * point out further capabilities briefly.
+ *
+ *
+ * <h4>3D meshes and codim-1 meshes in 3D</h4>
+ *
+ * An extension to 3D is quite straightforward. Both FE_SimplexP and
+ * QGaussSimplex are implemented in a dimensional-independent way so that simply
+ * replacing everywhere dim=2 with dim=3 should work out of the box.
+ *
+ * Furthermore, embedding of a 2D mesh consisting of triangles in 3D space is
+ * possible.
+ *
+ *
+ * <h4>Mixed meshes</h4>
+ *
+ * In step-3, we considered meshes only consisting of quadrilaterals. In this
+ * tutorial, we took a look at the case that the mesh only consists of
+ * triangles. In the general case (also known as mixed mesh), the mesh consists
+ * of both cell types. In 3D, meshes might even consist of more cell types, like
+ * wedges/prisms and pyramids. We consider such meshes in the tutorial
+ * step-3mixed.
+ *
+ *
+ * <h4>Alternative finite elements, quadrature rules, and mapping objects</h4>
+ *
+ * In this tutorial, we used the most basic finite-element, quadrature-rule, and
+ * mapping classes. However, more classes are compatible with simplices. The
+ * following list gives an overview of these classes:
+ * - finite elements: FE_SimplexP, FE_SimplexDGP, FE_SimplexP_Bubbles
+ * - quadrature rules: QGaussSimplex, QWitherdenVincentSimplex, QDuffy
+ * - mapping objects: MappingFE, MappingFEField
+ *
+ * It should be also pointed out that FESystems can also handle simplex finite
+ * elements which is crucial to solve vector-valued problems, as needed, e.g.,
+ * to solve elasticity and fluid problems (see also step-17).
+ *
+ *
+ * <h4>Alternative mesh generation approaches</h4>
+ *
+ * In this tutorial, we have created the mesh externally and read it with the
+ * help of GridIn. Since we believe that the main motivation to work on simplex
+ * meshes is that one has a complex geometry that can only be meshed with
+ * an external tool with simplices, deal.II does not have too many functions in
+ * the GridGenerator namespace, targeting simplex meshes. However, we would like
+ * to point out the following functions:
+ *  - GridGenerator::subdivided_hyper_cube_with_simplices() and
+ *    GridGenerator::subdivided_hyper_rectangle_with_simplices(), which fill a
+ *    hypercube and a hyperrectangle domain with simplices
+ *  - GridGenerator::convert_hypercube_to_simplex_mesh(), which converts meshes
+ *    consisting of hypercube cells to simplex meshes by replacing one
+ *    quadrilateral with 4 triangles and one hexahedron with 24 tetrahedrons
+ *
+ *
+ * <h4>hp-adaptivity</h4>
+ *
+ * Here, we considered a mesh without refinements and with all cells assigned
+ * the same type of element with the same polynomial degree. However, one is not
+ * restricted to this. For further details on hp-methods, see step-27.
+ *
+ *
+ * <h4>Parallelization</h4>
+ *
+ * To parallelize the code, one needs to replace the Triangulation object either
+ * with parallel::shared::Triangulation or
+ * parallel::fullydistributed::Triangulation and make some minor adjustments, as
+ * discussed in step-6.
+ *
+ *
+ * <h4>Face integrals and discontinuous Galerkin methods</h4>
+ *
+ * The classes FEFaceValues and FEInterfaceValues are also compatible with
+ * simplex meshes.
+ *
+ *
+ * <h4>Matrix-free operator evaluation</h4>
+ *
+ * In this tutorial, we showed a matrix-based approach. However, one could also
+ * rewrite the code using MatrixFree, FEEvaluation, and FEFaceEvaluation, which
+ * are also compatible with simplex meshes.
+ *
+ */

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.