Major revision of the tutorial after review.

doc/doxygen/references.bib

@@ -1710,6 +1710,20 @@ @article{cutfem_2015
   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

doc/doxygen/tutorial/tutorial.h.in

@@ -697,6 +697,11 @@
  *       <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>
  *
  *

examples/step-90/doc/intro.dox

@@ -1,95 +1,141 @@
 <i>
-This program was contributed by Vladimir Yushutin.
+This program was contributed by Vladimir Yushutin and Timo Heister, Clemson University, 2023.
 
-The material is based upon work partially supported by NSF.
+This material is based upon work partly supported by the National
+Science Foundation Award DMS-
 </i>
 
-<a name="Intro to step-90"></a>
+<a name="step-90-Intro"></a>
 <h1>Introduction</h1>
 
-In this example, we implement the trace finite element method (TraceFEM) in deal.II.
+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, but both should be based on the same bulk triangulation.
- Second, we make sure that 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.
+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>
-One of the possible bottlenecks in achieving the two-dimensional complexity while solving a surface PDE with TraceFEM
-is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create
-a mesh which is localized near the zero contour line, i.e near the surface, to restore the two-dimensional complexity.
-To demonstrate this 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}
- We note that the surface's curvature varies significantly.
+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, Laplace--Beltrami equation,
+We would like to solve the simplest possible problem  defined on a surface, namely the Laplace--Beltrami equation,
 @f{equation*}
-  -\Delta_\Gamma u + u = f \qquad  \text{in }\, \Gamma,
+ -\Delta_\Gamma u + c u = f \qquad  \text{in }\, \Gamma,
 @f}
- Here we added the term $u$ to the left-hand side so the problem becomes well-posed.
- The force $f$ is calculated analytically for a manufactured solution,
- @f{equation*}
- u=xy\,,
-  @f}
-using the exact normal $\mathbf{n}$ and the exact Hessian $\nabla^2\mathbf{n}$ of the surface computed analytically.
+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>The Trace Finite Element Method</h3>
+<h3>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.
 
-TraceFEM is an unfitted method meaning 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 an approximation, denoted by $\Gamma_h$,
-of $\Gamma$ in a form of quadrature points for the surface integration combined with the approximate surface normals.
+<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 the $\Gamma_h$ and that consists of the intersected cells $\mathcal{T}_\Omega^h$,
+that completely covers $\Gamma_h$ and that consists of the intersected cells $\mathcal{T}_\Gamma^h$,
 @f{equation*}
-  \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Gamma_h \neq \emptyset \}.
+ \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}_\Omega^h \},
+ 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}_\Omega^h} \overline{T}$.
+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}_\Omega^h$ or not,
+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.
-For this purpose, we will use the class NonMatching::MeshClassifier.
+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 bilinear forms
+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}.
+ 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 of $f$ from $\Gamma$ to $\Omega_h$.
+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
+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,
+ 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},
+ 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 worst case.
+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 \}.
+ \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\|$. This and the approximate quadrature rules for the integration over $\Gamma_h$ are often referred to as the "geometrical error".
-Luckily, one can show 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 us to independently choose a more accurate geometry representation
-with a higher-order level set function, compared to the function space for solving the Laplace--Beltrami equation.
+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.

examples/step-90/doc/results.dox

@@ -1,26 +1,32 @@
 <h1>Results</h1>
 
-The numerical solution for a very finest mesh is shown below.
+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
+@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 is shown in the table below.
+The results of the convergence study are shown in the following table.
 
-| Cycle |   DOFS   |   EOC | Iterations | $L^2$-Error |   EOC | $H^1$-Error |   EOC |Stabilization|   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  |
+| 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 term $s_h(u_h, u_h)$.
+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>
-In progress...
+<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%

examples/step-90/doc/tooltip

@@ -1 +1 @@
-Solving Laplace-Beltrami equation using the trace finite element method.
+Solving the Laplace-Beltrami equation on a surface using the trace finite element method.