Vertex mapping from a SubMesh of a BoundaryMesh back to the actual Mesh!

Question

So we can easily get the vertex mapping from a SubMesh back to a BoundaryMesh, as described here, utilizing

t = submesh.data().array('parent_vertex_indices', 0)

However, if submesh is a SubMesh of a BoundaryMesh, the BoundaryMesh does not have any 'parent_vertex_indices'. My question is, how do i map back the the original mesh:

from dolfin import *
from pylab  import *

mesh = UnitCubeMesh(10,10,10)          # original mesh
mesh.coordinates()[:,0] -= .5          # shift x-coords
mesh.coordinates()[:,1] -= .5          # shift y-coords
V    = FunctionSpace(mesh, "CG", 1)
u    = Function(V)

# apply expression over cube for clearer results :
u_i  = Expression('sqrt(pow(x[0],2) + pow(x[1], 2))')
u.interpolate(u_i)

bmesh  = BoundaryMesh(mesh, "exterior")   # surface boundary mesh

# mark the boundary of the bottom surface :
cellmap = bmesh.entity_map(2)
vertmap = bmesh.entity_map(0)
pb      = CellFunction("size_t", bmesh, 0)
for c in cells(bmesh):
  if Facet(mesh, cellmap[c.index()]).normal().z() < 0:
    pb[c] = 1

submesh = SubMesh(bmesh, pb, 1)           # bottom of boundary mesh

Vb  = FunctionSpace(bmesh,   "CG", 1)     # surface function space
Vs  = FunctionSpace(submesh, "CG", 1)     # submesh function space

ub  = Function(Vb)                        # boundary function
us  = Function(Vs)                        # surface function

ub.interpolate(u)                         # interpolate u onto boundary
us.interpolate(u)                         # interpolate u onto surface mesh

unb = Function(Vb)                        # new boundary function
un  = Function(V)                         # new whole function


# mappings we may need :
m    = vertex_to_dof_map(V)
b    = vertex_to_dof_map(Vb)
s    = vertex_to_dof_map(Vs)

mi   = dof_to_vertex_map(V)
bi   = dof_to_vertex_map(Vb)
si   = dof_to_vertex_map(Vs)

# mapping from submesh back to bmesh :
t = submesh.data().array('parent_vertex_indices', 0)

# get vertex-valued arrays :
us_a  = us.vector().array()
u_a   = u.vector().array()
ub_a  = ub.vector().array()

unb_a = unb.vector().array()
un_a  = un.vector().array()

# update the values of the new functions to be the values of the surface :
unb_a[b[t]]  = us_a[s]   # works

un_a[m[b[t]]] = us_a[s]  # need something to make this sort of thing work

un.vector().set_local(un_a)
unb.vector().set_local(unb_a)

# save for viewing :
File("output/u.pvd")      << u
File("output/ub_n.pvd")   << unb
File("output/un.pvd")     << un
File("output/us.pvd")     << us

so the function 'un' above will have 'us' on it's bottom side, and zero everywhere else.

Here is 'us' plotted on top of a wireframe 'u', for reference :

I know that you can apply a boundary condition to the new function like this:

ff = FacetFunction('size_t', mesh, 0)
for f in facets(mesh):
  n = f.normal()
  if n.z() < 1e-8 and abs(n.x()) < 1e-8 and abs(n.y()) < 1e-8:
    ff[f] = 1

File("output/ff.pvd") << ff
bc = DirichletBC(V, u, ff, 1)

bc.apply(un.vector())

but we are working with SubMeshes and would like to efficiently translate to grand-parent meshes.

maartent · Answer 1 · Sep 11, 2015

I encountered a very similar problem, and I actually used this question and its links to work out a solution. If I understand your question correctly, you can try the following (mesh, bmesh, submesh, V, Vs, u and un defined as in your example):

un = Function(V)
un.vector().array().fill(0.0)
for Vs_dof, val in enumerate(us.vector().array()):
    submesh_vertex = dof_to_vertex_map(Vs)[Vs_dof]
    boundary_vertex = submesh.data().array('parent_vertex_indices', 0)[submesh_vertex]
    mesh_vertex = boundarymesh.entity_map(0)[int(boundary_vertex)] # np.uint not accepted
    V_dof = vertex_to_dof_map(V)[mesh_vertex]
    un.vector()[V_dof] = val

Vertex mapping from a SubMesh of a BoundaryMesh back to the actual Mesh!

1 Answer