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

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.

Thanks for your time!

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[1], 0.0):
      self.f.eval(value, numpy.array([x[0]]))
    else:
      value[0] = 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.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[0].eval(u0, numpy.array([x[0], 0.0]))
    self.u[1].eval(u1, numpy.array([x[0], 0.0]))
    self.u[0].eval(u0m, numpy.array([x[0], eps]))
    self.u[1].eval(u1m, numpy.array([x[0], -eps]))
    dnu0 = (u0-u0m)/eps; dnu1 = (u1-u1m)/eps
    value[0] = .5*(u0 + u1 - (dnu0 + dnu1)/alpha)
    return

# Poisson bilinear form
def a(i, u, v):
  return inner(grad(u),grad(v))*dx \
    - inner(dot(grad(u),n[i]), v)*ds \
    - inner(u, dot(grad(v),n[i]))*ds \
    + 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[1] > 0.0):
      self.uh[0].eval ( value, x )
    if (x[1] <= 0.0):
      self.uh[1].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 ( )