# Implementation of mortar methods for coupling different meshes w/o overlap?

Asked by Christian Waluga on 2012-11-15

Dear all,

Is there a way to implement the coupling between two (possibly nonmatching) meshes via Nitsche-type methods or mortar methods based on Lagrange multipliers in FEniCS?

Moreover, is it possible to directly evaluate gradients of a finite element-solution at a given point (i.e., without solving additional problems or manually programming finite difference approximations)?

So far I only came up with quick hacks to implement such things, but I am curious, if there are natural ways to do this in FEniCS. My prototype basically implements (in an extremely sloppy and incredibly inefficient way) an iterative solver via a hybridized Nitsche-type mortar method as described in the following preprint:
http://www.aices.rwth-aachen.de:8080/aices/preprint/documents/aices-2009-2.pdf

Any comments on how to implement the coupling and the evaluation of the normal derivatives more elegantly are appreciated.

Best, Christian

Code:

from dolfin import *
import numpy

# domain and mesh parameters
lx = 1.0
ly = .5
nx = 30
ny = 15
scale = 3 # set to 1 for matching mesh

# polynomial degree
k = 1

# solver parameters
maxit = 500

# stabilization factor
alpha = 8.0*k**2/(lx/nx)

# subdomain meshes and interface mesh
mesh = [ Rectangle(0.0, 0.0, lx, ly, nx, ny), \
Rectangle(0.0, -ly, lx, 0.0, nx/scale, ny/scale) ]
meshi = Interval(nx, 0.0, lx)

mrange = range(0, len(mesh))

# normals
h = [ CellSize(mesh[i]) for i in mrange ]
n = [ FacetNormal(mesh[i]) for i in mrange ]

# function spaces
V = [ FunctionSpace(mesh[i], 'Lagrange', k) for i in mrange ]
M = FunctionSpace(meshi, 'DG', k)

# helpers for the coupling

class BCHelper(Expression):
def __init__(self, f):
self.f = f
def eval(self, value, x):
if near(x, 0.0):
self.f.eval(value, numpy.array([x]))
else:
value = 0.0 # add boundary conditions here
return

class FluxHelper(Expression):
def __init__(self, u):
self.u = u
def eval(self, value, x):
value = 0.0
eps = 1e-3
u0 = value*0.0 ; u0m = value*0.0; u1 = value*0.0; u1m = value*0.0
# evaluate solutions and approximate normal derivatives
self.u.eval(u0, numpy.array([x, 0.0]))
self.u.eval(u1, numpy.array([x, 0.0]))
self.u.eval(u0m, numpy.array([x, eps]))
self.u.eval(u1m, numpy.array([x, -eps]))
dnu0 = (u0-u0m)/eps; dnu1 = (u1-u1m)/eps
value = .5*(u0 + u1 - (dnu0 + dnu1)/alpha)
return

# Poisson bilinear form
def a(i, u, v):
+ alpha*dot(u, v)*ds

# right hand side
def l(i, g, v):
return 1.0*i*v*dx + alpha * dot(g, v)*ds - inner(g, dot(grad(v),n[i]))*ds

# variational problems
u = [ TrialFunction(V[i]) for i in mrange ]
v = [ TestFunction(V[i]) for i in mrange ]

# make discrete functions
uh = [ Function(V[i]) for i in mrange ]
lam = interpolate(Expression('0.0'), M)
err = 0.0; eps = 1e-5

# iteratively solve multi-domain problem
for j in range(0, maxit):

# save old interface value for comparison
lam0 = lam

# update interface flux
fh = FluxHelper(uh)
lam = interpolate(fh, M)
gi = BCHelper(lam)

# solve discrete problems
for i in mrange:
solve(a(i, u[i], v[i]) == l(i, gi, v[i]), uh[i], None)

# check correction
err = (lam.vector() - lam0.vector()).norm('l2')
print 'error({0}) = {1}'.format(j, err)
if j > 0 and err < eps: break

# visualization
class TwoDomainFunction ( Expression ):
def __init__(self, uh):
self.uh = uh
def eval(self, value, x):
if (x > 0.0):
self.uh.eval ( value, x )
if (x <= 0.0):
self.uh.eval ( value, x )
return

vismesh = Rectangle(0.0, -ly, lx, ly, nx, ny*2)
VV = FunctionSpace(vismesh, 'DG', k)
Uh = interpolate ( TwoDomainFunction ( uh ), VV )
plot ( Uh )
interactive ( )

## Question information

Language:
English Edit question
Status:
Expired
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Last query:
2012-11-15