Coupling to the interface

This is an example of how to use Dune-MMesh and solve coupled problems on the bulk and interface grid.

Grid creation

We use the horizontal grid file that contains an interface \(\Gamma = [0.25, 0.75] \times {0.5}\) embedded in a domain \(\Omega = [0,1]^2\). Grid creation from a mesh file works as follow.

[1]:
from dune.grid import reader
from dune.mmesh import mmesh
dim = 2
file = "grids/horizontal.msh"
gridView  = mmesh((reader.gmsh, file), dim)
igridView = gridView.hierarchicalGrid.interfaceGrid

Solve a problem on the bulk grid

Let us solve the Poisson equation

\begin{align} -\Delta u = f & \qquad \text{in}\ \Omega \end{align}

on the bulk grid. We use the manufactured solution \(\hat u = \sin(4 \pi x y)\) and therefore apply the source term \(f = -\operatorname{div}( \nabla \hat u )\). The weak form of the problem above reads \begin{align} \int_\Omega \nabla u \cdot \nabla v~dx &= \int_\Omega f v~dx \end{align} for all corresponding test functions \(v\). This can be implemented as follows.

[2]:
from ufl import *
from dune.ufl import DirichletBC
from dune.fem.space import lagrange
from dune.fem.scheme import galerkin
from dune.fem.function import integrate

space = lagrange(gridView, order=3)
u = TrialFunction(space)
v = TestFunction(space)

x = SpatialCoordinate(space)
exact = sin(x[0]*x[1]*4*pi)
f = -div(grad(exact))

a = inner(grad(u), grad(v)) * dx
b = f * v * dx

scheme = galerkin([a == b, DirichletBC(space, exact)], solver=("suitesparse", "umfpack"))
uh = space.interpolate(0, name="solution")
scheme.solve(target=uh)

def L2(u1, u2):
    return sqrt(integrate(u1.grid, (u1-u2)**2, order=5))

L2(uh, exact)
[2]:
5.628259763933402e-07

Solve a problem on the interface

We can solve similar problem on the interface \(\Gamma\) like \begin{align} -\Delta u_\Gamma = f & \qquad \text{in}\ \Gamma \end{align} with the weak form \begin{align} \int_\Gamma \nabla u_\Gamma \cdot \nabla v_\Gamma~dx &= \int_\Gamma f v_\Gamma~dx \end{align} for all corresponding test functions \(v_\Gamma\).

[3]:
ispace = lagrange(igridView, order=3)
iuh = ispace.interpolate(0, name="isolution")

iu = TrialFunction(ispace)
iv = TestFunction(ispace)

ix = SpatialCoordinate(ispace)
iexact = sin(0.5*ix[dim-2]*4*pi)
iF = -div(grad(iexact))

ia = inner(grad(iu), grad(iv)) * dx
ib = iF * iv * dx

ischeme = galerkin([ia == ib, DirichletBC(ispace, iexact)])
ischeme.solve(target=iuh)
L2(iuh, iexact)
[3]:
5.807742030532398e-08

We can use the plotPointData function to visualize the solution of both grids.

[4]:
import matplotlib.pyplot as plt
from dune.fem.plotting import plotPointData as plot

figure = plt.figure(figsize=(3,3))
plot(uh, figure=figure, gridLines=None)
plot(iuh, figure=figure, linewidth=0.04, colorbar=None)
../_images/examples_coupling_9_0.png

Couple bulk to surface

Dune-MMesh makes it possible to compute traces of discrete functions on \(\Omega\) along \(\Gamma\).

[5]:
from dune.mmesh import trace
tr = avg(trace(uh))
ib = inner(grad(tr), grad(iv)) * dx

iuh.interpolate(0)
ischeme = galerkin([ia == ib, DirichletBC(ispace, avg(trace(uh)))])
ischeme.solve(target=iuh)
L2(iuh, iexact)
[5]:
4.266679479547976e-08

Couple surface to bulk

Similarly, we can evaluate a discrete function on \(\Gamma\) at the skeleton of the triangulation of \(\Omega\).

[6]:
from dune.mmesh import skeleton

sk = skeleton(iuh)
b = avg(sk) * avg(v) * dS

uh.interpolate(0)
scheme = galerkin([a == b, DirichletBC(space, 0)])
scheme.solve(target=uh)

figure = plt.figure(figsize=(3,3))
plot(uh, figure=figure, gridLines=None)
../_images/examples_coupling_13_0.png

Compute jump of gradient of traces

We can use trace and skeleton within UFL expressions. This is useful when using jump and average operators, but also to compute gradients and more. Together with the utility function normals, which returns the positive normal of the underlying bulk facet, we can determine the orientation of the restricted values.

[7]:
from dune.mmesh import normals
inormal = normals(igridView)
jmp = jump(grad(trace(uh)), inormal)