How to find boundary nodes on the inner boundary of a doubly connected region?

1

Consider a doubly connected region, for example, an annulus with inner radius 1 and outer radius 2. After meshing the annulus, the coordinates of the nodes on the exterior and the interior boundary (collectively) can be found by

BoundaryMesh(mesh, "exterior", True)

I need to find the coordinates of the inner and the outer boundary nodes separately. I can do that easily for this simple geometry. But in case of a doubly connected region with general smooth analytic boundary curves, how to do it? I tried BoundaryMesh(mesh, "interior", True), but that didn’t work and returned a null set.

Consider the simple MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from dolfin import *
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C_1 = C2 - C1


class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class b2(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4
        
mesh = generate_mesh(C_1, 5)

boundary_nodes = BoundaryMesh(mesh, "exterior",  True)
boundary_coordinates = boundary_nodes.coordinates()
n1 = len(boundary_coordinates)
print('number of nodes on the boundary =',n1)
print('outer nodes coordinates =', boundary_coordinates)

boundary_nodes_inner = BoundaryMesh(mesh, "interior",  True)
boundary_coordinates_inner = boundary_nodes_inner.coordinates()
n2 = len(boundary_coordinates_inner)
print('number of nodes on the boundary =',n2)
print('inner nodes coordinates =', boundary_coordinates_inner)