FEnics mesh generation - Mark inner region

Question

I am solving the classical EM problem; scattering from a conducting sphere, so far i 2D. Assuming the sphere to be a perfect conductor, the problem can be solved on a mesh with a (circular) hole representing the sphere (along with appropriate BCs). This mesh is readily generated in FEniCS.

If a particular material (e.g. silver) is to be considered, the inner region of the hole must be meshed too. The "easy way out" would be to use a simple, rectangular mesh . However, a such mesh does not respect the boundary between the two regions, so the circular shape obtained is very poor.

I can construct a mesh externally, that enforces the boundary, but i was wondering if it is possible to construct the mesh from within FEniCS?

I know how to construct the circular "inner" domain as well as the "outer" domain (as demonstrated), so it seems somehow that the functionality is there.

Sample code for generating figures:

from dolfin import *
from dolfin.cpp.mesh import *
from mshr import *

class LambdaDomain(SubDomain):
def __init__(self, condition):
self.condition = condition
SubDomain.__init__(self)

def inside(self, x, on_boundary):
return self.condition(x, on_boundary)

lx = ly = 10
cx = cy = lx / 2.0
radius = 2.5
res = 25

# Generate mesh with hole.
base = Rectangle(Point(0, 0), Point(lx, ly))
hole = Circle(Point(cx, cy), radius)
plot(generate_mesh(base - hole, res))

# Mark hole in full mesh using mesh function.
mesh = generate_mesh(base, res)
mesh_func = MeshFunction("size_t", mesh, 2)
mesh_func.set_all(0)
LambdaDomain(lambda x, on: ((x[0] - cx) ** 2 + (x[1] - cy) ** 2) ** 0.5 <= radi    us).mark(mesh_func, 1)
plot(mesh_func)

interactive()

Answer 1

Yes. it's possible. Consider:

from dolfin import *
from mshr import *

set_log_level(TRACE)

cx, cy, radius = 0.5, 0.5, 0.25
lx, ly = 1.0, 1.0

class circle(SubDomain):
    def inside(self, x, on_boundary):
        return pow(x[0] - cx, 2) + pow(x[1] - cy, 2) <= pow(radius, 2)

# Define 2D geometry
domain = Rectangle(Point(0.0, 0.0), Point(lx, ly))
domain.set_subdomain(1, Circle(Point(cx, cy), radius))


# Generate and plot mesh
mesh2d = generate_mesh(domain, 10)
plot(mesh2d, "2D mesh")

# Convert subdomains to mesh function for plotting
mf = MeshFunction("size_t", mesh2d, 2)
mf.set_all(0)
circle = circle()

for c in cells(mesh2d):
  if circle.inside(c.midpoint(), True):
    mf[c.index()] = 1


plot(mf, "Subdomains")
interactive()

Comments

Use

mf = dolfin.MeshFunction("size_t", mesh2d, 2, mesh2d.domains())

instead of

mf = MeshFunction("size_t", mesh2d, 2)
mf.set_all(0)
circle = circle()

for c in cells(mesh2d):
  if circle.inside(c.midpoint(), True):
    mf[c.index()] = 1

also should work (reference: csg-subdomains-2D.py).

---

Thanks! The set_subdomain function was exactly what i was looking for. However, using it encounter a weird problem. Conder the MWE:

from dolfin import *
from mshr import *

lx = ly = 4.5
cx = cy = lx/2.0
radius = lx/4.0

# Define 2D geometry
domain = Rectangle(Point(0.0, 0.0), Point(lx, ly))
domain.set_subdomain(1, Circle(Point(cx, cy), radius))
# Generate and plot mesh
mesh2d = generate_mesh(domain, 10)
plot(mesh2d, "2D mesh")

interactive()

If i change the value of lx to more than 4.5, e.g. 10 as in my original problem, no mesh is generated. The problem occurs only when i use the set_subdomain function. It is not a major problem; i could just scale the dimensions of my problem, but it seems weird. Maybe a bug?

---

Yes, see Issue #33 for more information. Install the development version of the mshr library seems to be the solution.

---

I can confirm that switching to FEniCS 1.7.0dev fixed the problem.