Parallel Mesh Generation Using mshr and Legacy FEniCS

General
, ,
1

Hi everyone,

I’m using FEniCS 2019.1.0 with mshr. From what I understand, mshr is serial only - however, is there a way to modify meshes generated with mshr so that I can use MPI parallelization?

Consider as an example the rather trivial program below, where we solve the Poisson equation on a patterned surface:

import fenics as fe
import numpy as np
import mshr
 
comm = fe.MPI.comm_world
rank = fe.MPI.rank(comm)

fe.parameters["form_compiler"]["optimize"] = True
fe.parameters["form_compiler"]["cpp_optimize"] = True

L_x = 1
L_y = 1

surface_amplitude = 0.
surface_freq = L_x/(8*np.pi)
def surfaceExpr(x):
    
    return surface_amplitude*np.cos(surface_freq*(x - L_x/2) )

def surfExprDeriv(x):
    return - surface_amplitude*surface_freq*np.sin(surface_freq*(x-L_x/2))

domain_n_points = 60
domain_points = []
for n in range(domain_n_points + 1):
    x = n*L_x/domain_n_points
    domain_points.append(fe.Point(x, surfaceExpr(x) ))
domain_points.append(fe.Point(L_x,  surface_amplitude))
domain_points.append(fe.Point(L_x, L_y))
domain_points.append(fe.Point(0., L_y))
domain_points.append(fe.Point(0., surface_amplitude))
domain1 = mshr.Polygon(domain_points)
if rank == 0:
    mesh = mshr.generate_mesh(domain1, domain_n_points)

#mesh = mshr.generate_mesh(domain1, domain_n_points)

V = fe.FunctionSpace(mesh, "Lagrange", 1)


# Define boundary condition
u_D = fe.Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary
bc = fe.DirichletBC(V, u_D, boundary)
# Define variational problem
trialFn = fe.TrialFunction(V)
testFn = fe.TestFunction(V)
f = fe.Constant(-6.0)
a = fe.dot(fe.grad(trialFn), fe.grad(testFn))*fe.dx
L = f*testFn*fe.dx
# Compute solution
u = fe.Function(V)
fe.solve(a == L, u, bc)
2

You can save the mesh to disk on rank 0 with COMM_SELF using a parallel capable format (e.g., XDMF), and then reload it for parallel computation with all ranks using COMM_WORLD.