Multimesh interface integrals are not computed correctly when more than two meshes overlap

[Magne Nordaas]

The following example demonstrates the problem. The computed perimeter, obtained with an interface integral, is incorrect when three meshes overlap.

from dolfin import *

N = 3
meshes = [UnitSquareMesh(2*N, 2*N),
          RectangleMesh(Point(0.1, 0.15), Point(0.6, 0.65), N, N),
          RectangleMesh(Point(0.4, 0.35), Point(0.9, 0.85), N, N)]

multimesh = MultiMesh()
for mesh in meshes:
    multimesh.add(mesh)
multimesh.build()

V = MultiMeshFunctionSpace(multimesh, "CG", 1)
u, v = TrialFunction(V), TestFunction(V)


a = inner(u,v) * dX
b = a + inner(jump(u),jump(v)) * dI
A = assemble_multimesh(a)
B = assemble_multimesh(b)

# create multimesh function w such that
#     w(x) =  1  if x is in the last mesh
#          =  0  otherwise  

w = MultiMeshFunction(V); x = w.vector()
dofs = V.part(V.num_parts()-1).dofmap().dofs() \
     + sum(V.part(i).dim() for i in range(V.num_parts()-1))
w.vector()[dofs] = 1.

# Compute the area and perimeter of the last mesh
a = w.vector().inner(A * w.vector())
p = w.vector().inner(B * w.vector()) - a

print("Computed area (a) and perimeter (p) of last mesh:")
print("  a = {0:1.4f} (exact value is .25)".format(a))
print("  p = {0:1.4f} (exact value is 2.0)".format(p))

It seems that the interface is not correctly identified. For example, the dofs marked blue (on the first mesh) and red (on the second mesh) on the figure below give a non-zero contribution to the interface integral. This should not be the case since the relevant part of the interface between the first two meshes is covered by the third mesh.

While writing up the example, I also noticed that there is an issue with the interface integral when cell edges coincide with mesh interface, even with just two meshes. This may be a separate issue.

[August Johansson]

I've looked into the blue node case, and there are indeed two cells sharing this node that are classified as covered that should have been classified as cut. The two elements are northwest and southeast of the blue node. This classification is done in MultiMesh::_build_collision_maps() if you want to have a look.

I'm also getting the wrong interface area, so there's clearly something wrong.

Coinciding interface integrals is difficult from a geometric robustness perspective and needs special care from an algorithmic perspective. We've spent a lot of time on the former and hopefully that will make the latter easier.

[Anders Logg (Chalmers)]

Added new component multimesh. We should use this tag for all issues related to multimesh.

[August Johansson]

I have partly resolved this issue: I now get the correct interface area using the generated quadrature weights. When I get the python interface working I can look into the specifics of the nodal values.