Can I use CSG primitive to set boundary condition?

Question

Like the original poster, I have a 2D domain with holes that I can describe using the CSG primitives Rectangle and Polygon. Since I already have the boundaries as primitives, it would seem natural to use these for defining boundary conditions. My impression is that the way to do this is: first, create a MeshFunction that identifies each boundary, then use the function in the boundary conditions.

I'm having no success in even constructing the function. Based on code provided in previous answer [https://github.com/FEniCS/mshr/blob/master/demo/python/csg-subdomains-2D.py] I tried constructing a mesh function after setting a subdomain using the inner boundary.

from dolfin import *
from mshr import *

outer = Rectangle(Point(0,0), Point(10,10))
inner = Circle(Point(5,5), 2)
domain = outer - inner
domain.set_subdomain(42, inner)

mesh = generate_mesh(domain, 5)
plot(mesh)

mf = MeshFunction('size_t', mesh, 2, mesh.domains())
plot(mf)
File("mf.xml") << mf

interactive()

However, the mesh function is identically 0 everywhere. In this simple example, I can just enlarge the inner radius a bit and create a thin shell where the function is non-zero. But this is not going to work as easily in my real case where the inner boundary is a polygon.

Q1. How can I achieve my goal?

Q2. For a boundary condition, should I be using a 1D mesh function instead? I tried changing the third argument of MeshFunction() from "2" to "1", but it was no better. To my surprise, the output function seems to still be defined over triangles -- the cell index number ranges over 0-459, matching the number of mesh triangles . The value is 2^64-1.

Answer 1

In the end, I found no way to do this within fenics. But I was able to write some Python code that does what I need.

First, I should note that -- despite using "Circle" in my question -- all my input boundaries are polygons. After meshing and refinement, I know that the boundary edges of the final mesh are simply subdivided segments of the input boundary edges. So I coded up an edge walking function: given two vertices incident on an input edge, find all the vertices of the final mesh that lie on that line segment. From this, I get a list of mesh edges for each boundary.

To create the boundary conditions, I use an EdgeFunction that is assigned a different value for each of the boundaries. This function is set up as my edge walker discovers edges belonging to a given input boundary polygon.