Documentation and code structure revisited

examples/step-90/doc/intro.dox

@@ -7,83 +7,89 @@ The material is based upon work partially supported by NSF.
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
-<h3>The Trace Finite Element Method</h3>
-
 In this example, we implement the trace finite element method (TraceFEM) in deal.II.
 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 $Q_1$ solution, but both are based on the same triangulation. Second, we make sure
- that performance of TraceFEM corresponds to that of a classical fitted FEM for a two-dimensional problem.
- Both goals are achieved by using a combination of MeshWorkers, adaptive refinement and NonMatching capabilites,
-  and the example demonstrates the non-trivial details.
+ 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.
+
+<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.
+
+<h3>Model problem</h3>
 
-We would like to solve the simplest possible problem, Laplace--Beltrami equation, posed on a unit sphere $\Gamma$:
+We would like to solve the simplest possible problem, Laplace--Beltrami equation,
 @f{equation*}
-  -\Delta u + u &= f \qquad && \text{in }\, \Gamma,
+  -\Delta_\Gamma u + u = f \qquad  \text{in }\, \Gamma,
 @f}
-where we choose $f(x) = \sin\theta$. TraceFEM is an unfitted method meaning the surface is immersed into the background mesh.
-To solve this problem,
-we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$,
-that completely covers the $\Gamma$:
+ 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.
+
+<h3>The Trace Finite Element Method</h3>
+
+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.
+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$,
 @f{equation*}
-  \mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Gamma \neq \emptyset \}.
+  \mathcal{T}_\Omega^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(\mathcal{N}_\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}_\Omega^h \},
 @f}
-where
-@f{equation*}
-  \mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T},
-@f}
-and $\overline{T}$ denotes the closure of $T$.
+where $\Omega_h$ is the union of all intersected cells from $\bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T}$.
 
-This leads to the following weak formulation:
-Find $u_h \in V_\Omega^h$ such that
-@f{equation*}
-  a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
-@f}
-where
+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,
+we set the active_fe_index to either 0 or 1.
+For this purpose, we will use the class NonMatching::MeshClassifier.
+
+A natural candidate for a weak formulation  involves the following bilinear 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}
-  \\
-  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$.
 
-Thus, on each cell intersected by $\Gamma$,
-we need special quadrature rules that only integrate over these parts of the cell,
-that is, over $T \cap \Omega$ and $T \cap \Gamma$.
-
-This leads to the stabilized cut finite element method,
-which reads: Find $u_h \in V_h$ such that
-@f{equation*}
-  A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
-@f}
-where
+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) = a_h(u_h,v_h) + s_h(u_h, v_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}
-
-The normal-gradient stabilization acts on cells and reads
+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}(\bn_h\cdot\nabla u_h, \bn_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.
 
 <h3>Discrete Level Set Function</h3>
-This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function:
+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{dim}} : \psi_h(x) = 0 \}.
+  \Gamma_h &= \{x \in \mathbb{R}^{\text{3}} : \psi_h(x) = 0 \}.
 @f}
-This is often referred to as the "geometrical error".
-However, when the same element order, $p$, is used in $V^h$,
-one can often show that the method gives the same order of convergence
+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.
-
-<h3>The MeshClassifier Class</h3>
-To create $V_\Omega^h$, we first add an FE_Q and an
-FE_Nothing element to an hp::FECollection. We then iterate over each cell,
-$T$, and depending on whether $T$ belongs to $\mathcal{T}_\Omega^h$ or not,
-we set the active_fe_index to either 0 or 1.
-For this purpose, we will use the class NonMatching::MeshClassifier.

examples/step-90/doc/results.dox

@@ -1,20 +1,26 @@
 <h1>Results</h1>
 
-The numerical solution for one of the refinements is shown in the below figure.
-The zero-contour of the level set function is shown as a white line.
-On the intersected cells,
-we see that the numerical solution has a value also outside $\overline{\Omega}$.
-As mentioned earlier, this extension of the solution is artificial.
+The numerical solution for a very finest mesh is shown below.
+
+@image html step-90-solution.png
+
+<h3>Convergence test </h3>
 
 The results of the convergence study is shown in the table below.
-We see that the $L^2$ error decreases as we refine and that the estimated
-order of convergence, EOC, is close to 2.
 
-@image html step-85-solution.png
+| Cycle |   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  |
+
+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)$.
 
-| Cycle | Mesh size | $L^2$-Error | EOC  |
-|:-----:|:---------:|:-----------:|:----:|
-|   0   |  0.3025   |  8.0657e-02 |  -   |
-|   1   |  0.1513   |  1.8711e-02 | 2.11 |
-|   2   |  0.0756   |  4.1624e-03 | 2.17 |
-|   3   |  0.0378   |  9.3979e-04 | 2.15 |
+<h3>Parallel scalability</h3>
+In progress...