This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

SubMesh, MeshFunction, DirichletBC, MeshView etc.

+8 votes

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()
asked Mar 4, 2016 by maartent FEniCS User (3,910 points)
edited Mar 4, 2016 by maartent

1 Answer

0 votes

Unfortunately does not work with 2016.1.0

Error: Unable to create Dirichlet boundary condition.
Reason: User MeshFunction and FunctionSpace meshes are different.
Where: This error was encountered inside DirichletBC.cpp.

answered Aug 24, 2016 by meigel FEniCS User (1,520 points)

I am limited to version 1.4 (and occasionally 1.6) so I didn't notice yet. But good to know if I ever decide to update to the latest version... From what I read here and there, there is a less hackish solution on the roadmap for this kind of stuff, so I pin my hopes on that.

Any news on this new solution? I just ran into the same problem.

...