 Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
The step-60 tutorial program

This tutorial depends on step-6.

1. Introduction
2. The commented program
1. Results
2. The plain program

This program was contributed by Luca Heltai and Giovanni Alzetta, SISSA, Trieste.

Note
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: # Introduction

### Non-matching grid constraints through distributed Lagrange multipliers

In this tutorial we consider the case of two domains, $$\Omega$$ in $$R^{\text{spacedim}}$$ and $$\Gamma$$ in $$R^{\text{dim}}$$, where $$\Gamma$$ is embedded in $$\Omega$$ ( $$\Gamma \subseteq \Omega$$). We want to solve a partial differential equation on $$\Omega$$, enforcing some conditions on the solution of the problem on the embedded domain $$\Gamma$$.

There are two interesting scenarios:

• the geometrical dimension dim of the embedded domain $$\Gamma$$ is the same of the domain $$\Omega$$ (spacedim), that is, the spacedim-dimensional measure of $$\Gamma$$ is not zero, or
• the embedded domain $$\Gamma$$ has an intrinsic dimension dim which is smaller than that of $$\Omega$$ (spacedim), thus its spacedim-dimensional measure is zero; for example it is a curve embedded in a two dimensional domain, or a surface embedded in a three-dimensional domain.

In both cases define the restriction operator $$\gamma$$ as the operator that, given a continuous function on $$\Omega$$, returns its (continuous) restriction on $$\Gamma$$, i.e.,

$\gamma : C^0(\Omega) \mapsto C^0(\Gamma), \quad \text{ s.t. } \gamma u = u|_{\Gamma} \in C^0(\Gamma), \quad \forall u \in C^0(\Omega).$

It is well known that the operator $$\gamma$$ can be extended to a continuous operator on $$H^1(\Omega)$$, mapping functions in $$H^1(\Omega)$$ to functions in $$H^1(\Gamma)$$ when the intrinsic dimension of $$\Gamma$$ is the same of $$\Omega$$.

The same is true, with a less regular range space (namely $$H^{1/2}(\Gamma)$$), when the dimension of $$\Gamma$$ is one less with respect to $$\Omega$$, and $$\Gamma$$ does not have a boundary. In this second case, the operator $$\gamma$$ is also known as the trace operator, and it is well defined for Lipschitz co-dimension one curves and surfaces $$\Gamma$$ embedded in $$\Omega$$ (read this wikipedia article for further details on the trace operator).

The co-dimension two case is a little more complicated, and in general it is not possible to construct a continuous trace operator, not even from $$H^1(\Omega)$$ to $$L^2(\Gamma)$$, when the dimension of $$\Gamma$$ is zero or one respectively in two and three dimensions.

In this tutorial program we're not interested in further details on $$\gamma$$: we take the extension $$\gamma$$ for granted, assuming that the dimension of the embedded domain (dim) is always smaller by one or equal with respect to the dimension of the embedding domain $$\Omega$$ (spacedim).

We are going to solve the following differential problem: given a sufficiently regular function $$g$$ on $$\Gamma$$, find the solution $$u$$ to

\begin{eqnarray*} - \Delta u + \gamma^T \lambda &=& 0 \text{ in } \Omega\\ \gamma u &=& g \text{ in } \Gamma \\ u & = & 0 \text{ on } \partial\Omega. \end{eqnarray*}

This is a constrained problem, where we are looking for a harmonic function $$u$$ that satisfies homogeneous boundary conditions on $$\partial\Omega$$, subject to the constraint $$\gamma u = g$$ using a Lagrange multiplier.

This problem has a physical interpretation: harmonic functions, i.e., functions that satisfy the Laplace equation, can be thought of as the displacements of a membrane whose boundary values are prescribed. The current situation then corresponds to finding the shape of a membrane for which not only the displacement at the boundary, but also on $$\Gamma$$ is prescribed. For example, if $$\Gamma$$ is a closed curve in 2d space, then that would model a soap film that is held in place by a wire loop along $$\partial \Omega$$ as well as a second loop along $$\Gamma$$. In cases where $$\Gamma$$ is a whole area, you can think of this as a membrane that is stretched over an obstacle where $$\Gamma$$ is the contact area. (If the contact area is not known we have a different problem – called the "obstacle problem" – which is modeled in step-41.)

As a first example we study the zero Dirichlet boundary condition on $$\partial\Omega$$. The same equations apply if we apply zero Neumann boundary conditions on $$\partial\Omega$$ or a mix of the two.

The variational formulation can be derived by introducing two infinite dimensional spaces $$V(\Omega)$$ and $$Q^*(\Gamma)$$, respectively for the solution $$u$$ and for the Lagrange multiplier $$\lambda$$.

Multiplying the first equation by $$v \in V(\Omega)$$ and the second by $$q \in Q(\Gamma)$$, integrating by parts when possible, and exploiting the boundary conditions on $$\partial\Omega$$, we obtain the following variational problem:

Given a sufficiently regular function $$g$$ on $$\Gamma$$, find the solution $$u$$ to

\begin{eqnarray*} (\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& 0 \qquad \forall v \in V(\Omega) \\ (\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} \qquad \forall q \in Q(\Gamma), \end{eqnarray*}

where $$(\cdot, \cdot)_{\Omega}$$ and $$(\cdot, \cdot)_{\Gamma}$$ represent, respectively, $$L^2$$ scalar products in $$\Omega$$ and in $$\Gamma$$.

Inspection of the variational formulation tells us that the space $$V(\Omega)$$ can be taken to be $$H^1_0(\Omega)$$. The space $$Q(\Gamma)$$, in the co-dimension zero case, should be taken as $$H^1(\Gamma)$$, while in the co-dimension one case should be taken as $$H^{1/2}(\Gamma)$$.

The function $$g$$ should therefore be either in $$H^1(\Gamma)$$ (for the co-dimension zero case) or $$H^{1/2}(\Gamma)$$ (for the co-dimension one case). This leaves us with a Lagrange multiplier $$\lambda$$ in $$Q^*(\Gamma)$$, which is either $$H^{-1}(\Gamma)$$ or $$H^{-1/2}(\Gamma)$$.

There are two options for the discretization of the problem above. One could choose matching discretizations, where the Triangulation for $$\Gamma$$ is aligned with the Triangulation for $$\Omega$$, or one could choose to discretize the two domains in a completely independent way.

The first option is clearly more indicated for the simple problem we proposed above: it is sufficient to use a single Triangulation for $$\Omega$$ and then impose certain constraints depending $$\Gamma$$. An example of this approach is studied in step-40, where the solution has to stay above an obstacle and this is achieved imposing constraints on $$\Omega$$.

To solve more complex problems, for example one where the domain $$\Gamma$$ is time dependent, the second option could be a more viable solution. Handling non aligned meshes is complex by itself: to illustrate how is done we study a simple problem.

The technique we describe here is presented in the literature using one of many names: the immersed finite element method, the fictitious boundary method, the distributed Lagrange multiplier method, and others. The main principle is that the discretization of the two grids and of the two finite element spaces are kept completely independent. This technique is particularly efficient for the simulation of fluid-structure interaction problems, where the configuration of the embedded structure is part of the problem itself, and one solves a (possibly non-linear) elastic problem to determine the (time dependent) configuration of $$\Gamma$$, and a (possibly non-linear) flow problem in $$\Omega \setminus \Gamma$$, plus coupling conditions on the interface between the fluid and the solid.

In this tutorial program we keep things a little simpler, and we assume that the configuration of the embedded domain is given in one of two possible ways:

• as a deformation mapping $$\psi: \Gamma_0 \mapsto \Gamma \subseteq \Omega$$, defined on a continuous finite dimensional space on $$\Gamma_0$$ and representing, for any point $$x \in \Gamma_0$$, its coordinate $$\psi(x)$$ in $$\Omega$$;
• as a displacement mapping $$\delta \psi(x) = \psi(x)-x$$ for $$x\in \Gamma_0$$, representing for any point $$x$$ the displacement vector applied in order to deform $$x$$ to its actual configuration $$\psi(x) = x +\delta\psi(x)$$.

We define the embedded reference domain $$\Gamma_0$$ embedded_grid: on this triangulation we construct a finite dimensional space (embedded_configuration_dh) to describe either the deformation or the displacement through a FiniteElement system of FE_Q objects (embedded_configuration_fe). This finite dimensional space is used only to interpolate a user supplied function (embedded_configuration_function) representing either $$\psi$$ (if the parameter use_displacement is set to false) or $$\delta\psi$$ (if the parameter use_displacement is set to true).

The Lagrange multiplier $$\lambda$$ and the user supplied function $$g$$ are defined through another finite dimensional space embedded_dh, and through another FiniteElement embedded_fe, using the same reference domain. In order to take into account the deformation of the domain, either a MappingFEField or a MappingQEulerian object are initialized with the embedded_configuration vector.

In the embedding space, a standard finite dimensional space space_dh is constructed on the embedding grid space_grid, using the FiniteElement space_fe, following almost verbatim the approach taken in step-6.

We represent the discretizations of the spaces $$V$$ and $$Q$$ with

$V_h(\Omega) = \text{span} \{v_i\}_{i=1}^n$

and

$Q_h(\Gamma) = \text{span} \{q_i\}_{i=1}^m$

respectively, where $$n$$ is the dimension of space_dh, and $$m$$ the dimension of embedded_dh.

Once all the finite dimensional spaces are defined, the variational formulation of the problem above leaves us with the following finite dimensional system of equations:

$\begin{pmatrix} K & C^T \\ C & 0 \end{pmatrix} \begin{pmatrix} u \\ \lambda \end{pmatrix} = \begin{pmatrix} 0 \\ G \end{pmatrix}$

where

\begin{eqnarray*} K_{ij} &\dealcoloneq& (\nabla v_j, \nabla v_i)_\Omega \qquad i,j=1,\dots,n \\ C_{\alpha j} &\dealcoloneq& (v_j, q_\alpha)_\Gamma \qquad j=1,\dots,n, \alpha = 1,\dots, m \\\\ G_{\alpha} &\dealcoloneq& (g, q_\alpha)_\Gamma \qquad \alpha = 1,\dots, m. \end{eqnarray*}

While the matrix $$K$$ is the standard stiffness matrix for the Poisson problem on $$\Omega$$, and the vector $$G$$ is a standard right-hand-side vector for a finite element problem with forcing term $$g$$ on $$\Gamma$$, (see, for example, step-3), the matrix $$C$$ or its transpose $$C^T$$ are non-standard since they couple information on two non-matching grids.

In particular, the integral that appears in the computation of a single entry of $$C$$, is computed on $$\Gamma$$. As usual in finite elements we split this integral into contributions from all cells of the triangulation used to discretize $$\Gamma$$, we transform the integral on $$K$$ to an integral on the reference element $$\hat K$$, where $$F_{K}$$ is the mapping from $$\hat K$$ to $$K$$, and compute the integral on $$\hat K$$ using a quadrature formula:

$C_{\alpha j} \dealcoloneq (v_j, q_\alpha)_\Gamma = \sum_{K\in \Gamma} \int_{\hat K} \hat q_\alpha(\hat x) (v_j \circ F_{K}) (\hat x) J_K (\hat x) \mathrm{d} \hat x = \sum_{K\in \Gamma} \sum_{i=1}^{n_q} \big(\hat q_\alpha(\hat x_i) (v_j \circ F_{K}) (\hat x_i) J_K (\hat x_i) w_i \big)$

Computing this sum is non-trivial because we have to evaluate $$(v_j \circ F_{K}) (\hat x_i)$$. In general, if $$\Gamma$$ and $$\Omega$$ are not aligned, the point $$F_{K}(\hat x_i)$$ is completely arbitrary with respect to $$\Omega$$, and unless we figure out a way to interpolate all basis functions of $$V_h(\Omega)$$ on an arbitrary point on $$\Omega$$, we cannot compute the integral needed for an entry of the matrix $$C$$.

To evaluate $$(v_j \circ F_{K}) (\hat x_i)$$ the following steps needs to be taken (as shown in the picture below):

• For a given cell $$K$$ in $$\Gamma$$ compute the real point $$y_i \dealcoloneq F_{K} (\hat x_i)$$, where $$x_i$$ is one of the quadrature points used for the integral on $$K \subseteq \Gamma$$.
• Find the cell of $$\Omega$$ in which $$y_i$$ lies. We shall call this element $$T$$.
• To evaluate the basis function use the inverse of the mapping $$G_T$$ that transforms the reference element $$\hat T$$ into the element $$T$$: $$v_j(y_i) = \hat v_j \circ G^{-1}_{T} (y_i)$$. The three steps above can be computed by calling, in turn,

This is what the deal.II function VectorTools::point_value() does when evaluating a finite element field (not just a single shape function) at an arbitrary point; but this would be inefficient in this case.

A better solution is to use a convenient wrapper to perform the first three steps on a collection of points: GridTools::compute_point_locations(). If one is actually interested in computing the full coupling matrix, then it is possible to call the method NonMatching::create_coupling_mass_matrix(), that performs the above steps in an efficient way, reusing all possible data structures, and gathering expensive steps together. This is the function we'll be using later in this tutorial.

We solve the final saddle point problem by an iterative solver, applied to the Schur complement $$S$$ (whose construction is described, for example, in step-20), and we construct $$S$$ using LinearOperator classes.

### The testcase

The problem we solve here is identical to step-4, with the difference that we impose some constraints on an embedded domain $$\Gamma$$. The tutorial is written in a dimension independent way, and in the results section we show how to vary both dim and spacedim.

The tutorial is compiled for dim equal to one and spacedim equal to two. If you want to run the program in embedding dimension spacedim equal to three, you will most likely want to change the reference domain for $$\Gamma$$ to be, for example, something you read from file, or a closed sphere that you later deform to something more interesting.

In the default scenario, $$\Gamma$$ has co-dimension one, and this tutorial program implements the Fictitious Boundary Method. As it turns out, the same techniques are used in the Variational Immersed Finite Element Method, and the coupling operator $$C$$ defined above is the same in almost all of these non-matching methods.

The embedded domain is assumed to be included in $$\Omega$$, which we take as the unit square $$[0,1]^2$$. The definition of the fictitious domain $$\Gamma$$ can be modified through the parameter file, and can be given as a mapping from the reference interval $$[0,1]$$ to a curve in $$\Omega$$.

If the curve is closed, then the results will be similar to running the same problem on a grid whose boundary is $$\Gamma$$. The program will happily run also with a non-closed $$\Gamma$$, although in those cases the mathematical formulation of the problem is more difficult, since $$\Gamma$$ will have a boundary by itself that has co-dimension two with respect to the domain $$\Omega$$.