 |
Dune-Fufem 2.11-git
|
Table of Contents
- Introduction
- Mathematical background
- Setup and initialization
- Problem data
- Grid creation
- Setup of subdomains and boundary patches
- Selection of finite element basis
- Matrix, vector, and constraints data structures
- Global assembly of the problem
- Algebraic solution
- Writing the solution to VTK
- Error handling
Introduction
This example shows how to couple bulk and surface PDEs by presenting the coupling of elliptic problems in two bulk subdomains to an elliptic problem on the surface given by the interface between the subdomains.
Mathematical background
For a domain \(\Omega \subset \mathbb{R}^d\) let \(\Omega_A\) be an open subset with \(\overline{\Omega_A} \subset \Omega\) and set \(\Omega_B = \Omega \setminus \overline{\Omega_A}\). In the following we assume that the interface \(\Gamma = \overline{\Omega_A} \cap \overline{\Omega_B}\) is sufficiently smooth and consider bulk concentrations \(u_A : \Omega_A \to \mathbb{R}\) and \(u_B : \Omega_B \to \mathbb{R}\) and a surface concentration \(u_\Gamma : \Gamma \to \mathbb{R}\). In the bulk we impose the Poisson-problems
\[ \begin{align*} -\Delta u_A &= f_A && \text{ in } \Omega_A, \\ -\Delta u_B &= f_B && \text{ in } \Omega_B. \end{align*} \]
We assume that mass transfer to the interface is proportional to the concentration difference
\[ \begin{align*} -\frac{\partial u_A}{\partial n_A} &= \alpha (u_A - u_\Gamma) && \text{ on } \Gamma,\\ -\frac{\partial u_B}{\partial n_B} &= \alpha (u_B - u_\Gamma) && \text{ on } \Gamma, \end{align*} \]
where \(\alpha>0\) is a constant and \(n_A\) and \( n_B\) are the unit outward normals to \(\partial\Omega_A\) and \(\partial\Omega_B\), respectively. On the surface \(\Gamma\) we consider tangential diffusion for \(u_\Gamma\) with the outflow of the two bulk domains as source terms
\[ \begin{align*} -\delta\Delta_\Gamma u_\Gamma + \frac{\partial u_A}{\partial n_A} + \frac{\partial u_B}{\partial n_B} &= f_\Gamma && \text{ on } \Gamma. \end{align*} \]
On the outer boundary \(\partial \Omega = \partial \Omega_B\subset \Gamma\) we impose a Dirichlet boundary condition \(u_B=0\) as encoded in the space \( H_{\partial \Omega}^1(\Omega_B) := \{v \in H^1(\Omega_B)\,|\, v|_{\partial \Omega}=0\} \). Testing the equations with, \(v_A, v_B\), and \(v_\Gamma\), integrating over the respective domains, using partial integration, and inserting the boundary conditions we arrive at the following variational problem:
Find \( (u_A, u_B, u_\Gamma) \in \mathcal{H} := H^1(\Omega_A) \times H_{\partial \Omega}^1(\Omega_B) \times H^1(\Gamma)\) such that
\[ \small \begin{align*} \int_{\Omega_A}\nabla u_A \cdot \nabla v_A dx + \alpha\int_\Gamma (u_A - u_\Gamma) v_A d\sigma &= \int_{\Omega_A} f_A v_A dx && \forall v_A \in H_{\partial \Omega}^1(\Omega_A),\\ \int_{\Omega_B}\nabla u_B \cdot \nabla v_B dx + \alpha\int_\Gamma (u_B - u_\Gamma) v_B d\sigma &= \int_{\Omega_B} f_B v_B dx && \forall v_B \in H^1(\Omega_B),\\ \int_{\Gamma}\delta \nabla u_\Gamma \cdot \nabla v_\Gamma - \alpha(u_A - u_\Gamma) v_\Gamma - \alpha(u_B - u_\Gamma) v_\Gamma d\sigma &= \int_{\Gamma} f_\Gamma v_\Gamma d\sigma && \forall v_\Gamma \in H^1(\Gamma). \end{align*} \]
With the bilinear form
\[ \begin{multline} a(u,v) = \int_{\Omega_A}\nabla u_A \cdot \nabla v_A dx + \int_{\Omega_B}\nabla u_B \cdot \nabla v_B dx + \delta\int_{\Gamma} \nabla u_\Gamma \cdot \nabla v_\Gamma d\sigma \\ + \alpha\int_\Gamma (u_A - u_\Gamma) (v_A - v_\Gamma) + (u_B - u_\Gamma) (v_B - v_\Gamma) d\sigma \end{multline} \]
and linear form
\[ \begin{align*} \ell(v) &= \int_{\Omega_A}f_A v_B dx + \int_{\Omega_B}f_B v_B dx + \int_{\Gamma}f_\Gamma v_\Gamma d\sigma \end{align*} \]
for tuples \(u=(u_A,u_B,u_\Gamma)\) and \(v=(v_A,v_B,v_\Gamma)\) this can be compactly written as
\[ \begin{align*} u &\in \mathcal{H}: & a(u,v) &= \ell(v) && \forall v \in \mathcal{H}. \end{align*} \]
It is clear that the bilinear form is symmetric and continuous and that it satisfies
\[ a(v,v) \geq \|\nabla v_A\|_{L^2(\Omega_A)}^2 + \|\nabla v_B\|_{L^2(\Omega_B)}^2 + \delta \|\nabla_\Gamma v_\Gamma\|_{L^2(\Gamma)}^2. \]
Hence, by the classic Poincaré-inequality it is coercive on the orthogonal complement of the (component-wise) constant functions. On the other hand all functions in the kernel of \(a(\cdot,\cdot)\) must satisfy \(v_A=v_\Gamma=v_B\) on \(\Gamma\) and \(v_B=0\) on \(\partial \Omega\). Thus the kernel is trivial and the bilinear form is coercive (see, e.g., [9], Proposition 1) and the Lax-Milgram theorem provides well-posedness.
The discretization error of the finite element discretization use in the example has been investigated (for a slightly simpler problem with a single subdomain) in [6].
Setup and initialization
After including the required headers and defining a utility function that will be explained below, the main() function first creates a simple logger that allows to write log messages to std::cout (see the Poisson example for details).
Problem data
Next we parse the command line arguments into a Dune::ParameterTree and extract the parameters for the mesh file, the number of grid refinements, the tolerance, maximal iteration count and verbosity of the iterative solver, and the problem parameters \(\alpha\) and \(\delta\). Furthermore we define the functions for the right-hand side of the equations in the bulk domains (A) and (B) and on the interface (I), and the Dirichlet values on the boundary.
Grid creation
The grid is created by reading a Gmsh file into a GridFactory and extracting the element data that is used to define the two subdomains. Afterward the grid is refined and the element data is updated to follow the refinement. These steps as well as the handling of subdomains and interfaces and the restriction of finite element bases to (bulk) subdomains are described in more detail in the Stokes-Darcy example.
Setup of subdomains and boundary patches
The element data is used to setup two SubDomain objects representing the respective subdomains \(\Omega_A\) and \(\Omega_B\). For the surface equation we create an interface object representing the interface \(\Gamma\) between the subdomains. Notice that the interface is oriented from subDomain_A to subDomain_B which will be important later on.
Selection of finite element basis
Next we create the global finite element basis using first order finite elements in both bulk subdomains and on the interface.
For the interface we want to use the so called trace finite element approach which uses the trace of a bulk finite element space. To this end we create a basis that lives on an bulk domain, but will only assemble terms on the interface later on. This would make the assembled systems singular. Hence we will constrain the coefficients for all degrees of freedom that do not live on the interface to zero, similar to Dirichlet constraints. To reduce the number of DOFs to constrain, we start with a bulk basis defined on subDomain_A which contains all the DOFs we need for the trace space.
- Note
- Ideally one would be able to write
restrict(lagrange<1>(), interface)to restricting the bulk basis to the interface and only get the required DOFs. However, the dune-functions module does not support the restriction to lower dimensional subsets yet. Hence we have to use the 'trick' described above.
Matrix, vector, and constraints data structures
Since the basis' indices are plain integers we next define types for flat vectors and matrices of reals numbers and create a solution and a right-hand-side vector and a stiffness matrix.
To handle the Dirichlet and trace basis constraints we create a constraints object and compute the respective constraints.
While the `computeBoundaryConstraints()` function is provided by dune-fufem, computeTraceConstraints() is defined in the present example above the main() function as follows:
The function first marks all DOFs associated to intersections from the interface using Dune::Fufem::markIntersectionDofs. The vector for storing the respective flags is obtained by copying the one from the constraints. To access the vector with multi-indices of the basis (that may in general have multiple digits) it is wrapped into VectorBacked. Next the functions loops over all degrees of freedom of the respective (bulk) basis and marks the non-interface DOFs as constrained and sets the respective constraint value to zero.
Global assembly of the problem
To assemble the matrix and right hand side vector we use the Dune::Fufem::Forms framework to define the respective variational forms. For an explanation how to define such forms within this framework and how to assemble the associated matrices and vectors using the dune-assembler module refer to the Poisson example. Advanced features like the computation of subdomain and interface terms are explained in the Stokes-Darcy example.
Given the explanations provided in the other examples, the definition of the variational forms used here is essentially straight-forward. We only note that for the interface term we can use the expressions u_A, v_A, u_I, and v_I directly, while u_B and v_B have to be wrapped in outside(...). This is a consequence of the fact that the interface is oriented from subDomain_A (called inside) to subDomain_B (called outside). All terms that should be evaluated on the inside of the interface can be used directly, while one has to explicitly declares if terms should be evaluated on the outside.
The bilinear and linear form defined above can now be passed to the methods of the global Assembler class provided by the dune-assembler module to assemble the associated matrix and vector.
Once the system is assembled, we can constrain it to the affine subspace satisfying the constraints that we defined earlier. This includes the Dirichlet boundary constraints and the constraints to the trace for the interface basis.
Algebraic solution
Since the system is symmetric positive definite we use a preconditioned CG-method with an AMG preconditioner provided by the dune-istl module for the iterative solution.
Writing the solution to VTK
Finally the solution is written to a VTK file understood by the ParaView visualization program. We write the three components of the bulk and interface solutions independently and for convenience additionally the sum of both bulk solutions. All the written functions are obtained by binding the respective Dune::Fufem::Forms expressions to the solution vector using the ...|sol operator. Since the solutions are discontinuous across the interface, we use the DiscontinuousLagrangeDataCollector which writes the solution independently for each grid element.
Error handling
In case an exception is thrown within main(), the error message is printed and the program is terminated.
-
Legal Statements / Impressum
|
Hosted by TU Dresden & Uni Heidelberg |
Generated by
1.9.8