From: Vladimir Yushutin Date: Wed, 4 Oct 2023 21:20:25 +0000 (-0400) Subject: The tutorial step-90 demonstrates an efficient implementation of the stabilized Trace... X-Git-Tag: v9.6.0-rc1~196^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16878%2Fhead;p=dealii.git The tutorial step-90 demonstrates an efficient implementation of the stabilized TraceFEM using adaptive refinement and meshworkers. (Squashed after the review process) --- diff --git a/doc/doxygen/images/step-90-solution.png b/doc/doxygen/images/step-90-solution.png new file mode 100644 index 0000000000..0f90965953 Binary files /dev/null and b/doc/doxygen/images/step-90-solution.png differ diff --git a/doc/doxygen/images/step-90_mesh_cut.png b/doc/doxygen/images/step-90_mesh_cut.png new file mode 100644 index 0000000000..5367d05a96 Binary files /dev/null and b/doc/doxygen/images/step-90_mesh_cut.png differ diff --git a/doc/doxygen/images/step-90_prelim.png b/doc/doxygen/images/step-90_prelim.png new file mode 100644 index 0000000000..0b74116f65 Binary files /dev/null and b/doc/doxygen/images/step-90_prelim.png differ diff --git a/doc/doxygen/images/step-90_surface.png b/doc/doxygen/images/step-90_surface.png new file mode 100644 index 0000000000..0b9fb8a82b Binary files /dev/null and b/doc/doxygen/images/step-90_surface.png differ diff --git a/doc/doxygen/images/step-90_weak-vs-strong.png b/doc/doxygen/images/step-90_weak-vs-strong.png new file mode 100644 index 0000000000..0de7c87982 Binary files /dev/null and b/doc/doxygen/images/step-90_weak-vs-strong.png differ diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 85c081c077..cc963dd0b2 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1710,6 +1710,20 @@ pages = {472--501} } +%------------------------------------------------------------------------------- +% Step 90 +%------------------------------------------------------------------------------- +@InProceedings{traceFEM_review_2017, + author="Olshanskii, Maxim A. + and Reusken, Arnold", + editor="Bordas, St{\'e}phane P. A. and Burman, Erik and Larson, Mats G. and Olshanskii, Maxim A.", + title="Trace Finite Element Methods for PDEs on Surfaces", + booktitle="Geometrically Unfitted Finite Element Methods and Applications", + year="2017", + publisher="Springer International Publishing", + pages="211--258", + isbn="978-3-319-71431-8" +} %------------------------------------------------------------------------------- % Step 87 diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 89351c0cd5..d4b329efa9 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -697,6 +697,11 @@ *
Keywords: FERemoteEvaluation * * + * + * step-90 + * Solving the Laplace-Beltrami equation on a surface using the trace finite element method. + *
Keywords: MeshWorker::mesh_loop(), NonMatching::FEImmersedSurfaceValues + * * * * diff --git a/examples/step-90/CMakeLists.txt b/examples/step-90/CMakeLists.txt new file mode 100644 index 0000000000..8ead29f18f --- /dev/null +++ b/examples/step-90/CMakeLists.txt @@ -0,0 +1,59 @@ +## +# CMake script for the step-90 tutorial program: +## + +# Set the name of the project and target: +set(TARGET "step-90") + +# 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 +) + +# Define the output that should be cleaned: +set(CLEAN_UP_FILES *.vtu *.pvtu *.visit) + +# Usually, you will not need to modify anything beyond this point... + +cmake_minimum_required(VERSION 3.13.4) + +find_package(deal.II 9.6.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() + +# +# Are all dependencies fulfilled? +# +if(NOT DEAL_II_WITH_MPI OR NOT DEAL_II_WITH_P4EST OR NOT DEAL_II_WITH_TRILINOS) # keep in one line + message(FATAL_ERROR " +Error! This tutorial requires a deal.II library that was configured with the following options: + DEAL_II_WITH_MPI = ON + DEAL_II_WITH_P4EST = ON + DEAL_II_WITH_TRILINOS = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options: + DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI} + DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} + DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS} +This conflicts with the requirements." + ) +endif() + +deal_ii_initialize_cached_variables() +project(${TARGET}) +deal_ii_invoke_autopilot() diff --git a/examples/step-90/doc/builds-on b/examples/step-90/doc/builds-on new file mode 100644 index 0000000000..b1e52c0ace --- /dev/null +++ b/examples/step-90/doc/builds-on @@ -0,0 +1 @@ +step-85 diff --git a/examples/step-90/doc/intro.dox b/examples/step-90/doc/intro.dox new file mode 100644 index 0000000000..4318c4e16e --- /dev/null +++ b/examples/step-90/doc/intro.dox @@ -0,0 +1,141 @@ + +This program was contributed by Vladimir Yushutin and Timo Heister, Clemson University, 2023. + +This material is based upon work partly supported by the National +Science Foundation Award DMS- + + + +

Introduction

+ +In this tutorial, we implement the trace finite element method (TraceFEM) in deal.II. TraceFEM solves PDEs posed on a +possibly evolving $(dim-1)$-dimensional surface $\Gamma$ employing a fixed uniform background mesh of a $dim$-dimensional domain +in which the surface is embedded. Such surface PDEs arise in problems involving material films with complex +properties and in other situations in which a non-trivial condition is imposed on either a stationary or a moving interface. +Here we consider a steady, complex, non-trivial surface and the prototypical Laplace-Beltrami equation which is a counterpart of +the Poisson problem on flat domains. + +Being an unfitted method, TraceFEM allows to circumvent the need of remeshing of an evolving surface if it is implicitly +given by the zero contour of a level-set function. At the same time, it easily provides with an extension of the +discrete solution to a neighborhood of the surface which turns out to be very handy in case of non-stationary interfaces and films. +Certainly, this flexibility comes with a price: one needs to design the nodes and weights for a quadrature customized +for each implicit intersection of the zero level-set and the background mesh. Moreover, these intersections may be of +arbitrary shape and size manifesting in the so-called ``small cut" problem and requiring addition of a stabilization +form that restores well-conditioning of the problem. + +Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE, +e.g. a $Q_2$ discrete level-set and a $Q_1$ solution are possible on the same bulk triangulation. +Second, we make sure that the performance of TraceFEM in the parallel implementation corresponds to that of a classical +fitted FEM for a two-dimensional problem. We demonstrate how to achieve both goals by using a combination of MeshWorker +and NonMatching capabilities. + +A natural alternative to TraceFEM in solving surface PDEs is the parametric surface finite element method. The latter +method relies on an explicit parametrization of the surface which may be not feasible especially for evolving interfaces +with an unknown in advance shape - in this sense, TraceFEM is a technique inspired by the level-set description of +interfaces. However, the parametric surface finite element method, when applicable, enjoys many well-known properties +of fitted methods on flat domains provided the geometric errors - which a present for both methods - are taken under control. + + +

A non-trivial surface

+A fitted FEM on a flat two-dimensional domain, if discretized by piecewise linears with $N$ degrees of freedom, typically results in +$O(h)=O(N^{-1/2})$ convergence rate of the energy error; requires $O(N)$ storage for the degrees of freedom; and, +more importantly, takes $O(N)$ of construction time to create them, i.e. to mesh the domain. TraceFEM, +although solving a two-dimensional problem, relies on the inherently three-dimensional mesh on which the level-set +function must be defined and, if implemented naively, suffers from the increased storage and the increased construction +time in terms of the active degrees of freedom $N_a$ that actually +enters the scheme with, hopefully, $O(N_a^{-1/2})$ error. To combat these possible bottlenecks, we create iteratively +a mesh which is localized near the zero contour line of the level set function, i.e near the surface, to restore the aforementioned +two-dimensional performance typical for fitted FEM, see the first three typical iterations of this methodology below. + +@image html step-90_prelim.png "Iterative localization of the zero contour of a typical level set" width=60% + +The cells colored by red cary the active degrees of freedom (total number $N_a$) as the level set is not sign-definite +at support points. Notice also that the mesh is graded: any cell has at most 4 neighbors adjacent to each of 6 faces. + +Once a desired geometry approximation $\Gamma_h$ is achieved using the iterative approach above, we can start forming the linear system +using the constructed normals and quadratures. For the purposes of the tutorial we choose a non-trivial surface $\Gamma$ given by +@f{equation*} + \frac{x^2}{4}+ y^2 + \frac{4 z^2} {(1 + 0.5 \sin(\pi x))^{2}} = 1 +@f} +The OY and OX views of this tamarind-shaped, exact surface $\Gamma$ are shown below along with the mesh after +three iterations (the approximation $\Gamma_h$ is not shown). + +@image html step-90_surface.png "OY(left) and OZ(right) cross-sections of the background mesh along with the exact surface" width=80% + + +

Model problem

+ +We would like to solve the simplest possible problem defined on a surface, namely the Laplace--Beltrami equation, +@f{equation*} + -\Delta_\Gamma u + c u = f \qquad \text{in }\, \Gamma, +@f} +where we take $c=1$ for concreteness. We added the term $cu$ to the left-hand side so the problem becomes well-posed +in the absence of any boundary; an alternative could be to take $c=0$ but impose the zero mean condition. + +

Manufactured exact solution

+We choose the test solution and the right-hand side forcing + as the restriction to $\Gamma$ of +@f{equation*} + u(x,y,z)=xy\,,\quad + f(x,y,z)=xy + 2.0\,\mathbf{n}_x \mathbf{n}_y + \kappa (y \mathbf{n}_x + x\mathbf{n}_y), +@f} +where the latter is manufactured using the exact normal $\mathbf{n}$, the exact Hessian $\nabla^2\mathbf{n}$ and the mean curvature, +$\kappa=\mathrm{div} n$ of the surface. Note that we do not need to impose any boundary conditions as the surface $\Gamma$ is closed. + +

The Trace Finite Element Method

+TraceFEM is an unfitted method: the surface $\Gamma$ is immersed into a regular, uniform background mesh that +stays fixed even if the surface would be evolving. +To solve Laplace--Beltrami equation, we first construct a surface approximation $\Gamma_h$ by intersecting implicitly +the cells of the background mesh with the iso surface of an approximation of the level-set field. We note that we +never actually create any two-dimensional meshes for the surface but only compute approximate quadrature points and surface normals. +Next we distribute degrees of freedom over a thin subdomain $\Omega_h$ +that completely covers $\Gamma_h$ and that consists of the intersected cells $\mathcal{T}_\Gamma^h$, +@f{equation*} + \mathcal{T}_\Gamma^h = \{ T \in \mathcal{T}^{h} : T \cap \Gamma_h \neq \emptyset \}. +@f} +The finite element space where we want to find our numerical solution, $u_h$, is now +@f{equation*} + V_h = \{ v \in C(\Omega_h) : v \in Q_p(T), \, T \in \mathcal{T}_\Gamma^h \}, +@f} +where $\Omega_h$ is the union of all intersected cells from $\bigcup_{T \in \mathcal{T}_\Gamma^h} \overline{T}$. + +To create $V_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}_\Gamma^h$ or not, +we set the active_fe_index to either 0 or 1. +To determine whether a cell is intersected or not, we use the class NonMatching::MeshClassifier. + +A natural candidate for a weak formulation involves the following (bi)linear forms +@f{align*} + a_h(u_h, v_h) = (\nabla_{\Gamma_h} u_h, \nabla_{\Gamma_h} v_h)_{\Gamma_h}+(u_h, v_h)_{\Gamma_h}\,,\qquad + L_h(v_h) = (f^e,v_h)_{\Gamma_h}. +@f} +where $f^e$ is an extension (non-necessarily the the so-called normal extension) of $f$ from $\Gamma$ to $\Omega_h$. Note that the right-hand side $f$ of the Laplace-Beltrami +problem is defined on the exact surface $\Gamma$ only and we need to specify how to evaluate its action on the perturbed +approximate geometry $\Gamma_h$ which is immersed in $\Omega_h$. For the purposes of this test, the forcing $f$ is +manufactured using $u=xy$ and the level-set function and, therefore, is a function of Cartesian coordinates $x$,$y$,$z$. +The latter is identified with $f^e$ on $\Gamma_h$ and it is not the normal extension of the function $f$. + +However, the so-called "small-cut problem" may arise and one should +introduce the stabilized version of TraceFEM: Find $u_h \in V_h$ such that +@f{equation*} + a_h(u_h,v_h) + s_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h. +@f} +Here the normal-gradient stabilization $s_h$ involves the three-dimensional integration over whole (but intersected) cells and is given by +@f{equation*} + s_h(u_h,v_h) = h^{-1}(\mathbf{n}_h\cdot\nabla u_h, \mathbf{n}_h\cdot\nabla v_h)_{\Omega_h}, +@f} +Note that the $h^{-1}$ scaling may be relaxed for sufficiently smooth solutions such as the manufactured one, but we +choose the strong scaling to demonstrate the extreme case @cite traceFEM_review_2017. + +

Discrete Level Set Function

+In TraceFEM we construct the approximation $\Gamma_h$ using the interpolant $\psi_h$ of the exact level-set function on the bulk triangulation: +@f{align*} + \Gamma_h &= \{x \in \mathbb{R}^{\text{3}} : \psi_h(x) = 0 \}. +@f} +The exact normal vector $\mathbf{n}$ is approximated by $\mathbf{n}_h=\nabla\psi_h/\|\nabla\psi_h\|$ which, together +with approximate quadrature for the integration over $\Gamma_h$, leads to the so-called "geometrical error". +Luckily, one can show @cite traceFEM_review_2017 that the method converges optimally for the model problem +if the same element space $V_h$ is employed for the discrete functions and for the interpolation of the level set +function as if the exact domain would have been used. Furthermore, deal.II allows to choose independently the discrete +space for the solution and a higher-order discrete space for the level set function for a more accurate geometric approximation. diff --git a/examples/step-90/doc/kind b/examples/step-90/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-90/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-90/doc/results.dox b/examples/step-90/doc/results.dox new file mode 100644 index 0000000000..bbb954625c --- /dev/null +++ b/examples/step-90/doc/results.dox @@ -0,0 +1,37 @@ +

Results

+ +The numerical solution $u_h$ for a very fine mesh $\Gamma_h$ is shown below by plotting in Paraview the zero contour of +the approximate level set $\psi_h$ and restricting the discrete solution $u_h$ to the resulting surface approximation $\Gamma_h$. + +@image html step-90-solution.png width=50% + +Next, we demonstrate the corresponding set of intersected cells with active degrees of freedom. Note that not all cells +are of the same refinement level which is attributed to the insufficiently fine initial uniform grid. + +@image html step-90_mesh_cut.png width=50% + +

Convergence test

+ +The results of the convergence study are shown in the following table. The experimental orders of convergence (EOC) +are reported for the surface errors and the stabilization. + +| Cycle | DOFS | Rate | Iterations | $L^2$-Error | EOC | $H^1$-Error | EOC |$s_h^{1/2}(u_h)$| EOC | +|:-----:|:--------:|:----:|:----------:|:-----------:|:-----:|:-----------:|:----:|:--------------:|:-----:| +| 0 | 12370 | - | 15 | 7.6322e-02 | - | 3.6212e-01 | - | 2.2423e-01 | - | +| 1 | 49406 | 2.00 | 18 | 1.1950e-02 | 2.68 | 1.4752e-01 | 1.30 | 1.1238e-01 | 1.00 | +| 2 | 196848 | 1.99 | 19 | 1.7306e-03 | 2.79 | 7.4723e-02 | 0.98 | 6.1131e-02 | 0.88 | +| 3 | 785351 | 2.00 | 22 | 3.6276e-04 | 2.25 | 3.9329e-02 | 0.93 | 3.0185e-02 | 1.02 | +| 4 | 3136501 | 2.00 | 25 | 7.5910e-05 | 2.26 | 1.9694e-02 | 1.00 | 1.4875e-02 | 1.02 | +| 5 | 12536006 | 2.00 | 26 | 1.7279e-05 | 2.14 | 9.8443e-03 | 1.00 | 7.4067e-03 | 1.01 | +| 6 | 50122218 | 2.00 | 30 | 4.3891e-06 | 1.98 | 4.9219e-03 | 1.00 | 3.7042e-03 | 1.00 | + +In this test we refine the mesh near the surface and, as a result, the number of degrees of freedom scales in the two-dimensional fashion. +The optimal rates of error convergence in $L^2(\Gamma)$ and $H^1(\Gamma)$ norms are clearly observable. We also note +the first order convergence of the stabilization $s_h^{1/2}(u_h)=\sqrt{s_h(u_h, u_h)}$ evaluated at the solution $u_h$. + +

Parallel scalability

+ +The weak and strong scalability test results are shown in the following figure. Clearly, the refine() method is +responsible for the certain lack of parallel scalability. + +@image html step-90_weak-vs-strong.png width=100% diff --git a/examples/step-90/doc/tooltip b/examples/step-90/doc/tooltip new file mode 100644 index 0000000000..3f5c006a06 --- /dev/null +++ b/examples/step-90/doc/tooltip @@ -0,0 +1 @@ +Solving the Laplace-Beltrami equation on a surface using the trace finite element method. diff --git a/examples/step-90/step-90.cc b/examples/step-90/step-90.cc new file mode 100644 index 0000000000..f97e9c91b3 --- /dev/null +++ b/examples/step-90/step-90.cc @@ -0,0 +1,1243 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2024 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 Vladimir Yushutin and Timo Heister, Clemson + * University, 2023. + */ + +#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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; +using VectorType = TrilinosWrappers::MPI::Vector; +using MatrixType = TrilinosWrappers::SparseMatrix; +namespace Step90 +{ + // The parallization in this tutorial relies on the Trilinos library. We will + // grant to some cells empty finite element spaces FE_Nothing as done + // in step-85, but this time active DoFs will be only assigned to cell which + // are intersected by the surface approximation. + enum class ActiveFEIndex : types::fe_index + { + lagrange = 0, + nothing = 1 + }; + + // @sect3{Exact surface} + // The following class defines the surface using the implicit level set + // representation. The exact surface normal uses the Cartesian gradient of the + // level set function. The exact Hessian is needed for the construction of the + // test case only. + template + class TamarindShape : public Function + { + public: + TamarindShape() + : Function() + {} + double value(const Point &point, + const unsigned int component = 0) const override + { + AssertIndexRange(component, this->n_components); + (void)component; + Assert(dim == 3, ExcNotImplemented()); + + return 0.25 * Utilities::pow(point[0], 2) + Utilities::pow(point[1], 2) + + 4.0 * Utilities::pow(point[2], 2) * + std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2) - + 1.0; + } + + Tensor<1, dim> gradient(const Point &point, + const unsigned int component = 0) const override + { + AssertIndexRange(component, this->n_components); + (void)component; + Assert(dim == 3, ExcNotImplemented()); + + Tensor<1, dim> grad; + grad[0] = 0.5 * point[0] + + (-2.0) * 4.0 * Utilities::pow(point[2], 2) * + std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) * + (0.5 * numbers::PI * std::cos(numbers::PI * point[0])); + grad[1] = 2.0 * point[1]; + grad[2] = (2.0) * 4.0 * point[2] * + std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2); + + return grad; + } + + SymmetricTensor<2, dim> + hessian(const Point &point, + const unsigned int component = 0) const override + { + AssertIndexRange(component, this->n_components); + (void)component; + Assert(dim == 3, ExcNotImplemented()); + + SymmetricTensor<2, dim> hessian; + + hessian[0][0] = + 0.5 + + 8.0 * Utilities::pow(point[2], 2) * + (3.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -4) * + Utilities::pow(0.5 * numbers::PI * + std::cos(numbers::PI * point[0]), + 2) + + std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) * 0.5 * + numbers::PI * numbers::PI * std::sin(numbers::PI * point[0])); + hessian[0][1] = 0.0; + hessian[0][2] = + (-8.0) * point[2] * + std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) * + numbers::PI * std::cos(numbers::PI * point[0]); + + hessian[1][1] = 2.0; + hessian[1][2] = 0.0; + + hessian[2][2] = + 8.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2); + + return hessian; + } + }; + + // @sect3{Exact solution} + // The following class defines the chosen exact solution and its surface + // gradient. The exact solution we try to reproduce is $u=xy$ and it may be + // evaluated away from + // $\Gamma$ as any other function of Cartesian points. Also note that the + // gradient() method returns the surface gradient $\nabla_\Gamma u$ of the + // exact solution. + template + class AnalyticalSolution : public Function + { + private: + const TamarindShape tamarind; + + public: + AnalyticalSolution() + : Function() + {} + double value(const Point &point, + const unsigned int component = 0) const override; + + Tensor<1, dim> gradient(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 point[0] * point[1]; + } + + template + Tensor<1, dim> + AnalyticalSolution::gradient(const Point &point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + const Tensor<1, dim> grad = tamarind.gradient(point, component); + const Tensor<1, dim> normal = grad / grad.norm(); + + Tensor<1, dim> projector_first_column = -normal[0] * normal; + projector_first_column[0] += 1.0; + + Tensor<1, dim> projector_second_column = -normal[1] * normal; + projector_second_column[1] += 1.0; + + Tensor<1, dim> surface_gradient = + point[1] * projector_first_column + point[0] * projector_second_column; + + return surface_gradient; + } + + // @sect3{Exact forcing} + // We choose the right hand side equal to the evaluation of the surface + // Laplacian for a manufactured solution $u$. + // This corresponds to the exact forcing $f=-\Delta_\Gamma u+u$: + template + class RightHandSide : public Function + { + const TamarindShape tamarind; + + public: + RightHandSide() + : Function() + {} + + virtual double value(const Point &p, + const unsigned int component = 0) const override; + }; + + template + double RightHandSide::value(const Point &point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + Assert(dim == 3, ExcNotImplemented()); + + const Tensor<1, dim> grad = tamarind.gradient(point, component); + const Tensor<1, dim> normal = grad / grad.norm(); + const SymmetricTensor<2, dim> hessian = tamarind.hessian(point, component); + + double mean_curv = 0.0; + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + mean_curv += ((j == k ? 1 : 0) - normal[j] * normal[k]) * hessian[j][k]; + mean_curv /= grad.norm(); + + return point[0] * point[1] + 2.0 * normal[0] * normal[1] + + mean_curv * (point[1] * normal[0] + point[0] * normal[1]); + } + + // @sect3{Scratch and Copy objects for TraceFEM} + // Since the assembly procedure will be performed via MeshWorker, we need a + // Scratch object that handles the Non-Matching FEValues effectively. + // The input arguments of its constructor are discussed in the solver class + // below. + template + struct ScratchData + { + ScratchData(const Mapping &mapping, + const hp::FECollection &fe_collection, + const NonMatching::MeshClassifier &mesh_classifier, + const DoFHandler &level_set_dof_handler, + const VectorType &level_set, + const NonMatching::RegionUpdateFlags nonmatching_update_flags, + const Quadrature &quadrature, + const Quadrature<1> &quadrature_edge, + const UpdateFlags cell_update_flags = update_values | + update_gradients | + update_quadrature_points | + update_JxW_values) + : fe_values( + mapping, + fe_collection[static_cast(ActiveFEIndex::lagrange)], + quadrature, + cell_update_flags) + , region_update_flags(nonmatching_update_flags) + , quadrature_1D(quadrature_edge) + , fe_collection(fe_collection) + , mesh_classifier(mesh_classifier) + , level_set_dof_handler(level_set_dof_handler) + , level_set(level_set) + , level_set_fe_values(mapping, + level_set_dof_handler.get_fe(), + quadrature, + cell_update_flags) + , non_matching_fe_values(fe_collection, + quadrature_edge, + nonmatching_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set) + {} + + ScratchData(const ScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_mapping(), + scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + scratch_data.fe_values.get_update_flags()) + , region_update_flags(scratch_data.region_update_flags) + , quadrature_1D(scratch_data.quadrature_1D) + , fe_collection(scratch_data.fe_collection) + , mesh_classifier(scratch_data.mesh_classifier) + , level_set_dof_handler(scratch_data.level_set_dof_handler) + , level_set(scratch_data.level_set) + , level_set_fe_values(scratch_data.level_set_fe_values.get_mapping(), + scratch_data.level_set_fe_values.get_fe(), + scratch_data.level_set_fe_values.get_quadrature(), + scratch_data.level_set_fe_values.get_update_flags()) + , non_matching_fe_values(fe_collection, + quadrature_1D, + region_update_flags, + mesh_classifier, + level_set_dof_handler, + level_set) + {} + + // The following FEValues object is used for the standard quadrature on + // cells involving the FE space of the solution. In TraceFEM, we need this + // quadrature due to the stabilization term. In addition, a cell quadrature + // for the FE space of the level set is defined. + FEValues fe_values; + const NonMatching::RegionUpdateFlags region_update_flags; + const Quadrature<1> &quadrature_1D; + const hp::FECollection &fe_collection; + const NonMatching::MeshClassifier &mesh_classifier; + const DoFHandler &level_set_dof_handler; + const VectorType &level_set; + FEValues level_set_fe_values; + NonMatching::FEValues non_matching_fe_values; + }; + + // The MeshWorker framework also requires a "copy" data structure that is + // filled by the worker function working on a cell or face, and whose contents + // are then later copied into global matrices and vectors. This CopyData + // object is customized for TraceFEM. In particular, the implementation of the + // normal-gradient volume stabilization relies on it. + template + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + + void reinit(const typename DoFHandler::active_cell_iterator &cell) + { + 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->get_dof_indices(local_dof_indices); + } + }; + + template + struct CopyDataError + { + unsigned int cell_index; + double cell_L2_error_sqr; + double cell_H1_error_sqr; + double cell_stab_sqr; + + void reinit(const typename DoFHandler::active_cell_iterator &cell) + { + cell_index = cell->active_cell_index(); + cell_L2_error_sqr = 0.0; + cell_H1_error_sqr = 0.0; + cell_stab_sqr = 0.0; + } + }; + + // @sect3{Normal-gradient stabilization form of TraceFEM} + // The following class corresponds to the stabilization form, + // its contribution to the global matrix and to the error. More specifically, + // the method needs_cell_worker() indicates + // whether the bilinear form of the stabilization, unlike the main bilinear + // form of Laplace-Beltrami operator, needs the bulk cell quadratures. The + // cell worker which is useful in an accumulation by MeshWorkers is provided + // by the assemble_cell_worker() method. The remaining method + // evaluate_cell_worker() computes the stabilization error for the solution + // $u_h$, i.e $s_h(u_h,u_h)$. Also note that the method needs_cell_worker() + // indicates that the assembly and the evaluation of the form does require a + // bulk cell quadrature. This methodology may be utilized in the MeshWorker. + // The stabilization scaling is specified by + // $\mathrm{stabilization\_parameter}\cdot + // h^\mathrm{stabilization\_exponent}$. For elliptic problems with smooth + // solutions we can choose any + // $-1\leq \mathrm{stabilization\_exponent} \leq 1$ and + // a sufficiently large $\mathrm{stabilization\_parameter}$ that depends of + // $\Gamma$. + template + class NormalGradientVolumeStabilization + { + public: + NormalGradientVolumeStabilization() + : stabilization_parameter(1.0) + , stabilization_exponent(-1.0) + {} + + bool needs_cell_worker() const + { + return true; + } + + // We define the stabilization form here assuming that ScratchData and + // CopyData arguments are initialized properly. The local contribution of + // the stabilization from this cell to the global matrix is given in + // assemble_cell_worker() and, later in evaluate_cell_worker(), the + // local bilinear form of the stabilization is evaluated on the solution. + // Note the gradients of the discrete level set are computed + // in the bulk cell quadrature points, which, upon normalization, give the + // discrete normal vector in a bulk cell. + void assemble_cell_worker( + VectorType &level_set, + const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyData ©_data) const + { + const FEValues &fe_values = scratch_data.fe_values; + const FEValues &level_set_fe_values = + scratch_data.level_set_fe_values; + + const std::vector &JxW_cell = fe_values.get_JxW_values(); + + std::vector> grad_level_set( + level_set_fe_values.get_quadrature().size()); + level_set_fe_values.get_function_gradients(level_set, grad_level_set); + + const double factor = + stabilization_parameter * + std::pow(cell->minimum_vertex_distance(), stabilization_exponent); + for (const unsigned int q : fe_values.quadrature_point_indices()) + { + const Tensor<1, dim> &normal = + grad_level_set[q] / grad_level_set[q].norm(); + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + copy_data.cell_matrix(i, j) += + factor * (normal * fe_values.shape_grad(i, q)) * + (normal * fe_values.shape_grad(j, q)) * JxW_cell[q]; + } + } + + void evaluate_cell_worker( + VectorType &solution, + VectorType &level_set, + const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyDataError ©_data) const + { + double cell_stab_sqr = 0.0; + const FEValues &fe_values = scratch_data.fe_values; + const std::vector &JxW_cell = fe_values.get_JxW_values(); + const unsigned int n_q_points = fe_values.get_quadrature_points().size(); + const FEValues &level_set_fe_values = + scratch_data.level_set_fe_values; + + std::vector> level_set_grad(n_q_points); + level_set_fe_values.get_function_gradients(level_set, level_set_grad); + + std::vector> sol_grad(n_q_points); + fe_values.get_function_gradients(solution, sol_grad); + + const double factor = + stabilization_parameter * + std::pow(cell->minimum_vertex_distance(), stabilization_exponent); + + for (const unsigned int q : fe_values.quadrature_point_indices()) + { + const Tensor<1, dim> normal = + level_set_grad[q] / level_set_grad[q].norm(); + + const double stabilization_at_point = normal * sol_grad[q]; + cell_stab_sqr += + factor * Utilities::pow(stabilization_at_point, 2) * JxW_cell[q]; + } + copy_data.cell_stab_sqr = cell_stab_sqr; + } + + private: + const double stabilization_parameter; + const double stabilization_exponent; + }; + + // @sect3{Laplace--Beltrami solver} + // The main class whose method run() performs the computation. + // One may adjust main parameters of TraceFEM in the constructor. + // The other methods are discussed below. + template + class LaplaceBeltramiSolver + { + public: + LaplaceBeltramiSolver(); + void run(); + + private: + void make_grid(); + + void localize_surface(); + + void setup_discrete_level_set(); + + void distribute_dofs(); + + void initialize_matrices(); + + void assemble_system(); + + void solve(); + + void mark_intersected(); + + void refine_grid(); + + void compute_errors(); + + void output_level_set(unsigned int); + + void output_solution(); + + MPI_Comm mpi_communicator; + + // The surface of interest corresponds to the zero contour of the following + // exact level set function: + const TamarindShape tamarind; + + // The manufactured solution to the Laplace--Beltrami problem and the + // corresponding right-hand side: + const AnalyticalSolution analytical_solution; + const RightHandSide right_hand_side; + + // There is a single triangulation which is shared by the discretizations of + // the solution and of the level set. + parallel::distributed::Triangulation triangulation; + ConditionalOStream pcout; + TimerOutput computing_timer; + + // We need two separate FE spaces. + // The first manages the TraceFEM space which is active on intersected + // elements. The second manages the discrete + // level set function that describes the geometry of the surface. + // Also, the degrees of the FE spaces and the corresponding DoFHandler + // objects are given in the following: + const unsigned int fe_degree; + hp::FECollection fe_collection; + DoFHandler dof_handler; + + const unsigned int level_set_fe_degree; + const FE_Q level_set_fe; + DoFHandler level_set_dof_handler; + + const MappingQ1 mapping; + + // Since we will adaptively refine the bulk triangulation, two constraints + // are needed: one for the solution space and another for the level set + // space. + AffineConstraints constraints; + AffineConstraints level_set_constraints; + + // Discrete vectors initialized with dof_handler and level_set_dof_handler. + VectorType completely_distributed_solution; + VectorType locally_relevant_solution; + VectorType locally_relevant_exact; + VectorType level_set; + Vector active_fe_indicator; + + // The following NonMatching::MeshClassifier object is used to + // separate intersected elements and non-intersected ones. + // We will then use different finite elements from an hp::FECollection for + // these two categories: + NonMatching::MeshClassifier mesh_classifier; + + // The first bulk quadrature is required for the + // for TraceFEM stabilization, while the integration over implicit surface + // is based on the last, one-dimensional rule. + const QGauss cell_quadrature; + const QGauss<1> quadrature_1D; + + // Any TraceFEM needs a stabilization, and we choose the normal-gradient, + // volume stabilization. + const NormalGradientVolumeStabilization stabilization_scheme; + + // Discrete right-hand side and the final matrix corresponding to + // dof_handler. + VectorType global_rhs; + MatrixType global_matrix; + SparsityPattern sparsity_pattern; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + // Depending on the type of the quadrature, surface, face or volume, we need + // to define different update flags. + NonMatching::RegionUpdateFlags surface_update_flags; + + // The following variables are used to display the results of the + // convergence test: + ConvergenceTable convergence_table; + }; + + template + LaplaceBeltramiSolver::LaplaceBeltramiSolver() + : mpi_communicator(MPI_COMM_WORLD) + , tamarind() + , analytical_solution() + , right_hand_side() + , triangulation(mpi_communicator) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , computing_timer(mpi_communicator, + pcout, + TimerOutput::never, + TimerOutput::wall_times) + , fe_degree(1) + , fe_collection(FE_Q(fe_degree), FE_Nothing()) + , dof_handler(triangulation) + , level_set_fe_degree(1) + , level_set_fe(level_set_fe_degree) + , level_set_dof_handler(triangulation) + , mapping() + , mesh_classifier(level_set_dof_handler, level_set) + , cell_quadrature(fe_degree + 1) + , quadrature_1D(fe_degree + 1) + , stabilization_scheme() + { + surface_update_flags.surface = + update_values | update_gradients | update_JxW_values | + update_quadrature_points | update_normal_vectors; + } + + // @sect3{Geometric approximation} + // Let us start with a function that creates the background mesh, using a + // domain size chosen to avoid situations in which level set function vanishes + // at mesh vertices. The initial refinement helps the level set to approximate + // the surface meaningfully. + // + // In following next method we construct the discrete level set and determine + // which cells are intersected. Note that all cells, intersected and + // non-intersected, have a corresponding active_fe_indicator. + // Similarly, the exact level set function is approximated on the whole + // triangulation and postprocessed afterward resulting on a surface + // approximation with no gaps. + template + void LaplaceBeltramiSolver::make_grid() + { + pcout << "Creating background mesh..." + << "\n" + << std::flush; + const double cube_side = 2.008901281; + GridGenerator::hyper_cube(triangulation, -cube_side, cube_side); + triangulation.refine_global(3); + } + + template + void LaplaceBeltramiSolver::setup_discrete_level_set() + { + pcout + << "Setting up discrete level set function and reclassifying cells... " + << "\n" + << std::flush; + TimerOutput::Scope t(computing_timer, "setup_level_set"); + + active_fe_indicator.reinit(triangulation.n_active_cells()); + level_set_dof_handler.distribute_dofs(level_set_fe); + level_set_constraints.clear(); + const IndexSet level_set_locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(level_set_dof_handler); + level_set_constraints.reinit(level_set_locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(level_set_dof_handler, + level_set_constraints); + level_set_constraints.close(); + + // Here is where the geometric information enters the code. Next, using the + // discrete level set, we mark the cell which are intersected by its zero + // contour. Finally, once the triangulation's cells are classified, we + // determine which cells are active. + VectorType tmp_sol(level_set_dof_handler.locally_owned_dofs(), + mpi_communicator); + VectorTools::interpolate(level_set_dof_handler, tamarind, tmp_sol); + level_set_constraints.distribute(tmp_sol); + + level_set.reinit(level_set_locally_relevant_dofs, + level_set_dof_handler.locally_owned_dofs(), + mpi_communicator); + level_set = tmp_sol; + + mesh_classifier.reclassify(); + + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + if (mesh_classifier.location_to_level_set(cell) == + NonMatching::LocationToLevelSet::intersected) + cell->set_active_fe_index( + static_cast(ActiveFEIndex::lagrange)); + else + cell->set_active_fe_index( + static_cast(ActiveFEIndex::nothing)); + } + } + + // The method fills in the indicator telling which cells are intersected. It + // is used in the adaptive refinement near the surface. + template + void LaplaceBeltramiSolver::mark_intersected() + { + pcout << "Determining cells with active FE index..." + << "\n" + << std::flush; + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + if (mesh_classifier.location_to_level_set(cell) == + NonMatching::LocationToLevelSet::intersected) + active_fe_indicator[cell->active_cell_index()] = 1.0; + } + } + + + // We refine only intersected cells with active_fe_indicator=1. We are calling + // GridRefinement::refine_and_coarsen_fixed_fraction() instead of the + // GridRefinement::refine_and_coarsen_fixed_number() function called in most + // other tutorials because the number of non-intersected cells also grows + // interfering with the number of active, intersected cells. + template + void LaplaceBeltramiSolver::refine_grid() + { + TimerOutput::Scope t(computing_timer, "refine"); + pcout << "Refining near surface..." + << "\n" + << std::flush; + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction( + triangulation, active_fe_indicator, 1.0, 0.0); + + triangulation.execute_coarsening_and_refinement(); + } + + // As the surface is properly approximated by several adaptive steps, we may + // now distribute the degrees of + // freedom on cells which are intersected by the discrete approximation. + // Next, we initialize matrices for active DoFs and apply the constraints for + // the solution. + template + void LaplaceBeltramiSolver::distribute_dofs() + { + pcout << "Distributing degrees of freedom... " + << "\n" + << std::flush; + dof_handler.distribute_dofs(fe_collection); + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + completely_distributed_solution.reinit(dof_handler.locally_owned_dofs(), + mpi_communicator); + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + global_rhs.reinit(locally_owned_dofs, mpi_communicator); + + const unsigned int dof_handler_size = dof_handler.n_dofs(); + const unsigned int level_set_dof_handler_size = + level_set_dof_handler.n_dofs(); + + convergence_table.add_value("LevelSet dofs", level_set_dof_handler_size); + convergence_table.evaluate_convergence_rates( + "LevelSet dofs", ConvergenceTable::reduction_rate_log2); + + convergence_table.add_value("Active dofs", dof_handler_size); + convergence_table.evaluate_convergence_rates( + "Active dofs", ConvergenceTable::reduction_rate_log2); + } + + template + void LaplaceBeltramiSolver::initialize_matrices() + { + pcout << "Initializing the matrix... " + << "\n" + << std::flush; + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), + dof_handler.n_dofs(), + locally_relevant_dofs); + constraints.reinit(locally_owned_dofs, locally_relevant_dofs); + + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + + SparsityTools::distribute_sparsity_pattern(dsp, + locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + global_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + // @sect3{Assembly and surface accumulation} + // We use a MeshWorker to assemble the linear problem efficiently. + // This cell worker does not do anything for non-intersected cells. + template + void LaplaceBeltramiSolver::assemble_system() + { + pcout << "Assembling... " + << "\n" + << std::flush; + TimerOutput::Scope t(computing_timer, "assembly"); + + const auto cell_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch_data, + CopyData ©_data) { + if (mesh_classifier.location_to_level_set(cell) == + NonMatching::LocationToLevelSet::intersected && + cell->is_locally_owned()) + { + // Once we know that the cell is intersected, we construct the + // unfitted quadratures for the solutions FE space on the cell. + scratch_data.non_matching_fe_values.reinit(cell); + copy_data.reinit(cell); + copy_data.cell_matrix = 0; + copy_data.cell_rhs = 0; + const std::optional> + &surface_fe_values = + scratch_data.non_matching_fe_values.get_surface_fe_values(); + const std::vector &JxW_surface = + surface_fe_values->get_JxW_values(); + + // The accumulation of the surface integrals, including the forcing, + // is performed here. + for (unsigned int q : surface_fe_values->quadrature_point_indices()) + { + const Tensor<1, dim> &normal = + surface_fe_values->normal_vector(q); + + for (const unsigned int i : surface_fe_values->dof_indices()) + { + copy_data.cell_rhs(i) += + surface_fe_values->shape_value(i, q) * + right_hand_side.value( + surface_fe_values->quadrature_point(q)) * + JxW_surface[q]; + + for (const unsigned int j : + surface_fe_values->dof_indices()) + { + copy_data.cell_matrix(i, j) += + (surface_fe_values->shape_value(i, q) * + surface_fe_values->shape_value(j, q)) * + JxW_surface[q]; + copy_data.cell_matrix(i, j) += + (surface_fe_values->shape_grad(i, q) - + (normal * surface_fe_values->shape_grad(i, q)) * + normal) * + (surface_fe_values->shape_grad(j, q) - + (normal * surface_fe_values->shape_grad(j, q)) * + normal) * + JxW_surface[q]; + } + } + } + + // The normal-gradient volume stabilization form needs a bulk cell + // integration while other types of stabilization may need face + // quadratures, for example. So we check it first. + // The cell was provided by the solution's DoFHandler, + // so we recast it as a level set's DoFHandler cell. + // However, it is the same geometric entity of the common + // triangulation. + if (stabilization_scheme.needs_cell_worker()) + { + typename DoFHandler::active_cell_iterator level_set_cell = + cell->as_dof_handler_iterator(level_set_dof_handler); + scratch_data.fe_values.reinit(cell); + scratch_data.level_set_fe_values.reinit(level_set_cell); + stabilization_scheme.assemble_cell_worker(level_set, + cell, + scratch_data, + copy_data); + } + } + }; + + // Next, the copier worker distributes the local contributions from + // the CopyData taking into account the constraints. Finally, the + // MeshWorker goes over all cells provided by the solutions' + // DoFHandler. Note that this includes non-intersected cells as + // well, but the cell worker does nothing on them. + const auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_matrix, + c.cell_rhs, + c.local_dof_indices, + global_matrix, + global_rhs); + }; + + ScratchData scratch_data(mapping, + fe_collection, + mesh_classifier, + level_set_dof_handler, + level_set, + surface_update_flags, + cell_quadrature, + quadrature_1D); + + CopyData copy_data; + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + copy_data, + MeshWorker::assemble_own_cells); + + global_matrix.compress(VectorOperation::add); + global_rhs.compress(VectorOperation::add); + } + + // In the following, we solve the resulting linear system of equations. We + // either use a direct solver or AMG. + template + void LaplaceBeltramiSolver::solve() + { + TimerOutput::Scope t(computing_timer, "solve"); + bool apply_direct_solver = false; + const double relative_error = 1e-9 * global_rhs.l2_norm(); + unsigned int n_iterations = 0; + if (apply_direct_solver) + { + pcout << "Solving directly... " << '\n' << std::flush; + SolverControl solver_control(100, relative_error); + TrilinosWrappers::SolverDirect::AdditionalData data; + TrilinosWrappers::SolverDirect trilinos(solver_control, data); + trilinos.solve(global_matrix, + completely_distributed_solution, + global_rhs); + } + else + { + Timer timer; + pcout << "Solving with AMG... " + << "\n" + << std::flush; + const unsigned int max_iterations = 500; + SolverControl solver_control(max_iterations, relative_error); + std::vector> constant_modes; + DoFTools::extract_constant_modes(dof_handler, + ComponentMask(), + constant_modes); + TrilinosWrappers::PreconditionAMG preconditioner_stiffness; + TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data; + Amg_data.constant_modes = constant_modes; + Amg_data.elliptic = true; + Amg_data.higher_order_elements = false; + Amg_data.smoother_sweeps = 2; + Amg_data.aggregation_threshold = 0.02; + Amg_data.output_details = true; + preconditioner_stiffness.initialize(global_matrix); + + SolverCG cg(solver_control); + cg.solve(global_matrix, + completely_distributed_solution, + global_rhs, + preconditioner_stiffness); + n_iterations = solver_control.last_step(); + } + constraints.distribute(completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; + + convergence_table.add_value("Iterations", n_iterations); + } + + // Similarly to what we do in the assembly() function, + // a MeshWorker is used to accumulate errors + // including the stabilization term. At the end, we collect the results, + // and print them out. + template + void LaplaceBeltramiSolver::compute_errors() + { + pcout << "Evaluating errors on the surface..." + << "\n" + << std::flush; + TimerOutput::Scope t(computing_timer, "eval_errors"); + double error_L2_sqr = 0.0; + double error_H1_sqr = 0.0; + double error_stab_sqr = 0.0; + const auto cell_worker = [&](const auto &cell, + auto &scratch_data, + auto ©_data) { + if (mesh_classifier.location_to_level_set(cell) == + NonMatching::LocationToLevelSet::intersected && + cell->is_locally_owned()) + { + double cell_L2_error_sqr = 0.0; + double cell_H1_error_sqr = 0.0; + + copy_data.reinit(cell); + scratch_data.non_matching_fe_values.reinit(cell); + + const std::optional> + &surface_fe_values = + scratch_data.non_matching_fe_values.get_surface_fe_values(); + const std::vector &JxW_surface = + surface_fe_values->get_JxW_values(); + const unsigned int n_q_points = + surface_fe_values->n_quadrature_points; + + std::vector sol(n_q_points); + surface_fe_values->get_function_values(locally_relevant_solution, + sol); + + std::vector> sol_grad(n_q_points); + surface_fe_values->get_function_gradients(locally_relevant_solution, + sol_grad); + + 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); + const double error_at_point = + sol.at(q) - analytical_solution.value(point); + const Tensor<1, dim> grad_error_at_point = + (sol_grad.at(q) - (normal * sol_grad.at(q)) * normal - + analytical_solution.gradient(point)); + + cell_L2_error_sqr += + Utilities::pow(error_at_point, 2) * JxW_surface[q]; + cell_H1_error_sqr += + grad_error_at_point * grad_error_at_point * JxW_surface[q]; + } + copy_data.cell_L2_error_sqr = cell_L2_error_sqr; + copy_data.cell_H1_error_sqr = cell_H1_error_sqr; + + if (stabilization_scheme.needs_cell_worker()) + { + typename DoFHandler::active_cell_iterator level_set_cell = + cell->as_dof_handler_iterator(level_set_dof_handler); + scratch_data.fe_values.reinit(cell); + scratch_data.level_set_fe_values.reinit(level_set_cell); + stabilization_scheme.evaluate_cell_worker( + locally_relevant_solution, + level_set, + cell, + scratch_data, + copy_data); + } + } + }; + + const auto copier = [&](const auto ©_data) { + if (copy_data.cell_index < active_fe_indicator.size()) + { + error_L2_sqr += copy_data.cell_L2_error_sqr; + error_H1_sqr += copy_data.cell_H1_error_sqr; + error_stab_sqr += copy_data.cell_stab_sqr; + } + }; + + ScratchData scratch_data(mapping, + fe_collection, + mesh_classifier, + level_set_dof_handler, + level_set, + surface_update_flags, + cell_quadrature, + quadrature_1D); + + CopyDataError copy_data; + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + copy_data, + MeshWorker::assemble_own_cells); + + const double error_L2 = + std::sqrt(Utilities::MPI::sum(error_L2_sqr, mpi_communicator)); + const double error_semiH1 = + std::sqrt(Utilities::MPI::sum(error_H1_sqr, mpi_communicator)); + const double error_stab = + std::sqrt(Utilities::MPI::sum(error_stab_sqr, mpi_communicator)); + + 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); + + convergence_table.add_value("H1 error", error_semiH1); + convergence_table.evaluate_convergence_rates( + "H1 error", ConvergenceTable::reduction_rate_log2); + convergence_table.set_scientific("H1 error", true); + + convergence_table.add_value("Stab norm", error_stab); + convergence_table.evaluate_convergence_rates( + "Stab norm", ConvergenceTable::reduction_rate_log2); + convergence_table.set_scientific("Stab norm", true); + } + + // The following two methods perform VTK output of the preliminary mesh + // refinements for the geometry approximation and of the TraceFEM solution. + // The important difference between the two is that the non-intersected cells + // are excluded from the output saving considerable amount of time and + // storage. + template + void LaplaceBeltramiSolver::output_level_set(const unsigned int cycle) + { + pcout << "Writing vtu file for surface... " << '\n' << std::flush; + TimerOutput::Scope t(computing_timer, "output_level_set"); + DataOut data_out; + data_out.add_data_vector(level_set_dof_handler, level_set, "level_set"); + data_out.add_data_vector(active_fe_indicator, "ref_indicator"); + data_out.build_patches(); + + data_out.write_vtu_in_parallel("surface_" + std::to_string(cycle) + ".vtu", + mpi_communicator); + } + + template + void LaplaceBeltramiSolver::output_solution() + { + pcout << "Writing vtu file... " << std::flush; + TimerOutput::Scope t(computing_timer, "output_solution"); + Vector exact(dof_handler.locally_owned_dofs().size()); + + VectorTools::interpolate(dof_handler, analytical_solution, exact); + DataOut data_out; + data_out.add_data_vector(dof_handler, + locally_relevant_solution, + "solution"); + data_out.add_data_vector(dof_handler, exact, "exact"); + 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() && cell->is_locally_owned() && + mesh_classifier.location_to_level_set(cell) == + NonMatching::LocationToLevelSet::intersected; + }); + data_out.build_patches(); + + data_out.write_vtu_in_parallel("solution.vtu", mpi_communicator); + } + + // The method localize_surface() generates iteratively a surface approximation + // as described above. Once the surface approximation is constructed, the main + // logic of the solver is executed as presented in the method run(). + template + void LaplaceBeltramiSolver::localize_surface() + { + unsigned int preliminary_levels = 3; + for (unsigned int localization_cycle = 0; + localization_cycle < preliminary_levels; + ++localization_cycle) + { + pcout << std::endl + << "Preliminary refinement #" << localization_cycle << std::endl; + setup_discrete_level_set(); + mark_intersected(); + output_level_set(localization_cycle); + refine_grid(); + } + computing_timer.reset(); + } + + template + void LaplaceBeltramiSolver::run() + { + make_grid(); + localize_surface(); + const unsigned int convergence_levels = 3; + for (unsigned int cycle = 0; cycle < convergence_levels; ++cycle) + { + pcout << std::endl << "Convergence refinement #" << cycle << std::endl; + setup_discrete_level_set(); + distribute_dofs(); + initialize_matrices(); + assemble_system(); + solve(); + compute_errors(); + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + convergence_table.write_text(pcout.get_stream()); + + computing_timer.print_summary(); + computing_timer.reset(); + if (cycle < convergence_levels - 1) + { + mark_intersected(); + refine_grid(); + } + else + output_solution(); + + computing_timer.print_summary(); + computing_timer.reset(); + } + } +} // namespace Step90 + +int main(int argc, char *argv[]) +{ + try + { + using namespace dealii; + using namespace Step90; + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + LaplaceBeltramiSolver<3> LB_solver; + LB_solver.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}