I wrote a small example that was originally meant to ask a few questions here on the forum, but while doing so I figured the answers out myself. I will post it anyway as a reference for others. And perhaps to hear if there is any easier way than the workarounds I use
Included in the example are:
- How to define a SubMesh from a CellFunction?
- How to redefine MeshFunctions from the original mesh to the submesh?
- How to rewrite a CellFunction to a FacetFunction (to use in a DirichletBC)?
Hope it is helpful to anyone
On a related note: when will MeshView be released?
from dolfin import *
# define regions and boundaries
mesh = UnitSquareMesh(8, 16)
regions = CellFunction('size_t', mesh)
boundaries = FacetFunction('size_t', mesh)
class A(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5 and x[1] <= 0.75
a = A()
class B(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= 0.5 and x[1] <= 0.75
b = B()
class C(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.75
c = C()
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
l = LeftBoundary()
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
r = RightBoundary()
regions.set_all(0)
b.mark(regions, 1)
c.mark(regions, 2)
boundaries.set_all(0)
l.mark(boundaries, 1)
r.mark(boundaries, 2)
plot(regions, title="Regions")
plot(boundaries, title="Boundaries")
# define a submesh, composed of region 0 and 1
new_region = CellFunction('size_t', mesh)
new_region.set_all(0)
new_region.array()[regions.array() == 0] = 1
new_region.array()[regions.array() == 1] = 1
submesh = SubMesh(mesh, new_region, 1)
# define new meshfunctions on this submesh with the same values as their original mesh
submesh_regions = CellFunction('size_t', submesh)
submesh_regions.set_all(0)
submesh_boundaries = FacetFunction('size_t', submesh)
submesh_boundaries.set_all(0)
vmap = submesh.data().array('parent_vertex_indices', 0)
cmap = submesh.data().array('parent_cell_indices', 2)
for submesh_cell in cells(submesh):
parent_cell = cmap[submesh_cell.index()]
submesh_regions.array()[submesh_cell.index()] = regions.array()[parent_cell]
for submesh_facet in facets(submesh):
submesh_facet_vertices = vmap[submesh_facet.entities(0)]
for facet in facets(mesh):
if set(facet.entities(0)) == set(submesh_facet_vertices):
submesh_boundaries.array()[submesh_facet.index()] = boundaries.array()[facet.index()]
break
plot(submesh_regions, title="Submesh regions")
plot(submesh_boundaries, title="Submesh boundaries")
# solve some equation on the whole domain
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
a = inner(grad(u), grad(v)) * dx
L = Expression('x[0]*x[0] + 10*x[0]*x[1]*x[1]') * v * dx
bcs = [DirichletBC(V, Constant(0.0), boundaries, 1)]
solve(a == L, f, bcs)
plot(f, title="Function f on the whole mesh")
# use the result in a problem defined on the subdomain
Vsub = FunctionSpace(submesh, 'CG', 1)
u = TrialFunction(Vsub)
v = TestFunction(Vsub)
g = Function(Vsub)
dx = Measure('dx')[submesh_regions]
a = inner(grad(u), grad(v)) * dx
class F(Expression):
def eval(self, value, x):
value[0] = f(x)
L = F() * v * dx
bcs = [DirichletBC(Vsub, Constant(0.0), boundaries, 1)]
solve(a == L, g, bcs)
plot(g, title="Function g on the submesh")
# rewrite the cellfunction to a facetfunction
new_region_facets = FacetFunction('size_t', mesh)
new_region_facets.set_all(0)
mesh.init(1, 2)
for fac in facets(mesh):
if 1 in new_region.array()[fac.entities(2)]:
new_region_facets[fac.index()] = 1
class G(Expression):
def eval(self, value, x):
value[0] = g(x)
plot(new_region_facets, title="Facetfunction used to define a DirichletBC")
bcs = [DirichletBC(V, G(), new_region_facets, 1)]
# extrapolate the resulting function back to the whole domain
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
h = Function(V)
dx = Measure('dx')[new_region]
a = inner(grad(u), grad(v)) * dx
L = Constant(0.0) * v * dx
solve(a == L, h, bcs)
plot(h, title="Function g, extrapolated to the whole mesh")
interactive()
