From 26214d61a1f0386722246bfef2265fe6b363ef0e Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 9 May 2021 12:02:04 +0200 Subject: [PATCH] Add simplex/mixed mesh examples to simplex module page --- doc/doxygen/headers/simplex.h | 20 ++ doc/doxygen/options.dox.in | 3 +- examples/CMakeLists.txt | 18 ++ examples/doxygen/step_3_mixed.cc | 460 +++++++++++++++++++++++++++++ examples/doxygen/step_3_simplex.cc | 432 +++++++++++++++++++++++++++ 5 files changed, 932 insertions(+), 1 deletion(-) create mode 100644 examples/doxygen/step_3_mixed.cc create mode 100644 examples/doxygen/step_3_simplex.cc diff --git a/doc/doxygen/headers/simplex.h b/doc/doxygen/headers/simplex.h index b6a7be3b73..6cbcd744af 100644 --- a/doc/doxygen/headers/simplex.h +++ b/doc/doxygen/headers/simplex.h @@ -19,6 +19,26 @@ * * @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 diff --git a/doc/doxygen/options.dox.in b/doc/doxygen/options.dox.in index 56a50b7c6a..af3962aaa9 100644 --- a/doc/doxygen/options.dox.in +++ b/doc/doxygen/options.dox.in @@ -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 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 824abea40d..0a85323f31 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -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 index 0000000000..0d61428d75 --- /dev/null +++ b/examples/doxygen/step_3_mixed.cc @@ -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. + * + * --------------------------------------------------------------------- + * + *
+ * + * + * 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. + * + * + * + * + *

Introduction

+ * + *

Motivation

+ * + * 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. + * + * + *

Working with mixed meshes

+ * + * + * 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. + * + * + * 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. + * + * + * 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. + * + * + *

Mesh generation

+ * + * 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 + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +// Include files, as added in step-3simplex: +#include + +#include +#include + +#include + +// 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 +#include +#include +#include + +using namespace dealii; + +// @sect3{The Step3 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 system_matrix; + + Vector solution; + Vector 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 cell_matrix; + Vector cell_rhs; + std::vector 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 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> 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 main function} +// +// Nothing has changed here. +int main() +{ + deallog.depth_console(2); + + Step3 laplace_problem; + laplace_problem.run(); + + return 0; +} + +/** + *

Results

+ * + * The following figures show the mesh and the result obtained by executing this + * program: + * + * + * + * + * + * + *
+ * @image html step_3_mixed_0.png + * + * @image html step_3_mixed_1.png + *
+ * + * Not surprisingly, the result looks as expected. + * + * + *

Possibilities for extensions

+ * + * 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. + * + * + *

Pure hypercube and simplex meshes

+ * + * 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. + * + * + *

3D meshes

+ * + * 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 + * + *

Parallelization, face integrals, discontinuous Galerkin methods, and + * matrix-free operator evaluation

+ * + * 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 index 0000000000..c695fe654c --- /dev/null +++ b/examples/doxygen/step_3_simplex.cc @@ -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. + * + * --------------------------------------------------------------------- + * + *
+ * + * + * 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. + * + * + * + * + *

Introduction

+ * + *

Motivation

+ * + * 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. + * + * + *

Working on simplex meshes

+ * + * 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. + * + * + *

Mesh generation

+ * + * 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 + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +// Include files that contain appropriate quadrature rules, finite elements, +// and mapping objects for simplex meshes. +#include + +#include +#include + +// The following file contains the class GridIn, which allows us to read +// external meshes. +#include + +using namespace dealii; + +// @sect3{The Step3 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 system_matrix; + + Vector solution; + Vector 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector 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 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> 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 main function} +// +// Nothing has changed here. +int main() +{ + deallog.depth_console(2); + + Step3 laplace_problem; + laplace_problem.run(); + + return 0; +} + +/* + *

Results

+ * + * The following figures show the mesh and the result obtained by executing this + * program: + * + * + * + * + * + * + *
+ * @image html step_3_simplex_0.png + * + * @image html step_3_simplex_1.png + *
+ * + * Not surprisingly, the result looks as expected. + * + * + *

Possibilities for extensions

+ * + * 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. + * + * + *

3D meshes and codim-1 meshes in 3D

+ * + * 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. + * + * + *

Mixed meshes

+ * + * 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. + * + * + *

Alternative finite elements, quadrature rules, and mapping objects

+ * + * 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). + * + * + *

Alternative mesh generation approaches

+ * + * 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 + * + * + *

hp-adaptivity

+ * + * 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. + * + * + *

Parallelization

+ * + * 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. + * + * + *

Face integrals and discontinuous Galerkin methods

+ * + * The classes FEFaceValues and FEInterfaceValues are also compatible with + * simplex meshes. + * + * + *

Matrix-free operator evaluation

+ * + * 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. + * + */ -- 2.39.5