]> https://gitweb.dealii.org/ - dealii.git/commitdiff
The tutorial step-90 demonstrates an efficient implementation of the stabilized Trace... 16878/head
authorVladimir Yushutin <vyushut@clemson.edu>
Wed, 4 Oct 2023 21:20:25 +0000 (17:20 -0400)
committerVladimir Yushutin <vyushut@clemson.edu>
Mon, 10 Jun 2024 19:35:52 +0000 (15:35 -0400)
14 files changed:
doc/doxygen/images/step-90-solution.png [new file with mode: 0644]
doc/doxygen/images/step-90_mesh_cut.png [new file with mode: 0644]
doc/doxygen/images/step-90_prelim.png [new file with mode: 0644]
doc/doxygen/images/step-90_surface.png [new file with mode: 0644]
doc/doxygen/images/step-90_weak-vs-strong.png [new file with mode: 0644]
doc/doxygen/references.bib
doc/doxygen/tutorial/tutorial.h.in
examples/step-90/CMakeLists.txt [new file with mode: 0644]
examples/step-90/doc/builds-on [new file with mode: 0644]
examples/step-90/doc/intro.dox [new file with mode: 0644]
examples/step-90/doc/kind [new file with mode: 0644]
examples/step-90/doc/results.dox [new file with mode: 0644]
examples/step-90/doc/tooltip [new file with mode: 0644]
examples/step-90/step-90.cc [new file with mode: 0644]

diff --git a/doc/doxygen/images/step-90-solution.png b/doc/doxygen/images/step-90-solution.png
new file mode 100644 (file)
index 0000000..0f90965
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 (file)
index 0000000..5367d05
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 (file)
index 0000000..0b74116
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 (file)
index 0000000..0b9fb8a
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 (file)
index 0000000..0de7c87
Binary files /dev/null and b/doc/doxygen/images/step-90_weak-vs-strong.png differ
index 85c081c07751d4c2db88adb76e609cd0be106aec..cc963dd0b28cda3396bb591f1355dd6675b23214 100644 (file)
   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
index 89351c0cd5c6350c134d0a23e6a49d7528e51a12..d4b329efa90886c36590ed8baf513ce85ca24c4d 100644 (file)
  *       <br/> Keywords: FERemoteEvaluation
  *       </td></tr>
  *
+ *   <tr valign="top">
+ *       <td>step-90</td>
+ *       <td> Solving the Laplace-Beltrami equation on a surface using the trace finite element method.
+ *       <br/> Keywords: MeshWorker::mesh_loop(), NonMatching::FEImmersedSurfaceValues
+ *       </td></tr>
  * </table>
  *
  *
diff --git a/examples/step-90/CMakeLists.txt b/examples/step-90/CMakeLists.txt
new file mode 100644 (file)
index 0000000..8ead29f
--- /dev/null
@@ -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 (file)
index 0000000..b1e52c0
--- /dev/null
@@ -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 (file)
index 0000000..4318c4e
--- /dev/null
@@ -0,0 +1,141 @@
+<i>
+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-
+</i>
+
+<a name="step-90-Intro"></a>
+<h1>Introduction</h1>
+
+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.
+
+
+<h3>A non-trivial surface</h3>
+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%
+
+
+<h3>Model problem</h3>
+
+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.
+
+<h3>Manufactured exact solution</h3>
+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.
+
+<h3>The Trace Finite Element Method</h3>
+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.
+
+<h3>Discrete Level Set Function</h3>
+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 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-90/doc/results.dox b/examples/step-90/doc/results.dox
new file mode 100644 (file)
index 0000000..bbb9546
--- /dev/null
@@ -0,0 +1,37 @@
+<h1>Results</h1>
+
+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%
+
+<h3>Convergence test </h3>
+
+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$.
+
+<h3>Parallel scalability</h3>
+
+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 (file)
index 0000000..3f5c006
--- /dev/null
@@ -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 (file)
index 0000000..f97e9c9
--- /dev/null
@@ -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 <deal.II/base/convergence_table.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+#include <deal.II/non_matching/fe_immersed_values.h>
+#include <deal.II/non_matching/fe_values.h>
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+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 <int dim>
+  class TamarindShape : public Function<dim>
+  {
+  public:
+    TamarindShape()
+      : Function<dim>()
+    {}
+    double value(const Point<dim>  &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<dim>  &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<dim>  &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 <int dim>
+  class AnalyticalSolution : public Function<dim>
+  {
+  private:
+    const TamarindShape<dim> tamarind;
+
+  public:
+    AnalyticalSolution()
+      : Function<dim>()
+    {}
+    double value(const Point<dim>  &point,
+                 const unsigned int component = 0) const override;
+
+    Tensor<1, dim> gradient(const Point<dim>  &point,
+                            const unsigned int component = 0) const override;
+  };
+
+  template <int dim>
+  double AnalyticalSolution<dim>::value(const Point<dim>  &point,
+                                        const unsigned int component) const
+  {
+    AssertIndexRange(component, this->n_components);
+    (void)component;
+    return point[0] * point[1];
+  }
+
+  template <int dim>
+  Tensor<1, dim>
+  AnalyticalSolution<dim>::gradient(const Point<dim>  &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 <int dim>
+  class RightHandSide : public Function<dim>
+  {
+    const TamarindShape<dim> tamarind;
+
+  public:
+    RightHandSide()
+      : Function<dim>()
+    {}
+
+    virtual double value(const Point<dim>  &p,
+                         const unsigned int component = 0) const override;
+  };
+
+  template <int dim>
+  double RightHandSide<dim>::value(const Point<dim>  &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 <int dim>
+  struct ScratchData
+  {
+    ScratchData(const Mapping<dim>                     &mapping,
+                const hp::FECollection<dim>            &fe_collection,
+                const NonMatching::MeshClassifier<dim> &mesh_classifier,
+                const DoFHandler<dim>                  &level_set_dof_handler,
+                const VectorType                       &level_set,
+                const NonMatching::RegionUpdateFlags nonmatching_update_flags,
+                const Quadrature<dim>               &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<types::fe_index>(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<dim> &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<dim>                           fe_values;
+    const NonMatching::RegionUpdateFlags    region_update_flags;
+    const Quadrature<1>                    &quadrature_1D;
+    const hp::FECollection<dim>            &fe_collection;
+    const NonMatching::MeshClassifier<dim> &mesh_classifier;
+    const DoFHandler<dim>                  &level_set_dof_handler;
+    const VectorType                       &level_set;
+    FEValues<dim>                           level_set_fe_values;
+    NonMatching::FEValues<dim>              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 <int dim>
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+
+    void reinit(const typename DoFHandler<dim>::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 <int dim>
+  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<dim>::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 <int dim>
+  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<dim>::active_cell_iterator &cell,
+      ScratchData<dim>                                     &scratch_data,
+      CopyData<dim>                                        &copy_data) const
+    {
+      const FEValues<dim> &fe_values = scratch_data.fe_values;
+      const FEValues<dim> &level_set_fe_values =
+        scratch_data.level_set_fe_values;
+
+      const std::vector<double> &JxW_cell = fe_values.get_JxW_values();
+
+      std::vector<Tensor<1, dim>> 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<dim>::active_cell_iterator &cell,
+      ScratchData<dim>                                     &scratch_data,
+      CopyDataError<dim>                                   &copy_data) const
+    {
+      double                     cell_stab_sqr = 0.0;
+      const FEValues<dim>       &fe_values     = scratch_data.fe_values;
+      const std::vector<double> &JxW_cell      = fe_values.get_JxW_values();
+      const unsigned int n_q_points = fe_values.get_quadrature_points().size();
+      const FEValues<dim> &level_set_fe_values =
+        scratch_data.level_set_fe_values;
+
+      std::vector<Tensor<1, dim>> level_set_grad(n_q_points);
+      level_set_fe_values.get_function_gradients(level_set, level_set_grad);
+
+      std::vector<Tensor<1, dim>> 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 <int dim>
+  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<dim> tamarind;
+
+    // The manufactured solution to the Laplace--Beltrami problem and the
+    // corresponding right-hand side:
+    const AnalyticalSolution<dim> analytical_solution;
+    const RightHandSide<dim>      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<dim, dim> 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<dim> fe_collection;
+    DoFHandler<dim>       dof_handler;
+
+    const unsigned int level_set_fe_degree;
+    const FE_Q<dim>    level_set_fe;
+    DoFHandler<dim>    level_set_dof_handler;
+
+    const MappingQ1<dim> 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<double> constraints;
+    AffineConstraints<double> 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<float> 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<dim> 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<dim> cell_quadrature;
+    const QGauss<1>   quadrature_1D;
+
+    // Any TraceFEM needs a stabilization, and we choose the normal-gradient,
+    // volume stabilization.
+    const NormalGradientVolumeStabilization<dim> 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 <int dim>
+  LaplaceBeltramiSolver<dim>::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<dim>(fe_degree), FE_Nothing<dim>())
+    , 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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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<types::fe_index>(ActiveFEIndex::lagrange));
+        else
+          cell->set_active_fe_index(
+            static_cast<types::fe_index>(ActiveFEIndex::nothing));
+      }
+  }
+
+  // The method fills in the indicator telling which cells are intersected. It
+  // is used in the adaptive refinement near the surface.
+  template <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::assemble_system()
+  {
+    pcout << "Assembling... "
+          << "\n"
+          << std::flush;
+    TimerOutput::Scope t(computing_timer, "assembly");
+
+    const auto cell_worker =
+      [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+          ScratchData<dim>                                     &scratch_data,
+          CopyData<dim>                                        &copy_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<NonMatching::FEImmersedSurfaceValues<dim>>
+              &surface_fe_values =
+                scratch_data.non_matching_fe_values.get_surface_fe_values();
+            const std::vector<double> &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<dim>::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<dim> &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             global_matrix,
+                                             global_rhs);
+    };
+
+    ScratchData<dim> scratch_data(mapping,
+                                  fe_collection,
+                                  mesh_classifier,
+                                  level_set_dof_handler,
+                                  level_set,
+                                  surface_update_flags,
+                                  cell_quadrature,
+                                  quadrature_1D);
+
+    CopyData<dim> 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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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<std::vector<bool>> 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<VectorType> 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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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       &copy_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<NonMatching::FEImmersedSurfaceValues<dim>>
+            &surface_fe_values =
+              scratch_data.non_matching_fe_values.get_surface_fe_values();
+          const std::vector<double> &JxW_surface =
+            surface_fe_values->get_JxW_values();
+          const unsigned int n_q_points =
+            surface_fe_values->n_quadrature_points;
+
+          std::vector<double> sol(n_q_points);
+          surface_fe_values->get_function_values(locally_relevant_solution,
+                                                 sol);
+
+          std::vector<Tensor<1, dim>> 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<dim> &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<dim>::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 &copy_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<dim> scratch_data(mapping,
+                                  fe_collection,
+                                  mesh_classifier,
+                                  level_set_dof_handler,
+                                  level_set,
+                                  surface_update_flags,
+                                  cell_quadrature,
+                                  quadrature_1D);
+
+    CopyDataError<dim> 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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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<dim>       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 <int dim>
+  void LaplaceBeltramiSolver<dim>::output_solution()
+  {
+    pcout << "Writing vtu file... " << std::flush;
+    TimerOutput::Scope t(computing_timer, "output_solution");
+    Vector<double>     exact(dof_handler.locally_owned_dofs().size());
+
+    VectorTools::interpolate(dof_handler, analytical_solution, exact);
+    DataOut<dim> 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<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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 <int dim>
+  void LaplaceBeltramiSolver<dim>::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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.