Issues when running in parallel (mixed-dimensional branch)

1

Dear all,

When using the mixed-dimensional branch, some errors appear in the building of mappings between parent meshes when running in parallel. The problem seems related to how the mesh is being partitioned, i.e., if the partition in the interface is coherent, the code runs without troubles, if not then the message Error in building the mapping (21, 15): Index not found appears freezing the execution.

For instance, when running using 1, 3, or 5 processes everything works perfectly, but fails when running with 2 or 4.

A MWE to reproduce this behavior is the following:

from dolfin import *

class Dirichlet(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 2) \
        or near(x[1], 0) \
        or near(x[1], 1)

class Neumann(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0)

class Interface(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 1)

# Mesh generation
n = 20
mesh = RectangleMesh(Point(0,0,0),Point(2,1,0),2*n,n,"crossed")

# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
  domain[c] = c.midpoint().x() > 1.0

# Interface markers
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)

# Meshes
mesh_L = MeshView.create(domain, 0)
mesh_R = MeshView.create(domain, 1)
mesh_I = MeshView.create(facets, 1)

# Function spaces
FE = VectorElement("CG", mesh_I.ufl_cell(), 1, dim=2) # Lagrange multiplier
VE = VectorElement("CG", mesh.ufl_cell(), 1, dim=2)   # Displacement
V_L  = FunctionSpace(mesh_L,VE)    # Displacement on left domain
V_R  = FunctionSpace(mesh_R,VE)    # Displacement on right domain
V_I_disp    = FunctionSpace(mesh_I,FE)   # Lagrange multiplier on interface
V_I_stress  = FunctionSpace(mesh_I,FE)   # Lagrange multiplier on interface
W = MixedFunctionSpace(V_L, V_R, V_I_disp)

# Markers for Dirichlet and Neuman bcs
ff_L = MeshFunction("size_t", mesh_L, mesh_L.geometry().dim()-1, 0)
Dirichlet().mark(ff_L, 1)
Neumann().mark(ff_L, 2)
Interface().mark(ff_L, 3)
ff_R = MeshFunction("size_t", mesh_R, mesh_R.geometry().dim()-1, 0)
Dirichlet().mark(ff_R, 1)
Interface().mark(ff_R, 3)

# Boundary conditions
bc_L = DirichletBC(W.sub_space(0), Constant((0,0)), ff_L, 1)
bc_R = DirichletBC(W.sub_space(1), Constant((0,0)), ff_R, 1)
bcs = [bc_L, bc_R]

# Elasticity parameters
rho   = Constant(2200)
lmbda = Constant(7.428e+04)
mu    = Constant(7.428e+04)

# Neumann condition
g = Constant((1e+5, 0))

# Measures
dx_L  = Measure("dx", domain=W.sub_space(0).mesh())
ds_L  = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=ff_L)
dx_R  = Measure("dx", domain=W.sub_space(1).mesh())
ds_R  = Measure("ds", domain=W.sub_space(1).mesh(), subdomain_data=ff_R)
ds_I  = Measure("dx", domain=W.sub_space(2).mesh())

# Function and test functions
(ul, ur, l_disp) = TrialFunctions(W)
(wl, wr, e_disp) = TestFunctions(W)

# Stress tensor
def sigma(r):
    return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))

# Variational problem
a_left = inner(sigma(ul),sym(grad(wl)))*dx_L
L_left = inner(g, wl)*ds_L(2)

a_right = inner(sigma(ur),sym(grad(wr)))*dx_R
L_right = inner(Constant((0,0)), wr)*dx_R

a = a_left + a_right \
  + inner(l_disp, wr-wl)*ds_I + inner(ur-ul, e_disp)*ds_I
L = L_left + L_right

# Solve
u_mixed = Function(W)
solve(a == L, u_mixed, bcs, solver_parameters={"linear_solver":"direct"})

Someone knows why this happens?

Thanks in advance!

2

Dear Hernan,

Sorry for my late reply. I can reproduce your issue, it might be related with the partitioner.
The error you get happens when building the mapping between two submeshes, (at least) one cell of the mesh with id 21 cannot be found in the submesh with id 15.
Since it doesn’t happen with 1,3,5 procs, my guess is that it could be related with the partitioner, but I don’t see (yet) how to fix it.

1 Like