[copied from ALLANSWERED]
Hi, my question is related to this old forum post: https://fenicsproject.org/qa/6810/vertex-mapping-from-submesh-boundarymesh-back-actual-mesh/
The link above shows how to map function values between a mesh, its boundarymesh, and a submesh of the boundarymesh. For convenience, the code below is an updated and cleaned up version of the question and answer from that old post. This (working) code maps functions from the submesh to the boundarymesh, as well as from the submesh to the grandparent mesh (though a little less elegantly).
from dolfin import *
mesh = UnitCubeMesh(10, 10, 10) # this will be the "grandparent" mesh
mesh.coordinates()[:,0] -= .5 # shift x-coords
mesh.coordinates()[:,1] -= .5 # shift y-coords
bmesh = BoundaryMesh(mesh, "exterior") # surface boundary mesh
# mark the cells on the bottom of the bmesh surface by iterating over bmesh cells.
# Look up ccorresponding facet in mesh and mark if facet normal points in -z direction
cellmap = bmesh.entity_map(2)
vertmap = bmesh.entity_map(0)
pb = MeshFunction("size_t", bmesh, dim =1)
for c in cells(bmesh):
if Facet(mesh, cellmap[c.index()]).normal().z() < 0:
pb[c] = 1
# use the marked bottom of the bmesh to create a Submesh
submesh = SubMesh(bmesh, pb, 1) # bottom of boundary mesh
# FunctionSpaces on main mesh, bmesh, submesh
V = FunctionSpace(mesh, "CG", 1) # mesh function space
Vb = FunctionSpace(bmesh, "CG", 1) # surface function space
Vs = FunctionSpace(submesh, "CG", 1) # submesh function space
# 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)
t = submesh.data().array('parent_vertex_indices', 0) # mapping from submesh back to bmesh
# Functions on main mesh, bmesh, and submesh
u = Function(V)
ub = Function(Vb) # boundary function
us = Function(Vs) # surface function
# Interpolate the following expr onto u, ub, us
expr = Expression('sqrt(pow(x[0],2) + pow(x[1], 2))', degree=2)
u.interpolate(expr)
ub.interpolate(expr)
us.interpolate(expr)
# Some empty function to test the mappings
ub_test = Function(Vb) # empty bmesh function
u_test = Function(V) # empty mesh function
# Mapping from submesh to bmesh (works)... Any way to avoid calls to get_local/set_local??
us_a = us.vector().get_local() # origin array
ub_test_a = ub_test.vector().get_local() # destination array
ub_test_a[b[t]] = us_a[s] # transfer
ub_test.vector().set_local(ub_test_a) # set destination function values
# Mapping from submesh to grandparent mesh (less than satisfying solution to question in fenics forum link)
# Bonus points for getting this kind of map composition to work:
# u_test_a = u_test.vector().get_local() # destination array
# u_test_a[m[b[t]]] = us_a[s] # transfer ( not working )
# u_test.vector().set_local(u_test_a)
for Vs_dof, val in enumerate(us.vector().get_local()):
submesh_vertex = si[Vs_dof]
boundary_vertex = t[submesh_vertex]
mesh_vertex = vertmap[int(boundary_vertex)] # np.uint not accepted
V_dof = m[mesh_vertex]
u_test.vector()[V_dof] = val
u.rename('u','u')
ub_test.rename('ub_test','ub_test')
u_test.rename('u_test','u_test')
us.rename('us','us')
File("output/u.pvd") << u
File("output/ub_test.pvd") << ub_test
File("output/u_test.pvd") << u_test
File("output/us.pvd") << us​
How can I do the same when V, Vb, and Vs are vector-valued?
# How to do the above when V, Vb, Vs are vector valued?
# FunctionSpaces on main mesh, bmesh, submesh
V = VectorFunctionSpace(mesh, "CG", 1, dim=2) # mesh function space
Vb = VectorFunctionSpace(bmesh, "CG", 1, dim=2) # surface function space
Vs = VectorFunctionSpace(submesh, "CG", 1, dim=2) # submesh function space
Obviously, there is no longer a one-to-one correspondence between vertex and dof. I’ve heard it’s possible to split the functions into components and use FunctionAssigners, but I cannot find any relevant documentation.
I tried this approach, but it is not robust (fails for certain meshes). https://www.allanswered.com/post/jewrb/why-does-this-mapping-of-vector-valued-function-dofs-between-boundarymesh-and-parent-mesh-fail/
Thank you!