Assembling over internal surfaces in a triangular domain

1

I am trying to integrate over internal surfaces, evaluating the total content of a function across vertical slices within the domain. I’m using dS (cf previous topics). I find in my real code, a zero value is (incorrectly) returned. But my solution function is Mixed and if I project the target component into a fresh Function, then it assembles into a nonzero value as expected. But it shouldn’t be necessary to have to project to get the job done.

My domain is a triangle. A minimum working example suggests this geometry is part of the problem. I’ve set up a constant Function with value π everywhere. Operating over a square unit mesh, integration over vertical slices returns 3.14159 as expected. But after changing to a triangular mesh, assembly over the VerticalSampleLines returns 0.

My MWE is

from dolfin import *
import mshr
import numpy

# define internal surface subdomains as vertical lines at position X
class VerticalSampleLine(SubDomain):
  def __init__(self, X, **kwargs):
    self.X = X
    super().__init__(**kwargs)

  def inside(self, x, on_boundary):
    return near(x[0],self.X)


# simple 2D mesh
useSquareMesh = False
if useSquareMesh:
  mesh = UnitSquareMesh(10,10)
else:
  unitRightTriangle = mshr.Polygon([Point([0,0]),Point([1,0]),Point([1,1])])
  mesh = mshr.generate_mesh(unitRightTriangle, 10)
V = FunctionSpace(mesh, "CG", 1)

# use a constant function with value π
f = project(Constant(numpy.pi),V)

verticalBoundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# set up a set of 5 internal surfaces to test
xset = numpy.linspace(0.1,0.9, 5)
verticalBoundarySet=[]
for i, x in enumerate(xset):
    lineIndex=i+1
    verticalBoundarySet.append( VerticalSampleLine(x) )
    verticalBoundarySet[lineIndex-1].mark(verticalBoundaries, lineIndex)

dS = Measure('dS', domain=mesh, subdomain_data=verticalBoundaries)

for i, verticalBoundary in enumerate(verticalBoundarySet):
    lineIndex=i+1
    lineValue = assemble(f*dS(lineIndex))
    
    # vertical line integrals should = π
    print("lineValue {} at x={:.1} is {}".format(lineIndex,verticalBoundary.X,lineValue))

With the triangular mesh (useSquareMesh = False) it prints

lineValue 1 at x=0.1 is 0.0
lineValue 2 at x=0.3 is 0.0
lineValue 3 at x=0.5 is 0.0
lineValue 4 at x=0.7 is 0.0
lineValue 5 at x=0.9 is 0.0

Set useSquareMesh = True to see the behaviour on the square mesh:

lineValue 1 at x=0.1 is 3.1415926535897936
lineValue 2 at x=0.3 is 3.141592653589793
lineValue 3 at x=0.5 is 3.141592653589794
lineValue 4 at x=0.7 is 3.141592653589794
lineValue 5 at x=0.9 is 3.1415926535897944

Why does lineValue get set to 0 when the triangular mesh is used (useSquareMesh = False) ?

2

Hi, on a triangular mesh there seem to be no interior facets that get marked by VerticalSampleLine, i.e. dS(i) is a zero measure for every i. Consider

# Will fail for triangular mesh
for lineIndex, x in enumerate(xset, 1):
    verticalBoundarySet.append( VerticalSampleLine(x) )
    verticalBoundarySet[lineIndex-1].mark(verticalBoundaries, lineIndex)

    tagged, = np.where(verticalBoundaries.array() == lineIndex)
    assert len(tagged)
3

Thanks @MiroK. That is deeply disturbing if the internal surface is arbitrarily not integrable like that. Is it a problem in principle with any triangular mesh or does it just mean mshr is buggy? I could use dmsh instead if the problem is mshr.

4

dS integrates over internal facets of the mesh. If the intersection of the surface where you wish to integrate with the facets is empty there is no wonder that the integration is zero. So this is not bug of mshr or gmsh but a consequence of the fact that for no facets in the mshr generated mesh VerticalSampleLine.inside returns True. You can either make the mesh conform to the line or consider the following

# Get the surface as independent mesh
surface = BoundaryMesh(UnitSquareMesh(10, 10), 'exterior')

# The curve is contained in some larger mesh
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 20, 20)
V = FunctionSpace(mesh, 'CG', 1)
f = interpolate(Expression('x[0]+x[1]', degree=1), V)

TV = FunctionSpace(surface, 'CG', 1)
Tf = interpolate(f, TV)

print(assemble(inner(Tf, Tf)*dx))
5

I think I see what you mean. It means assembly over a surface is performed strictly over the facets that cross that surface.

What I was expecting (or hoping) was that interpolation would be done over the intersection of the surface with the internal volume of the cells. We can evaluate a function pointwise at any given Point in the volume, so I was expecting integration over an intersecting surface which generalises on that pointwise operation.

I wonder if any such intersecting surface integration is interesting as a topic for study in finite elements. On the other hand, perhaps the computational cost of determining the intersection of surfaces is sufficiently high that there’s no advantage over your suggestion of just creating a second FunctionSpace to interpolate onto.