MeshView dofs distribution on parallel

[javijv4]

Hello,

I have a problem that requires to define a submesh of different topology and then relate the dofs of each mesh with each other (considering that the same vertex at each mesh has the same dofs). This works fine in serial, however, when I moved to parallel I found some issues on how MeshView deals with the shared vertex dofs. Here is an MWE of a rectangle mesh, where the submesh is just the upper and lower boundary:

import dolfin as d
import numpy as np
nx, ny = 8, 6

#%% Parent mesh
mesh = d.RectangleMesh(d.Point(0., -1), d.Point(10., 1.), nx, ny)
d.File('mesh/mesh.pvd') << mesh

comm = mesh.mpi_comm()
rank = comm.rank

V = d.VectorFunctionSpace(mesh, 'CG', 1)
dofmap = V.dofmap()
u = d.Function(V)

dofs = d.Function(V)
dofs.rename('dofs', 'dofs')
dofs.vector().set_local(np.arange(dofs.vector().local_size())+ dofmap.ownership_range()[0])
d.File('dofs.pvd') << dofs
#%% Submesh
class Surface(d.SubDomain):
    def inside(self, x, on_boundary):
        return (d.near(x[1], -1.) or d.near(x[1], 1)) and on_boundary
    
surf = Surface()
marker = d.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
surf.mark(marker, 1)

submesh = d.MeshView.create(marker, 1)
SV = d.VectorFunctionSpace(submesh, 'CG', 1)
sdofmap = SV.dofmap()
us = d.Function(SV)

sdofs = d.Function(SV)
sdofs.rename('sdofs', 'sdofs')
sdofs.vector().set_local(np.arange(sdofs.vector().local_size()) + 
                         sdofmap.ownership_range()[0])
d.File('sdofs.pvd') << sdofs

#%% Unowned dofs
print(rank, 'parent', dofmap.local_to_global_unowned())
print(rank, 'submesh', sdofmap.local_to_global_unowned())

If I run these MWE with two processors, I have 7 shared vertex in the parent mesh and 2 in the submesh (one in the upper boundary, one in the lower). I would expect that, if the dofs of the shared upper node on the parent mesh belongs to the partition 0, then the dofs of this node on the submesh belongs to this partition too. This way I can relate all the dofs of one parent partition with the dofs of the submesh partition. However, this is not the case. If I plot the dofs at each partition, I get something like these (if the value is 0 it means that the dofs doesn’t belong to the partition):

Here you can see that the X dof of the upper shared node belongs to the left partition in the submesh whereas it belongs to the right partition in the parent mesh. This doesn’t happen always, for example, if I use a coarser mesh (nx=ny=4) I get the desired behavior:

In summary, I want to be able to relate the dofs of a submesh to the parent mesh inside a partition. Right now, I can’t do that because the shared vertex dofs are not distributed equally in the parent and the submesh. It wasn’t so easy to explain, so I hope I made it clear.

Thank you!

[Fvergnet]

I encounter the same issue… Is this a bug in meshview ? Does anybody has a solution ?

[nate]

This problem you’re having isn’t immediately clear to me.

Are you sure your issue isn’t just that you haven’t updated the ghost values?

E.g.

dofs.vector().set_local(np.arange(dofs.vector().local_size())+ dofmap.ownership_range()[0])
dofs.vector().update_ghost_values()

[Fvergnet]

Hi @nate,

I already do what you suggest (i think) with u.vector().apply("insert"). The problem as I see it is that the mesh partitioning is not the same between the meshview and the parent mesh.

May be it will be more explicit with an example. Below I construct two P1 functions on the mehview and the parent mesh and for each dof I insert the proc rank value. When running the code in parallel with 3 proc, you should see that the “boundaries” between proc do not match (see screenshots attached).

To me, this is an issue if one wants to construct a dofmap between the meshview and the parent mesh (but may be you know a workaround for that).

PS: I run the code in quay.io/fenicsproject/dev docker container.

from dolfin import *

# Parameters
N = 30
xc = .5
yc = .5
rc = .4
rank = MPI.comm_world.rank

mesh = UnitSquareMesh(N,N)

# Create meshfunction to tag subdomain (a disk of center (xc,yc) and radius rc)
mf = MeshFunction("size_t",mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
    if (c.midpoint()[0]-xc)**2 + (c.midpoint()[1]-yc)**2 <= rc**2:
        mf.set_value(c.index(), 1)

# Create meshview
submesh = MeshView.create(mf, 1)

# Create rank function on mesh
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)
rank_func_mesh = Function(V)
rank_func_mesh.rename("rank on mesh","rank on mesh")
ndofs = rank_func_mesh.vector().vec().getArray().size
for i in range(ndofs):
    rank_func_mesh.vector().vec().getArray(readonly=0)[i] = rank
rank_func_mesh.vector().update_ghost_values()
#rank_func_mesh.vector().apply("insert") 

# Create rank function on submesh
P1sub = FiniteElement("CG", submesh.ufl_cell(), 1)
Vsub = FunctionSpace(submesh, P1sub)
rank_func_submesh = Function(Vsub)
rank_func_submesh.rename("rank on submesh","rank on submesh")
subndofs = rank_func_submesh.vector().vec().getArray().size
for i in range(subndofs):
    rank_func_submesh.vector().vec().getArray(readonly=0)[i] = rank
rank_func_submesh.vector().update_ghost_values()
#rank_func_submesh.vector().apply("insert")

# Save
with XDMFFile(MPI.comm_world, "rank_on_mesh.xdmf") as file:
    file.write(rank_func_mesh,XDMFFile.Encoding.HDF5)

with XDMFFile(MPI.comm_world, "rank_on_submesh.xdmf") as file:
    file.write(rank_func_submesh,XDMFFile.Encoding.HDF5)

[Fvergnet]

Any clues ? I am a bit stuck for the construction of the dofmap between the meshview and the parent mesh in parallel…