How to map dofs of VECTOR functions between meshes, boundarymeshes, submeshes?

[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!

I tried modifying your earlier attempt from the AllAnswered post, to make it more robust against epsilons:

import dolfin as d

def make_mapping(sub_space, super_space):
   # Degrees of freedom
    dofs = super_space.dofmap().dofs()
    d_ = d.Function(super_space)
    d_.vector()[:] = dofs
    mapping = d.interpolate(d_, sub_space)

    # Removed this line #######################################
    #mapping.vector()[:] = mapping.vector()[:].astype('int')
 
    return mapping


# test the mapping for various mesh discretizations
for discr in range(1,30):
    
    # make a mesh, boundarymesh
    mesh = d.UnitSquareMesh(discr, discr)
    bmesh = d.BoundaryMesh(mesh, 'exterior')

    
    # Some function spaces corresponding to the above meshes    
    M = d.VectorFunctionSpace(mesh, 'CG', 1, dim=2)
    B = d.VectorFunctionSpace(bmesh, 'CG', 1, dim=2)
    
    # Expression with which to test mappings
    coord_expression = d.Expression(('x[0]','x[1]'), degree=1)
    
    coord_Mfun = d.Function(M)
    coord_Mfun.assign(d.interpolate(coord_expression, M))
    
    coord_Bfun = d.Function(B)
    coord_Bfun.assign(d.interpolate(coord_expression, B))   
    
    M_2_B = make_mapping(B, M)
    
    # Test M to B    
    for i, dof in enumerate(M_2_B.vector().get_local()):

        # Changed conversion to index here #################################
        
        # Add 0.5, then round down for a robust conversion from float to
        # non-negative integer:
        index = int(dof+0.5)
        
        try:
            assert(coord_Mfun.vector()[index]
                   == coord_Bfun.vector().get_local()[i])
        except:
            print('=='*16)
            print('i: ', i, 'dof: ', index, 'discr: ' , discr)
            print(coord_Mfun.vector()[index] ,
                  coord_Bfun.vector().get_local()[i])

This seems to pass all of the assertions.

3 Likes

I tried almost everything I could think of, besides checking my basic understanding of type conversions…
Thank you!!

On a side note, I’m a bit surprised there isn’t built in support for type of operation. Maybe I’ll get around to making a feature request or at least cleaning this up and posting as an undocumented demo.

As far as I can see, this implementation does not work in parallel.
I’ll just post an alternative solution for conversion between BoundaryMesh and its parent mesh, for scalar and vector “CG” functions. This version uses the FunctionAssigner to go from Scalar to Vector spaces. This function has parallel support.

import numpy
import dolfin

def SyncSum(vec):
    """ Returns sum of vec over all mpi processes.
    
    Each vec vector must have the same dimension for each MPI process """

    comm = dolfin.MPI.comm_world
    NormalsAllProcs = numpy.zeros(comm.Get_size()*len(vec), dtype=vec.dtype)
    comm.Allgather(vec, NormalsAllProcs)
        
    out = numpy.zeros(len(vec))
    for j in range(comm.Get_size()):
        out += NormalsAllProcs[len(vec)*j:len(vec)*(j+1)]
    return out


def mesh_to_boundary(v, b_mesh):
    """
    Returns a boundary interpolation of the CG-1 function v
    """

    # Extract the underlying volume and boundary meshes
    mesh = v.function_space().mesh()

    # We use a Dof->Vertex mapping to create a global 
    # array with all DOF values ordered by mesh vertices
    DofToVert = dolfin.dof_to_vertex_map(v.function_space())
    VGlobal = numpy.zeros(v.vector().size())

    vec = v.vector().get_local()
    for i in range(len(vec)):
        Vert = dolfin.MeshEntity(mesh, 0, DofToVert[i])
        globalIndex = Vert.global_index()
        VGlobal[globalIndex] = vec[i]
    VGlobal = SyncSum(VGlobal)

    # Use the inverse mapping to se the DOF values of a boundary
    # function
    surface_space = dolfin.FunctionSpace(b_mesh, "CG", 1)
    surface_function = dolfin.Function(surface_space)
    mapa = b_mesh.entity_map(0)
    DofToVert = dolfin.dof_to_vertex_map(dolfin.FunctionSpace(b_mesh, "CG", 1))

    LocValues = surface_function.vector().get_local()
    for i in range(len(LocValues)):
        VolVert = dolfin.MeshEntity(mesh, 0, mapa[int(DofToVert[i])])
        GlobalIndex = VolVert.global_index()
        LocValues[i] = VGlobal[GlobalIndex]

    surface_function.vector().set_local(LocValues)
    surface_function.vector().apply('')
    return surface_function



def boundary_to_mesh(f, mesh):
    """ Take a CG1 function f defined on a surface mesh and return a
    volume vector with same values on boundary but zero in volume
    """
    b_mesh = f.function_space().mesh()
    SpaceV = dolfin.FunctionSpace(mesh, "CG", 1)
    SpaceB = dolfin.FunctionSpace(b_mesh, "CG", 1)
    
    F = dolfin.Function(SpaceV)
    LocValues = numpy.zeros(F.vector().local_size())
    GValues = numpy.zeros(F.vector().size())


    map = b_mesh.entity_map(0) # Vertex map from boundary mesh to parent mesh
    d2v = dolfin.dof_to_vertex_map(SpaceB)
    v2d = dolfin.vertex_to_dof_map(SpaceV)

    dof = SpaceV.dofmap()
    imin, imax = dof.ownership_range()


    for i in range(f.vector().local_size()):
        GVertID = dolfin.Vertex(b_mesh, d2v[i]).index() # Local Vertex ID for given dof on boundary mesh
        PVertID = map[GVertID] # Local Vertex ID of parent mesh
        PDof = v2d[PVertID] # Dof on parent mesh
        value = f.vector()[i] # Value on local processor
        GValues[dof.local_to_global_index(PDof)] = value
    GValues = SyncSum(GValues)


    F.vector().set_local(GValues[imin:imax])
    F.vector().apply("")
    return F



mesh = dolfin.UnitSquareMesh(50,50)
b_mesh = dolfin.BoundaryMesh(mesh, "exterior")
def mesh_to_bmesh():
    V = dolfin.VectorFunctionSpace(mesh, "CG", 1)
    v_org = dolfin.project(dolfin.Expression(("x[0]+2", "pow(x[1]-x[0], 2)"), element=V.ufl_element()), V)
    vbs = []
    v_is = v_org.split(deepcopy=True)
    for v_i in v_is:
        vbs.append(mesh_to_boundary(v_i, b_mesh))

    Vb = dolfin.VectorFunctionSpace(b_mesh, "CG", 1)
    split_to_vec = dolfin.FunctionAssigner(Vb, [vi.function_space() for vi in vbs])
    vb_new = dolfin.Function(Vb)
    split_to_vec.assign(vb_new,vbs)
    
    outfile_1 = dolfin.XDMFFile("v_org.xdmf")
    outfile_1.write(v_org)
    outfile_1.close()
    
    outfile_15 = dolfin.XDMFFile("vb.xdmf")
    outfile_15.write(vb_new)
    outfile_15.close()
    assert(numpy.isclose(dolfin.assemble(dolfin.inner(vb_new, vb_new)*dolfin.dx),
                         dolfin.assemble(dolfin.inner(v_org, v_org)*dolfin.ds)))
    return v_org,vb_new

mesh_to_bmesh()

def bmesh_to_mesh():
    v_origin, vb_org = mesh_to_bmesh()
    vb_is = vb_org.split(deepcopy=True)
    v_vol = []
    for vb_i in vb_is:
        v_vol.append(boundary_to_mesh(vb_i, v_origin.function_space().mesh()))
    split_to_vec = dolfin.FunctionAssigner(v_origin.function_space(), [v_i.function_space() for v_i in v_vol])
    v_new = dolfin.Function(v_origin.function_space())
    split_to_vec.assign(v_new, v_vol)
    outfile_2 = dolfin.XDMFFile("v_new.xdmf")
    outfile_2.write(v_new)
    outfile_2.close()
    assert(dolfin.assemble(dolfin.inner(v_origin, v_origin)*dolfin.ds) ==
           dolfin.assemble(dolfin.inner(v_new, v_new)*dolfin.ds))
    
    assert(numpy.isclose(dolfin.assemble(dolfin.inner(vb_org, vb_org)*dolfin.dx),
                         dolfin.assemble(dolfin.inner(v_new, v_new)*dolfin.ds)))
    
bmesh_to_mesh()
3 Likes

@dokken Hey thanks so much for this solution and your help! Would you recommend using the same approach for mapping between the boundary mesh and a submesh (of the boundary mesh)?

The current implementation of SubMesh does not support parallel and from what I know, that’s not planned. So I am afraid you cannot use the same approach with a submesh.
However, we are developing a new MeshView class to create submeshes, which holds the mappings that you would need to use the approach suggested by @dokken with a submesh.
We are currently finalizing the development of this feature and trying to make it available as soon as possible.

4 Likes

@cdaversin Thank you for your message. When you say that “SubMesh does not support parallel”, do you mean that any code using SubMesh cannot be reliably “run” on multiple cores, or just that submeshes can only be assigned to one process?

I am solving an FSI problem using ALE methods in which a (small) portion of the domain boundary is a membrane structure. Hence, pressure forces are applied to the submesh from a function defined on the main mesh (as discussed above). I have been moving the boundarymesh with displacements defined on the submesh, and then calling ALE.move(mesh, boundarymesh). Except for the mapping issues discussed above, this seems to work fine in parallel. Is there another way to achieve this without the use of SubMesh? Thanks!

@niewiarowski When constructed, a SubMesh object builds a parent-to-submesh mapping for its vertex indices. This mapping is needed if you want to use the approach previously discussed (it corresponds to the entity_map(0) in the code sent by @dokken ).
The SubMesh implementation specifies explicitly that this is not working in parallel (see SubMesh.cpp, l. 195) :

    if (MPI::size(mesh.mpi_comm()) > 1)
      not_working_in_parallel("SubMesh::init");

It might work if your submesh is owned by a single process, but I am not quite sure.
As I mentioned, another way to do this would be to use the new MeshView class (still under development) instead of SubMesh. This has not been merged yet in the master branch, but you can look into my branch cecile/mixed-dimensional in Dolfin if you’re interested.

2 Likes

This does not help, but I stumbled across a similar mapping problem and I found this other way of implementing make_mapping that I found being extremely fast even for large 3D problems (still serial only though). I am posting this just in case you find this helpful.

from scipy.spatial import cKDTree
def make_mapping(sub_space, super_space):
    super_dof_coor = super_space.tabulate_dof_coordinates()
    sub_dof_coor = sub_space.tabulate_dof_coordinates()

    tree = cKDTree(super_dof_coor)
    _,mapping = tree.query(sub_dof_coor, k=1)
    return mapping
1 Like

@cdaversin Is this new MeshView will suport creating submeshes on the same domain but with coarser mesh (still aligned with parent)? Similar to built-in refine(), but as far as I know, refine() does not support parallel. Thanks!

really cool, thank you!

MeshView creates a new Mesh object built as a set of entities of the parent mesh, whose topology contains a mapping giving the relationship between cells numbering in the submesh and in the parent mesh. I don’t know much about the refine() function, but the dedicated demo demo/undocumented/refinement/demo_refinement.py seems to work in parallel.
So I guess you could use refine on the submesh built from MeshView, and the mapping submesh - parent can be used if you want to refine a specific region. But I haven’t tried, so this needs to be checked.