Reading a mesh in serial in order to access ordered contour coordinates

mesh
1

Hello

I have a slightly strange issue at the moment regarding parallelism within Fenics. I’m working on an adaptive re-meshing scheme and I’m using the coordinates of a contour (generated through matplotlib) to control the adaptivity. When run in serial, the code below returns an ordered list of coordinates that I can then use to generate a mesh using gmsh. However, when run in parallel this doesn’t work as the mesh is distributed across cores.

from fenics import *
from mpi4py import MPI
import matplotlib.pyplot as plt

comm = MPI.COMM_WORLD
rank = comm.rank

if rank == 0:

    nx = 5

    mesh= RectangleMesh(Point(0,0),Point(1,2),nx,2*nx) 

    V = FunctionSpace(mesh, 'CG', 2)

    eps = Constant(0.5*mesh.hmin()**0.9)

    dist = Expression('sqrt((pow((x[0]-A),2))+(pow((x[1]-B),2)))-r', degree=2, A=0.5,B=0.5,r=0.25)

    dist2 = Expression('(1/(1+exp((dist/eps))))', degree=2, eps=eps, dist=dist)

    phi0 = interpolate(dist2, V)

    cs = plot(phi0 , mode='contour', levels=[0.5])
    plt.close()

    a=[]
    for item in cs.collections:
        for i in item.get_paths():
            a.append(i.vertices)
else:

    a = None

coords = comm.bcast(a, root = 0)

The code above is what I would like to have working, but it just seems to freeze when ran in the terminal. The contour coordinates are calculated on one processor so they are in an ordered list, then broadcast to all other cores so the rest of the code can run in parallel.

Thanks! :slight_smile:

2

If you’re solving the problem on process 0, then you should pass the appropriate communicator (e.g. COMM_SELF) to the Mesh and any other objects that are innately parallelised. Otherwise process 0 will just be waiting for all other processes to take part in the parallelisation whilst they do nothing.

3

Hi Nate,

Thanks for the help, I did not know that I could use COMM_SELF in the argument of a mesh. I will post the correct code here if anyone else needs help.

from fenics import *
from mpi4py import MPI
import matplotlib.pyplot as plt

comm = MPI.COMM_WORLD
rank = comm.rank

nx = 5

mesh = RectangleMesh(MPI.COMM_SELF, Point(0,0),Point(1,2),nx,2*nx) 

V = FunctionSpace(mesh, 'CG', 2)

eps = Constant(0.5*mesh.hmin()**0.9)

dist = Expression('sqrt((pow((x[0]-A),2))+(pow((x[1]-B),2)))-r', degree=2, A=0.5,B=0.5,r=0.25)

dist2 = Expression('(1/(1+exp((dist/eps))))', degree=2, eps=eps, dist=dist)

phi0 = interpolate(dist2, V)

if rank == 0:

    cs = plot(phi0 , mode='contour', levels=[0.5])
    plt.close()

    a=[]
    for item in cs.collections:
        for i in item.get_paths():
            a.append(i.vertices)
else:
    
    a = None

b = comm.bcast(a, root = 0)
print(rank, b)

Thanks again! :slight_smile:

1 Like