Define boundaries for intersection of two cylinders

variational formulation
2

Please consider the following:

class cyl1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and np.isclose(np.sqrt((x[0]-c1)**2+(x[1]-c1)**2), 0.075, rtol=3e-2)

class cyl2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and np.isclose(np.sqrt((x[0]-c2)**2+(x[1]-c1)**2), 0.075, rtol=3e-2)


mf = MeshFunction("size_t", mesh, 2, 0)
cyl1().mark(mf, 2)
cyl2().mark(mf, 3)
File("mf.pvd") << mf

Then you can use the DirichletBC interface to enforce the boundary condition