Parallel implementation using mesh.geometry.dofmap

dolfinx
, ,
1

Hi everyone,

I’m encountering an issue while implementing parallel processing in my code. It seems that when I try to compute an internal area enclosed by a surface mesh, there’s a problem when the function calls the indexes.

Below is my MWE:

import os

from mpi4py import MPI # type: ignore
from dolfinx.io import XDMFFile, gmshio # type: ignore

import basix.ufl # type: ignore
import dolfinx # type: ignore
import gmsh  # type: ignore
import numpy as np


def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.
    """

    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", h)

    circle = model.occ.addDisk(0, 0, 0, R, R, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0, gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)

def order_dofs(dof):
    """Sort the dofs in order.
    """

    len,_= dof.shape
    # print(type(dof))
    dof_ordered = np.empty_like(dof)
    dof_ordered[0,:] = dof[0,:]
    
    for i in range(len):
        pos = 0
        for j in range(len):
            if dof[j,0] == dof_ordered[i,1]:
                pos = j
                break
        if pos != 0:
            dof_ordered[i+1,:] = dof[pos,:]     
            
    return dof_ordered

def area(dofs,coordinates):
    """Compute the area with the dofs and coordinates.
    """

    num_ele, _ = dofs.shape

    area = 0
    count = 0

    for i in range(num_ele):
        x1 = dofs[i,0]
        x2 = dofs[i,2]
        x3 = dofs[i,1]        
        area += coordinates[x1,0]*coordinates[x2,1] - coordinates[x1,1]*coordinates[x2,0]
        area += coordinates[x2,0]*coordinates[x3,1] - coordinates[x2,1]*coordinates[x3,0]

    area  = 1/2 * abs(area)

    return area

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

h = 1e-2
R = 1

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

dof_ordered = order_dofs(mesh.geometry.dofmap)  # Order dofs only once
area_init = area(dof_ordered, mesh.geometry.x)   # Calculate area_init only once

The code works well when run sequentially, but when running in parallel, there is a conflict with the index handling across processes and I also tried implementing the following:

if MPI.COMM_WORLD.rank == 0
    dof_ordered = order_dofs(mesh.geometry.dofmap)  # Order dofs only once
    area_init = area(dof_ordered, mesh.geometry.x)   # Calculate area_init only once

However, this approach didn’t resolve the issue

I would greatly appreciate any help or suggestions on how to resolve this issue.

Best regards,

Best,