![]() |
deal.II version GIT relicensing-2652-g41e7496fd4 2025-02-17 19:10:00+00:00
|
This tutorial depends on step-85.
Table of contents | |
---|---|
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-2028346, OAC-2015848, EAR-1925575, and by the Computational Infrastructure in Geodynamics initiative (CIG), through the NSF under Award EAR-0949446, EAR-1550901, EAR-2149126 via the University of California – Davis.
In this tutorial, we implement the trace finite element method (TraceFEM) in deal.II. TraceFEM solves PDEs posed on a possibly evolving \((dim-1)\)-dimensional surface \(\Gamma\) employing a fixed uniform background mesh of a \(dim\)-dimensional domain in which the surface is embedded. Such surface PDEs arise in problems involving material films with complex properties and in other situations in which a non-trivial condition is imposed on either a stationary or a moving interface. Here we consider a steady, complex, non-trivial surface and the prototypical Laplace-Beltrami equation which is a counterpart of the Poisson problem on flat domains.
Being an unfitted method, TraceFEM allows to circumvent the need of remeshing of an evolving surface if it is implicitly given by the zero contour of a level-set function. At the same time, it easily provides with an extension of the discrete solution to a neighborhood of the surface which turns out to be very handy in case of non-stationary interfaces and films. Certainly, this flexibility comes with a price: one needs to design the nodes and weights for a quadrature customized for each implicit intersection of the zero level-set and the background mesh. Moreover, these intersections may be of arbitrary shape and size manifesting in the so-called "small cut" problem and requiring addition of a stabilization form that restores well-conditioning of the problem.
Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE, e.g., a \(Q_2\) discrete level-set and a \(Q_1\) solution are possible on the same bulk triangulation. Second, we make sure that the performance of TraceFEM in the parallel implementation corresponds to that of a classical fitted FEM for a two-dimensional problem. We demonstrate how to achieve both goals by using a combination of MeshWorker and NonMatching capabilities.
A natural alternative to TraceFEM in solving surface PDEs is the parametric surface finite element method. The latter method relies on an explicit parametrization of the surface which may be not feasible especially for evolving interfaces with an unknown in advance shape - in this sense, TraceFEM is a technique inspired by the level-set description of interfaces. However, the parametric surface finite element method, when applicable, enjoys many well-known properties of fitted methods on flat domains provided the geometric errors - which a present for both methods - are taken under control.
A 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.
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
\begin{equation*} \frac{x^2}{4}+ y^2 + \frac{4 z^2} {(1 + 0.5 \sin(\pi x))^{2}} = 1 \end{equation*}
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).
We would like to solve the simplest possible problem defined on a surface, namely the Laplace–Beltrami equation,
\begin{equation*} -\Delta_\Gamma u + c u = f \qquad \text{in }\, \Gamma, \end{equation*}
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.
We choose the test solution and the right-hand side forcing as the restriction to \(\Gamma\) of
\begin{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), \end{equation*}
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: 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\),
\begin{equation*} \mathcal{T}_\Gamma^h = \{ T \in \mathcal{T}^{h} : T \cap \Gamma_h \neq \emptyset \}. \end{equation*}
The finite element space where we want to find our numerical solution, \(u_h\), is now
\begin{equation*} V_h = \{ v \in C(\Omega_h) : v \in Q_p(T), \, T \in \mathcal{T}_\Gamma^h \}, \end{equation*}
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
\begin{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}. \end{align*}
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
\begin{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. \end{equation*}
Here the normal-gradient stabilization \(s_h\) involves the three-dimensional integration over whole (but intersected) cells and is given by
\begin{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}, \end{equation*}
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 [171].
In TraceFEM we construct the approximation \(\Gamma_h\) using the interpolant \(\psi_h\) of the exact level-set function on the bulk triangulation:
\begin{align*} \Gamma_h &= \{x \in \mathbb{R}^{\text{3}} : \psi_h(x) = 0 \}. \end{align*}
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 [171] 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.
The parallelization 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.
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.
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.
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\):
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.
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.
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.
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\).
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.
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.
The surface of interest corresponds to the zero contour of the following exact level set function:
The manufactured solution to the Laplace–Beltrami problem and the corresponding right-hand side:
There is a single triangulation which is shared by the discretizations of the solution and of the level set.
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:
Since we will adaptively refine the bulk triangulation, two constraints are needed: one for the solution space and another for the level set space.
Discrete vectors initialized with dof_handler and level_set_dof_handler.
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:
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.
Any TraceFEM needs a stabilization, and we choose the normal-gradient, volume stabilization.
Discrete right-hand side and the final matrix corresponding to dof_handler.
Depending on the type of the quadrature, surface, face or volume, we need to define different update flags.
The following variables are used to display the results of the convergence test:
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.
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.
The method fills in the indicator telling which cells are intersected. It is used in the adaptive refinement near the surface.
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.
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.
We use a MeshWorker to assemble the linear problem efficiently. This cell worker does not do anything for non-intersected cells.
Once we know that the cell is intersected, we construct the unfitted quadratures for the solutions FE space on the cell.
The accumulation of the surface integrals, including the forcing, is performed here.
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.
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.
In the following, we solve the resulting linear system of equations. We either use a direct solver or AMG.
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.
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.
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().
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\).
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.
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\).
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.