Obtain connected mesh of boundary

I am implementing a boundary condition that is based on precomputing quantities
based on geometric information on the boundary (2D domain, 1D boundary).
I want to implement this by extracting a submesh of the boundary, which I can then do computations on. I have tried this post: Extract parts of the boundary

which is written for 3D. Translated to 2D, my code looks like this:

import dolfin as dl
import mshr as ms

domain = ms.Rectangle(dl.Point(0,0),dl.Point(10,10))
sphere = ms.Circle(dl.Point(5, 5), 3)
domain = domain - sphere
mesh   = ms.generate_mesh(domain,10)
exterior = dl.BoundaryMesh(mesh, 'exterior')

from numpy import array
from numpy.linalg import norm
def inSphere(x):
    v = x - array([5,5])
    return norm(v)<3+1e2*dl.DOLFIN_EPS
class SphereDomain(dl.SubDomain):
    def inside(self,x,on_boundary):
        return inSphere(x)
class BoxDomain(dl.SubDomain):
    def inside(self,x,on_boundary):
        return not inSphere(x)
    
print(dl.__version__)
dl.plot(dl.SubMesh(exterior, SphereDomain()))

With the output:

bild

(This is dolfin v 2019.1.0)

But as you see, the resulting mesh is not connected in the way you’d expect. Any help on this would be appreciated. How do I make the resulting submesh agree with the boundary?

Edit:

Thanks to user dokken, I now know that the mesh is connected, it’s the plotting that misrepresents the mesh structure. The following manual plot:

bmesh = dl.SubMesh(exterior, SphereDomain())
for f in dl.cells(bmesh):
    pts = [v.point() for v in dl.vertices(f)]
    plt.plot([pts[0].x(), pts[1].x()], [pts[0].y(), pts[1].y()], 'r.-')

Shows a correct visualisation of the boundary mesh.

This is more likely an artifact of the plotting with matplotlib, rather than a bug with the mesh connectivity. I would suggest writing the submesh to xdmf and open it with paraview for proper inspection.

You could Also assemble print(assemble(1*dx(domain=submesh_variable_name)))
And check what value you get

1 Like

Thanks a bunch! I tried plotting edge by edge and sure enough, the mesh is actually connected.