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:

(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.