From: Simon Sticko Date: Sun, 2 Jan 2022 10:27:50 +0000 (+0100) Subject: Add a CutFEM tutorial as step-85 X-Git-Tag: v9.4.0-rc1~549^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=18b20aac6be1827ce47b9da1d30cb9e525293bae;p=dealii.git Add a CutFEM tutorial as step-85 This step shows how the NonMatching classes: MeshClassifier, FEImmersedSurfaceValues and NonMatching::FEValues are intended to work together. This is done by solving the Poisson equation on a hypercube mesh where the domain is described by a level set function --- diff --git a/doc/doxygen/images/step-85-active-mesh.svg b/doc/doxygen/images/step-85-active-mesh.svg new file mode 100644 index 0000000000..05f00cae1a --- /dev/null +++ b/doc/doxygen/images/step-85-active-mesh.svg @@ -0,0 +1,377 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-background-mesh.svg b/doc/doxygen/images/step-85-background-mesh.svg new file mode 100644 index 0000000000..8387329a38 --- /dev/null +++ b/doc/doxygen/images/step-85-background-mesh.svg @@ -0,0 +1,411 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-ghost-faces.svg b/doc/doxygen/images/step-85-ghost-faces.svg new file mode 100644 index 0000000000..396725c343 --- /dev/null +++ b/doc/doxygen/images/step-85-ghost-faces.svg @@ -0,0 +1,239 @@ + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/doxygen/images/step-85-level-set.png b/doc/doxygen/images/step-85-level-set.png new file mode 100644 index 0000000000..68cb21ea40 Binary files /dev/null and b/doc/doxygen/images/step-85-level-set.png differ diff --git a/doc/doxygen/images/step-85-solution.png b/doc/doxygen/images/step-85-solution.png new file mode 100644 index 0000000000..3b8932a25b Binary files /dev/null and b/doc/doxygen/images/step-85-solution.png differ diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index ceb8ad7ec6..ca27997f69 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1375,6 +1375,59 @@ url = {https://arxiv.org/abs/2106.13877} } +% ------------------------------------ +% Step 85 +% ------------------------------------ + +@article{saye_2015, + title = {High-{Order} {Quadrature} {Methods} for {Implicitly} {Defined} {Surfaces} and {Volumes} in {Hyperrectangles}}, + volume = {37}, + issn = {1064-8275, 1095-7197}, + url = {http://epubs.siam.org/doi/10.1137/140966290}, + doi = {10.1137/140966290}, + language = {en}, + number = {2}, + urldate = {2016-01-27}, + journal = {SIAM Journal on Scientific Computing}, + author = {Saye, R. I.}, + month = jan, + year = {2015}, + pages = {A993--A1019} +} + +@article{burman_hansbo_2012, + title = {Fictitious domain finite element methods using cut elements: {II}. {A} stabilized {Nitsche} method}, + volume = {62}, + issn = {01689274}, + url = {http://linkinghub.elsevier.com/retrieve/pii/S0168927411000298}, + doi = {10.1016/j.apnum.2011.01.008}, + language = {en}, + number = {4}, + urldate = {2015-12-17}, + journal = {Applied Numerical Mathematics}, + author = {Burman, Erik and Hansbo, Peter}, + month = apr, + year = {2012}, + pages = {328--341} +} + +@article{cutfem_2015, + title = {{CutFEM}: {Discretizing} geometry and partial differential equations}, + volume = {104}, + issn = {00295981}, + url = {http://doi.wiley.com/10.1002/nme.4823}, + doi = {10.1002/nme.4823}, + language = {en}, + number = {7}, + urldate = {2016-01-27}, + journal = {International Journal for Numerical Methods in Engineering}, + author = {Burman, Erik and Claus, Susanne and Hansbo, Peter and Larson, Mats G. and Massing, André}, + month = nov, + year = {2015}, + pages = {472--501} +} + + % ------------------------------------ % References used elsewhere % ------------------------------------ diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 851bda1476..fed5538315 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -664,6 +664,11 @@ * Solving the fourth-order biharmonic equation using a lifting operator approach. * * + * + * step-85 + * Solving the Poisson equation using the cut finite element method. + * + * * * * @@ -1004,7 +1009,8 @@ * * * step-60, - * step-70 + * step-70, + * step-85 * * * diff --git a/examples/step-85/CMakeLists.txt b/examples/step-85/CMakeLists.txt new file mode 100644 index 0000000000..dbc0f7dfb4 --- /dev/null +++ b/examples/step-85/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-85 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-85") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 10.0.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-85/doc/builds-on b/examples/step-85/doc/builds-on new file mode 100644 index 0000000000..a22675975f --- /dev/null +++ b/examples/step-85/doc/builds-on @@ -0,0 +1 @@ +step-12 step-46 diff --git a/examples/step-85/doc/intro.dox b/examples/step-85/doc/intro.dox new file mode 100644 index 0000000000..dbd5feca4c --- /dev/null +++ b/examples/step-85/doc/intro.dox @@ -0,0 +1,251 @@ + +This program was contributed by Simon Sticko. + +The material is based upon work partially supported by +eSSENCE of e-Science and the Swedish Research Council +under grants 2014-6088 (Kreiss) and 2017-05038 (Massing). + + + +

Introduction

+ +

The Cut Finite Element Method

+ +In this example, we show how to use the cut finite element method (CutFEM) in deal.II. +For illustration, we want to solve the simplest possible problem, +so we again consider Poisson's equation: +@f{align*} + -\Delta u &= f \qquad && \text{in }\, \Omega, + \\ + u &= u_D \qquad && \text{on }\, \Gamma = \partial \Omega, +@f} +where we choose $f(x) = 4$ and $u_D(x) = 1$. +CutFEM is an immersed method. +In this context, +"immersed" means that the mesh is unfitted to the geometry of the domain, $\Omega$. +Instead, $\Omega$ floats freely on top of a uniform background mesh, $\mathcal{T}^h$. +@image html step-85-background-mesh.svg +Since we no longer use the mesh to describe the geometry of the domain, +we need some other way to represent it. +This can be done in several ways but here we assume that $\Omega$ is described by a level set function, +$\psi : \mathbb{R}^{\text{dim}} \to \mathbb{R}$ such that +@f{align*} + \Omega &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) < 0 \}, \\ + \Gamma &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) = 0 \}. +@f} +For simplicity, we choose $\Omega$ to be a unit disk, so that +@f{equation*} + \psi(x) = \| x \| - 1. +@f} +As can be seen from the figure below, +the level set function is negative for points in $\Omega$, +zero on the boundary, and positive everywhere else. +@image html step-85-level-set.png +To solve this problem, +we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$, +that completely covers the domain: +@f{equation*} + \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Omega \neq \emptyset \}. +@f} +This is usually referred to as the "active mesh". +@image html step-85-active-mesh.svg +The finite element space where we want to find our numerical solution, $u_h$, is now +@f{equation*} + V_\Omega^h = \{ v \in C(\mathcal{N}_\Omega^h) : v \in Q_p(T), \, T \in \mathcal{T}_\Omega^h \}, +@f} +where +@f{equation*} + \mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T}, +@f} +and $\overline{T}$ denotes the closure of $T$. +The set $\mathcal{N}_\Omega^h$ is sometimes referred to as the "fictitious domain". +Since $\Omega \subset \mathcal{N}_\Omega^h$, +we see that the numerical solution is defined over a slightly larger region than the analytical solution. + +In this type of immersed finite element method, +the standard way to apply boundary conditions is using Nitsche's method. +Multiplying the PDE with a test function, $v_h \in V_\Omega^h$, +and integrating by parts over $\Omega$, as usual, gives us +@f{equation*} + (\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma = (f,v)_\Omega. +@f} +Let $\gamma_D > 0$ be a scalar penalty parameter and let $h$ be some measure of the local cell size. +We now note that the following terms are consistent with the Dirichlet boundary condition: +@f{align*} + -(u_h, \partial_n v_h)_\Gamma &= -(u_D, \partial_n v_h)_\Gamma, \\ + \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma &= \left (\frac{\gamma_D}{h}u_D, v_h \right )_\Gamma. +@f} +Thus, we can add these to the weak formulation to enforce the boundary condition. +This leads to the following weak formulation: +Find $u_h \in V_\Omega^h$ such that +@f{equation*} + a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, +@f} +where +@f{align*} + a_h(u_h, v_h) &= (\nabla u_h, \nabla v_h)_\Omega + - (\partial_n u_h, v_h)_\Gamma + - (u_h, \partial_n v_h)_\Gamma + + \left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma, + \\ + L_h(v_h) &= (f,v)_\Omega + + \left (u_D, \frac{\gamma_D}{h} v_h -\partial_n v_h \right )_\Gamma. +@f} +In this formulation, there is one big difference, +compared to a standard boundary-fitted finite element method. +On each cell, +we need to integrate over the part of the domain and the part of the boundary that falls within the cell. +Thus, on each cell intersected by $\Gamma$, +we need special quadrature rules that only integrate over these parts of the cell, +that is, over $T \cap \Omega$ and $T \cap \Gamma$. +@image html immersed_quadratures.svg +Since $\Omega \cap T$ is the part of the cell that lies inside the domain, +we shall refer to the following regions +@f{align*} + \{x \in T : \psi(x) < 0 \}, \\ + \{x \in T : \psi(x) > 0 \}, \\ + \{x \in T : \psi(x) = 0 \}, +@f} +as the "inside", "outside" and the "surface region" of the cell $T$. + +The above finite element method that uses the bilinear form $a_h(\cdot, \cdot)$ +is sometimes referred to as the "naive weak formulation" +because it suffers from the so-called "small cut problem". +Depending on how $\Omega$ is located relative to $\mathcal{T}_h$, +a cut between a cell, $T \in \mathcal{T}_h$, and $\Omega$ can become arbitrarily small: +$|\Omega \cap T | \rightarrow 0$. +For Neumann boundary conditions, +the consequence is that the stiffness matrix can become arbitrarily ill-conditioned +as the cut-size approaches zero. +For a Dirichlet condition, the situation is even worse. +For any finite choice of Nitsche constant, $\gamma_D$, +the bilinear form $a_h(\cdot,\cdot)$ loses coercivity as the size of a cell cut approaches zero. +This makes the above weak formulation essentially useless +because as we refine we typically can not control how the cells intersect $\Gamma$. +One way to avoid this problem is to add a so-called ghost penalty term, $g_h$, +to the weak formulation (see e.g. @cite burman_hansbo_2012 and @cite cutfem_2015). +This leads to the stabilized cut finite element method, +which reads: Find $u_h \in V_\Omega^h$ such that +@f{equation*} + A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h, +@f} +where +@f{equation*} + A_h(u_h,v_h) = a_h(u_h,v_h) + g_h(u_h, v_h). +@f} +The point of this ghost penalty is that it makes the numerical method essentially independent +of how $\Omega$ relates to the background mesh. +In particular, $A_h$ can be shown to be continuous and coercive, +with constants that do not depend on how $\Omega$ intersects $\mathcal{T}^h$. +To define the ghost penalty, let $\mathcal{T}_\Gamma^h$ be the set of intersected cells: +@f{equation*} + \mathcal{T}_{\Gamma}^h = \{ T \in \mathcal{T}_{\Omega}^{h} : T \cap \Gamma \neq \emptyset \}, +@f} +and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in the active mesh: +@f{equation*} + \mathcal{F}_h = \{ F = \overline{T}_+ \cap \overline{T}_- : \, + T_+ \in \mathcal{T}_{\Gamma}^h, \, + T_- \in \mathcal{T}_{\Omega}^h + \}. +@f} +@image html step-85-ghost-faces.svg +The ghost penalty acts on these faces and reads +@f{equation*} + g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} g_F(u_h, v_h), +@f} +where $g_F$ is the face-wise ghost penalty: +@f{equation*} + g_F(u_h, v_h) = \gamma_A \sum_{k=0}^p \left(\frac{h_F^{2k-1}}{k!^2}[\partial_n^k u_h], [\partial_n^k v_h] \right)_F. +@f} +Here, $\gamma_A$ is a penalty parameter and $h_F$ is some measure of the face size. +We see that $g_F$ penalizes the jumps in the face-normal derivatives, $\partial_n^k$, +over $F = \overline{T}_+ \cap \overline{T}_-$. +Since we include all normal derivatives up to the polynomial degree, +we weakly force the piecewise polynomial to behave as a single polynomial over $\overline{T}_+ \cup \overline{T}_-$. +Hand-wavingly speaking, +this is the reason why we obtain a cut-independent method when we enforce $g_F(u_h, v_h) = 0$ over the faces in $\mathcal{F}_h$. +Here, we shall use a continuous space of $Q_1$-elements, +so the ghost penalty is reduced to +@f{equation*} + g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} (h_F [\partial_n u_h], [\partial_n v_h])_F. +@f} + +

Discrete Level Set Function

+A typical use case of a level set method is a problem where the domain is advected in a velocity field, +such that the domain deforms with time. +For such a problem, +one would typically solve for an approximation of the level set function, $\psi_h \in V^h$, +in a separate finite element space over the whole background mesh: +@f{equation*} + V^h = \{ v \in C(\mathcal{N}^h) : v \in Q_p(T), \, T \in \mathcal{T}^h \}, +@f} +where $\mathcal{N}^h = \bigcup_{T \in \mathcal{T}^h} \overline{T}$. +Even if we solve a much simpler problem with a stationary domain in this tutorial, +we shall, just to illustrate, still use a discrete level set function for the Poisson problem. +Technically, +this is a so-called "variational crime" because we are actually not using the bilinear form $a_h$ but instead +@f{equation*} + a_h^\star(u_h, v_h) = (\nabla u_h, \nabla v_h)_{\Omega_h} + - (\partial_n u_h, v_h)_{\Gamma_h} + \ldots +@f} +This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function: +@f{align*} + \Omega_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) < 0 \}, \\ + \Gamma_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) = 0 \}. +@f} +Using $\Omega_h$ instead of $\Omega$ in the method will give rise to a larger error in the numerical solution. +This is often referred to as the "geometrical error". +However, when the same element order, $p$, is used in $V^h$ and $V_\Omega^h$, +one can often show that the method gives the same order of convergence +as if the exact domain would have been used. +Furthermore, deal.II allows us to independently choose a more accurate geometry representation +with a higher-order level set function, compared to the function space for solving the Poisson equation. + +

The MeshClassifier Class

+Even if we have used $\mathcal{T}_\Omega^h$ to define the finite element space, +we will not create this submesh in practice. As in step-46, we shall +instead use the hp-framework. To create $V_\Omega^h$, we first add an FE_Q and an +FE_Nothing element to an hp::FECollection. We then iterate over each cell, +$T$, and depending on whether $T$ belongs to $\mathcal{T}_\Omega^h$ or not, +we set the active_fe_index to either 0 or 1. To do so, we need to +determine if a given cell is in $\mathcal{T}_\Omega^h$ or not. +For this purpose, we will use the class NonMatching::MeshClassifier. +The MeshClassifier takes the discrete level set function, +described as a (DoFHandler, Vector)-pair, as arguments to its constructor: +@code + MeshClassifier(const DoFHandler &level_set_dof_handler, + const VectorType & level_set); +@endcode +When we call the reclassify() function on an object of this class, +each active cell and face is associated with one of the values +{inside, outside, intersected} of the enum NonMatching::LocationToLevelSet. +Here, "inside" means that the level set function is negative over the whole cell +so that it lies completely inside the domain. +Analogously, "outside" means that $\psi$ is positive over the whole cell, +and "intersected" means that $\psi$ varies in sign over $T$ +so that the zero-contour of $\psi$ goes through $T$. + +| LocationToLevelSet | $\psi(x)$ for $x \in T$ | Relation to $\Omega$ | +|:------------------:|:------------------------:|:------------------------------:| +| inside | $\psi(x) < 0$ | $T \cap \Omega = T$ | +| outside | $0 < \psi(x)$ | $T \cap \Omega = \emptyset$ | +| intersected | $\psi(x)$ varies in sign | $T \cap \Gamma \neq \emptyset$ | + +Each active face is classified in the same way, according to how the sign of $\psi$ varies over the face. +MeshClassifier lets you query this information for a given cell/face via its MeshClassifier::location_to_level_set() methods: +@code + NonMatching::MeshClassifier mesh_classifier(dof_handler, level_set); + mesh_classifier.reclassify(); + + for (const auto &cell : triangulation.active_cell_iterators()) + { + NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + for (const unsigned int f : cell->face_indices()) + { + NonMatching::LocationToLevelSet face_location = + mesh_classifier.location_to_level_set(cell, f); + } + } +@endcode diff --git a/examples/step-85/doc/kind b/examples/step-85/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-85/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-85/doc/results.dox b/examples/step-85/doc/results.dox new file mode 100644 index 0000000000..b8cfe758ff --- /dev/null +++ b/examples/step-85/doc/results.dox @@ -0,0 +1,20 @@ +

Results

+ +The numerical solution for one of the refinements is shown in the below figure. +The zero-contour of the level set function is shown as a white line. +On the intersected cells, +we see that the numerical solution has a value also outside $\overline{\Omega}$. +As mentioned earlier, this extension of the solution is artificial. + +The results of the convergence study is shown in the table below. +We see that the $L^2$ error decreases as we refine and that the estimated +order of convergence, EOC, is close to 2. + +@image html step-85-solution.png + +| Cycle | Mesh size | $L^2$-Error | EOC | +|:-----:|:---------:|:-----------:|:----:| +| 0 | 0.3025 | 8.0657e-02 | - | +| 1 | 0.1513 | 1.8711e-02 | 2.11 | +| 2 | 0.0756 | 4.1624e-03 | 2.17 | +| 3 | 0.0378 | 9.3979e-04 | 2.15 | diff --git a/examples/step-85/doc/tooltip b/examples/step-85/doc/tooltip new file mode 100644 index 0000000000..927b4dc20e --- /dev/null +++ b/examples/step-85/doc/tooltip @@ -0,0 +1 @@ +Solving Poisson's equation using the cut finite element method. diff --git a/examples/step-85/step-85.cc b/examples/step-85/step-85.cc new file mode 100644 index 0000000000..a3a8e8be15 --- /dev/null +++ b/examples/step-85/step-85.cc @@ -0,0 +1,714 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2021 - 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. + * + * --------------------------------------------------------------------- + */ + +// @sect3{Include files} + +// The first include files have all been treated in previous examples. + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +// The first new header contains some common level set functions. +// For example, the spherical geometry that we use here. +#include + +// We also need 3 new headers from the NonMatching namespace. +#include +#include +#include + +// @sect3{The LaplaceSolver class Template} +// We then define the main class that solves the Laplace problem. + +namespace Step85 +{ + using namespace dealii; + + template + class LaplaceSolver + { + public: + LaplaceSolver(); + + void run(); + + private: + void make_grid(); + + void setup_discrete_level_set(); + + void distribute_dofs(); + + void initialize_matrices(); + + void assemble_system(); + + void solve(); + + void output_results() const; + + double compute_L2_error() const; + + bool face_has_ghost_penalty( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index) const; + + const unsigned int fe_degree; + + const Functions::ConstantFunction rhs_function; + const Functions::ConstantFunction boundary_condition; + + Triangulation triangulation; + + // We need two separate DoFHandlers. The first manages the DoFs for the + // discrete level set function that describes the geometry of the domain. + const FE_Q fe_level_set; + DoFHandler level_set_dof_handler; + Vector level_set; + + // The second DoFHandler manages the DoFs for the solution of the Poisson + // equation. + hp::FECollection fe_collection; + DoFHandler dof_handler; + Vector solution; + + NonMatching::MeshClassifier mesh_classifier; + + SparsityPattern sparsity_pattern; + SparseMatrix stiffness_matrix; + Vector rhs; + }; + + + + template + LaplaceSolver::LaplaceSolver() + : fe_degree(1) + , rhs_function(4.0) + , boundary_condition(1.0) + , fe_level_set(fe_degree) + , level_set_dof_handler(triangulation) + , dof_handler(triangulation) + , mesh_classifier(level_set_dof_handler, level_set) + {} + + + + // @sect3{Setting up the Background Mesh} + // We generate a background mesh with perfectly Cartesian cells. Our domain is + // a unit disc centered at the origin, so we need to make the background mesh + // a bit larger than $[-1, 1]^{\text{dim}}$ to completely cover $\Omega$. + template + void LaplaceSolver::make_grid() + { + std::cout << "Creating background mesh" << std::endl; + + GridGenerator::hyper_cube(triangulation, -1.21, 1.21); + triangulation.refine_global(2); + } + + + + // @sect3{Setting up the Discrete Level Set Function} + // The discrete level set function is defined on the whole background mesh. + // Thus, to set up the DoFHandler for the level set function, we distribute + // DoFs over all elements in $\mathcal{T}_h$. We then set up the discrete + // level set function by interpolating onto this finite element space. + template + void LaplaceSolver::setup_discrete_level_set() + { + std::cout << "Setting up discrete level set function" << std::endl; + + level_set_dof_handler.distribute_dofs(fe_level_set); + level_set.reinit(level_set_dof_handler.n_dofs()); + + const Functions::SignedDistance::Sphere signed_distance_sphere; + VectorTools::interpolate(level_set_dof_handler, + signed_distance_sphere, + level_set); + } + + + + // @sect3{Setting up the Finite Element Space} + // To set up the finite element space $V_\Omega^h$, we will use 2 different + // elements: FE_Q and FE_Nothing. For better readability we define an enum for + // the indices in the order we store them in the hp::FECollection. + enum ActiveFEIndex + { + lagrange = 0, + nothing = 1 + }; + + // We then use the MeshClassifier to check LocationToLevelSet for each cell in + // the mesh and tell the DoFHandler to use FE_Q on elements that are inside or + // intersected, and FE_Nothing on the elements that are outside. + template + void LaplaceSolver::distribute_dofs() + { + std::cout << "Distributing degrees of freedom" << std::endl; + + fe_collection.push_back(FE_Q(fe_degree)); + fe_collection.push_back(FE_Nothing()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + const NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + if (cell_location == NonMatching::LocationToLevelSet::outside) + cell->set_active_fe_index(ActiveFEIndex::nothing); + else + cell->set_active_fe_index(ActiveFEIndex::lagrange); + } + + dof_handler.distribute_dofs(fe_collection); + } + + + + // @sect3{Sparsity Pattern} + // The added ghost penalty results in a sparsity pattern similar to a DG + // method with a symmetric-interior-penalty term. Thus, we can use the + // make_flux_sparsity_pattern() function to create it. However, since the + // ghost-penalty terms only act on the faces in $\mathcal{F}_h$, we can pass + // in a lambda function that tells make_flux_sparsity_pattern() over which + // faces the flux-terms appear. This gives us a sparsity pattern with minimal + // number of entries. When passing a lambda function, + // make_flux_sparsity_pattern requires us to also pass cell and face coupling + // tables to it. If the problem was vector-valued, these tables would allow us + // to couple only some of the vector components. This is discussed in step-46. + template + void LaplaceSolver::initialize_matrices() + { + std::cout << "Initializing matrices" << std::endl; + + const auto face_has_flux_coupling = [&](const auto & cell, + const unsigned int face_index) { + return this->face_has_ghost_penalty(cell, face_index); + }; + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + + const unsigned int n_components = fe_collection.n_components(); + Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components); + Table<2, DoFTools::Coupling> face_coupling(n_components, n_components); + cell_coupling[0][0] = DoFTools::always; + face_coupling[0][0] = DoFTools::always; + + const AffineConstraints constraints; + const bool keep_constrained_dofs = true; + + DoFTools::make_flux_sparsity_pattern(dof_handler, + dsp, + constraints, + keep_constrained_dofs, + cell_coupling, + face_coupling, + numbers::invalid_subdomain_id, + face_has_flux_coupling); + sparsity_pattern.copy_from(dsp); + + stiffness_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + rhs.reinit(dof_handler.n_dofs()); + } + + + + // The following function describes which faces are part of the set + // $\mathcal{F}_h$. That is, it returns true if the face of the incoming cell + // belongs to the set $\mathcal{F}_h$. + template + bool LaplaceSolver::face_has_ghost_penalty( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index) const + { + if (cell->at_boundary(face_index)) + return false; + + const NonMatching::LocationToLevelSet cell_location = + mesh_classifier.location_to_level_set(cell); + + const NonMatching::LocationToLevelSet neighbor_location = + mesh_classifier.location_to_level_set(cell->neighbor(face_index)); + + if (cell_location == NonMatching::LocationToLevelSet::intersected && + neighbor_location != NonMatching::LocationToLevelSet::outside) + return true; + + if (neighbor_location == NonMatching::LocationToLevelSet::intersected && + cell_location != NonMatching::LocationToLevelSet::outside) + return true; + + return false; + } + + + + // @sect3{Assembling the System} + template + void LaplaceSolver::assemble_system() + { + std::cout << "Assembling" << std::endl; + + const unsigned int n_dofs_per_cell = fe_collection[0].dofs_per_cell; + FullMatrix local_stiffness(n_dofs_per_cell, n_dofs_per_cell); + Vector local_rhs(n_dofs_per_cell); + std::vector local_dof_indices(n_dofs_per_cell); + + const double ghost_parameter = 0.5; + const double nitsche_parameter = 5 * (fe_degree + 1) * fe_degree; + + // Since the ghost penalty is similar to a DG flux term, the simplest way to + // assemble it is to use an FEInterfaceValues object. + const QGauss face_quadrature(fe_degree + 1); + FEInterfaceValues fe_interface_values(fe_collection[0], + face_quadrature, + update_gradients | + update_JxW_values | + update_normal_vectors); + + + // As we iterate over the cells in the mesh, we would in principle have to + // do the following on each cell, $T$, + // + // 1. Construct one quadrature rule to integrate over the intersection with + // the domain, $T \cap \Omega$, and one quadrature rule to integrate over + // the intersection with the boundary, $T \cap \Gamma$. + // 2. Create FEValues-like objects with the new quadratures. + // 3. Assemble the local matrix using the created FEValues-objects. + // + // To make the assembly easier, we use the class NonMatching::FEValues, + // which does the above steps 1 and 2 for us. The algorithm @cite saye_2015 + // that is used to generate the quadrature rules on the intersected cells + // uses a 1-dimensional quadrature rule as base. Thus, we pass a 1D + // Gauss--Legendre quadrature to the constructor of NonMatching::FEValues. + // On the non-intersected cells, a tensor product of this 1D-quadrature will + // be used. + // + // As stated in the introduction, each cell has 3 different regions: inside, + // surface, and outside, where the level set function in each region is + // negative, zero, and positive. We need an UpdateFlags variable for each + // such region. These are stored on an object of type + // NonMatching::RegionUpdateFlags, which we pass to NonMatching::FEValues. + const QGauss<1> quadrature_1D(fe_degree + 1); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + region_update_flags.surface = update_values | update_gradients | + update_JxW_values | update_quadrature_points | + update_normal_vectors; + + NonMatching::FEValues non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set); + + // As we iterate over the cells, we don't need to do anything on the cells + // that have FENothing elements. To disregard them we use an iterator + // filter. + for (const auto &cell : + dof_handler.active_cell_iterators() | + IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange)) + { + local_stiffness = 0; + local_rhs = 0; + + const double cell_side_length = cell->minimum_vertex_distance(); + + // First, we call the reinit function of our NonMatching::FEValues + // object. In the background, NonMatching::FEValues uses the + // MeshClassifier passed to its constructor to check if the incoming + // cell is intersected. If that is the case, NonMatching::FEValues calls + // the NonMatching::QuadratureGenerator in the background to create the + // immersed quadrature rules. + non_matching_fe_values.reinit(cell); + + // After calling reinit, we can retrieve a dealii::FEValues object with + // quadrature points that corresponds to integrating over the inside + // region of the cell. This is the object we use to do the local + // assembly. This is similar to how hp::FEValues builds dealii::FEValues + // objects. However, one difference here is that the dealii::FEValues + // object is returned as an optional. This is a type that wraps an + // object that may or may not be present. This requires us to add an + // if-statement to check if the returned optional contains a value, + // before we use it. This might seem odd at first. Why does the function + // not just return a reference to a const FEValues? The reason is + // that in an immersed method, we have essentially no control of how the + // cuts occur. Even if the cell is formally intersected: $T \cap \Omega + // \neq \emptyset$, it might be that the cut is only of floating point + // size $|T \cap \Omega| \sim \epsilon$. When this is the case, we can + // not expect that the algorithm that generates the quadrature rule + // produces anything useful. It can happen that the algorithm produces 0 + // quadrature points. When this happens, the returned optional will not + // contain a value, even if the cell is formally intersected. + const std_cxx17::optional> &inside_fe_values = + non_matching_fe_values.get_inside_fe_values(); + + if (inside_fe_values) + for (const unsigned int q : + inside_fe_values->quadrature_point_indices()) + { + const Point &point = inside_fe_values->quadrature_point(q); + for (const unsigned int i : inside_fe_values->dof_indices()) + { + for (const unsigned int j : inside_fe_values->dof_indices()) + { + local_stiffness(i, j) += + inside_fe_values->shape_grad(i, q) * + inside_fe_values->shape_grad(j, q) * + inside_fe_values->JxW(q); + } + local_rhs(i) += rhs_function.value(point) * + inside_fe_values->shape_value(i, q) * + inside_fe_values->JxW(q); + } + } + + // In the same way, we can use NonMatching::FEValues to retrieve an + // FEFaceValues-like object to integrate over $T \cap \Gamma$. The only + // thing that is new here is the type of the object. The transformation + // from quadrature weights to JxW-values is different for surfaces, so + // we need a new class: NonMatching::FEImmersedSurfaceValues. In + // addition to the ordinary functions shape_value(..), shape_grad(..), + // etc., one can use its normal_vector(..)-function to get an outward + // normal to the immersed surface, $\Gamma$. In terms of the level set + // function, this normal reads + // @f{equation*} + // n = \frac{\nabla \psi}{\| \nabla \psi \|}. + // @f} + // An additional benefit of std::optional is that we do not need any + // other check for whether we are on intersected cells: In case we are + // on an inside cell, we get an empty object here. + const std_cxx17::optional> + &surface_fe_values = non_matching_fe_values.get_surface_fe_values(); + + if (surface_fe_values) + { + for (const unsigned int q : + surface_fe_values->quadrature_point_indices()) + { + const Point &point = + surface_fe_values->quadrature_point(q); + const Tensor<1, dim> &normal = + surface_fe_values->normal_vector(q); + for (const unsigned int i : surface_fe_values->dof_indices()) + { + for (const unsigned int j : + surface_fe_values->dof_indices()) + { + local_stiffness(i, j) += + (-normal * surface_fe_values->shape_grad(i, q) * + surface_fe_values->shape_value(j, q) + + -normal * surface_fe_values->shape_grad(j, q) * + surface_fe_values->shape_value(i, q) + + nitsche_parameter / cell_side_length * + surface_fe_values->shape_value(i, q) * + surface_fe_values->shape_value(j, q)) * + surface_fe_values->JxW(q); + } + local_rhs(i) += + boundary_condition.value(point) * + (nitsche_parameter / cell_side_length * + surface_fe_values->shape_value(i, q) - + normal * surface_fe_values->shape_grad(i, q)) * + surface_fe_values->JxW(q); + } + } + } + + cell->get_dof_indices(local_dof_indices); + + stiffness_matrix.add(local_dof_indices, local_stiffness); + rhs.add(local_dof_indices, local_rhs); + + // The assembly of the ghost penalty term is straight forward. As we + // iterate over the local faces, we first check if the current face + // belongs to the set $\mathcal{F}_h$. The actual assembly is simple + // using FEInterfaceValues. Assembling in this we will traverse each + // internal face in the mesh twice, so in order to get the penalty + // constant we expect, we multiply the penalty term with a factor 1/2. + for (unsigned int f : cell->face_indices()) + if (face_has_ghost_penalty(cell, f)) + { + const unsigned int invalid_subface = + numbers::invalid_unsigned_int; + + fe_interface_values.reinit(cell, + f, + invalid_subface, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + invalid_subface); + + const unsigned int n_interface_dofs = + fe_interface_values.n_current_interface_dofs(); + FullMatrix local_stabilization(n_interface_dofs, + n_interface_dofs); + for (unsigned int q = 0; + q < fe_interface_values.n_quadrature_points; + ++q) + { + const Tensor<1, dim> normal = fe_interface_values.normal(q); + for (unsigned int i = 0; i < n_interface_dofs; ++i) + for (unsigned int j = 0; j < n_interface_dofs; ++j) + { + local_stabilization(i, j) += + .5 * ghost_parameter * cell_side_length * normal * + fe_interface_values.jump_in_shape_gradients(i, q) * + normal * + fe_interface_values.jump_in_shape_gradients(j, q) * + fe_interface_values.JxW(q); + } + } + + const std::vector + local_interface_dof_indices = + fe_interface_values.get_interface_dof_indices(); + + stiffness_matrix.add(local_interface_dof_indices, + local_stabilization); + } + } + } + + + // @sect3{Solving the System} + template + void LaplaceSolver::solve() + { + std::cout << "Solving system" << std::endl; + + const unsigned int max_iterations = solution.size(); + SolverControl solver_control(max_iterations); + SolverCG<> solver(solver_control); + solver.solve(stiffness_matrix, solution, rhs, PreconditionIdentity()); + } + + + + // @sect3{Data Output} + // Since both DoFHandler instances use the same triangulation, we can add both + // the level set function and the solution to the same vtu-file. Further, we + // do not want to output the cells that have LocationToLevelSet value outside. + // To disregard them, we write a small lambda function and use the + // set_cell_selection function of the DataOut class. + template + void LaplaceSolver::output_results() const + { + std::cout << "Writing vtu file" << std::endl; + + DataOut data_out; + data_out.add_data_vector(dof_handler, solution, "solution"); + data_out.add_data_vector(level_set_dof_handler, level_set, "level_set"); + + data_out.set_cell_selection( + [this](const typename Triangulation::cell_iterator &cell) { + return cell->is_active() && + mesh_classifier.location_to_level_set(cell) != + NonMatching::LocationToLevelSet::outside; + }); + + data_out.build_patches(); + std::ofstream output("step-85.vtu"); + data_out.write_vtu(output); + } + + + + // @sect3{$L^2$-Error} + // To test that the implementation works as expected, we want to compute the + // error in the solution in the $L^2$-norm. The analytical solution to the + // Poisson problem stated in the introduction reads + // @f{align*} + // u(x) = 1 - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in + // \overline{\Omega}. + // @f} + // We first create a function corresponding to the analytical solution: + template + class AnalyticalSolution : public Function + { + public: + double value(const Point & point, + const unsigned int component = 0) const override; + }; + + + + template + double AnalyticalSolution::value(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + return 1. - 2. / dim * (point.norm_square() - 1.); + } + + + + // Of course, the analytical solution, and thus also the error, is only + // defined in $\overline{\Omega}$. Thus, to compute the $L^2$-error we must + // proceed in the same way as when we assembled the linear system. We first + // create an NonMatching::FEValues object. + template + double LaplaceSolver::compute_L2_error() const + { + std::cout << "Computing L2 error" << std::endl; + + const QGauss<1> quadrature_1D(fe_degree + 1); + + NonMatching::RegionUpdateFlags region_update_flags; + region_update_flags.inside = + update_values | update_JxW_values | update_quadrature_points; + + NonMatching::FEValues non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set); + + // We then iterate iterate over the cells that have LocationToLevelSetValue + // value inside or intersected again. For each quadrature point, we compute + // the pointwise error and use this to compute the integral. + const AnalyticalSolution analytical_solution; + double error_L2_squared = 0; + + for (const auto &cell : + dof_handler.active_cell_iterators() | + IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange)) + { + non_matching_fe_values.reinit(cell); + + const std_cxx17::optional> &fe_values = + non_matching_fe_values.get_inside_fe_values(); + + if (fe_values) + { + std::vector solution_values(fe_values->n_quadrature_points); + fe_values->get_function_values(solution, solution_values); + + for (const unsigned int q : fe_values->quadrature_point_indices()) + { + const Point &point = fe_values->quadrature_point(q); + const double error_at_point = + solution_values.at(q) - analytical_solution.value(point); + error_L2_squared += + std::pow(error_at_point, 2) * fe_values->JxW(q); + } + } + } + + return std::sqrt(error_L2_squared); + } + + + + // @sect3{A Convergence Study} + // Finally, we do a convergence study to check that the $L^2$-error decreases + // with the expected rate. We refine the background mesh a few times. In each + // refinement cycle, we solve the problem, compute the error, and add the + // $L^2$-error and the mesh size to a ConvergenceTable. + template + void LaplaceSolver::run() + { + ConvergenceTable convergence_table; + const unsigned int n_refinements = 3; + + make_grid(); + for (unsigned int cycle = 0; cycle <= n_refinements; cycle++) + { + std::cout << "Refinement cycle " << cycle << std::endl; + triangulation.refine_global(1); + setup_discrete_level_set(); + std::cout << "Classifying cells" << std::endl; + mesh_classifier.reclassify(); + distribute_dofs(); + initialize_matrices(); + assemble_system(); + solve(); + if (cycle == 1) + output_results(); + const double error_L2 = compute_L2_error(); + const double cell_side_length = + triangulation.begin_active()->minimum_vertex_distance(); + + convergence_table.add_value("Cycle", cycle); + convergence_table.add_value("Mesh size", cell_side_length); + convergence_table.add_value("L2-Error", error_L2); + + convergence_table.evaluate_convergence_rates( + "L2-Error", ConvergenceTable::reduction_rate_log2); + convergence_table.set_scientific("L2-Error", true); + + std::cout << std::endl; + convergence_table.write_text(std::cout); + std::cout << std::endl; + } + } + +} // namespace Step85 + + + +// @sect3{The main() function} +int main() +{ + const int dim = 2; + + Step85::LaplaceSolver laplace_solver; + laplace_solver.run(); +}