Hello,
I am currently having some issues reading in a mesh with about 400,000 elements in parallel inside of fenicsX using XDMFFile. My code is as follows:
from petsc4py import PETSc
import numpy as np
from ufl import (TestFunction, TrialFunction, div, dot, dx, grad, nabla_grad, inner, outer, split, Identity, derivative, CellDiameter)
from basix.ufl import element, mixed_element
from dolfinx import la, default_scalar_type, log
from dolfinx.fem import (
Constant,
Expression,
Function,
dirichletbc,
DirichletBC,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile, VTXWriter
from mpi4py import MPI
from pathlib import Path
from petsc4py.PETSc import ScalarType
import adios4dolfinx
################# Read in Mesh
with XDMFFile(MPI.COMM_WORLD, "Aneur_400k_X.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
cd = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "Aneur_400k_X_facet.xdmf", "r") as xdmf:
fd = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
tdim = mesh.topology.dim
fdim = tdim - 1
################# Define Finite Elements
# Taylor-Hood Elements
K = 1
U_el = element("CG", mesh.basix_cell(), K + 1, shape = (mesh.geometry.dim,)) # Velocity Element
P_el = element("CG", mesh.basix_cell(), K) # Pressure Element
################# Define Function Space
W_el = mixed_element([U_el, P_el]) # Define Mixed element
W = functionspace(mesh, W_el) # Define Mixed Function Space
################# Define Trial and Test Functions
du = TrialFunction(W) # For Jacobian
up = Function(W)
u,p = split(up)
vq = TestFunction(W)
v,q = split(vq)
################# Define Function Spaces for Runge-Kutta Terms
V_sub, _ = W.sub(0).collapse()
u_1 = Function(V_sub)
################ Set up Initial Conditions using initial value for each RK Term
U_initial = 0.1 # Initial Speed
n_in = [-0.3419, 0.9251, 0.1651]
for i, vel in enumerate([-U_initial*n_in[0], -U_initial*n_in[1], -U_initial*n_in[2]]):
u_1.sub(i).interpolate(lambda x, vel_in= vel: np.full(x.shape[1], vel_in, dtype=np.float64))
#### Check if Velocity was assigned properly
folder = Path("../scratch")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "Vel_check.bp", [u_1], engine="BP4")
vtx_u.write(0.0)
My mesh was converted from a gmsh file to xdmf using:
from dolfinx.io import gmshio
from mpi4py import MPI
import meshio
import numpy as np
mesh, cell_markers, facet_markers = gmshio.read_from_msh("Aneur_400k_X.msh", MPI.COMM_WORLD, gdim=3)
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :3] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
return out_mesh
#if proc == 0:
# Read in mesh
msh = meshio.read("Aneur_400k_X.msh")
# Create and save one file for the mesh, and one file for the facets
tetra_mesh = create_mesh(msh, "tetra", prune_z=True) # Mesh domain
triangle_mesh = create_mesh(msh, "triangle", prune_z=True) # Mesh 2D surface
meshio.write("Aneur_400k_X", tetra_mesh)
meshio.write("Aneur_400k_X_facet.xdmf", triangle_mesh)
MPI.COMM_WORLD.barrier()When I read in the mesh in serial, everything is fine and the mesh look correct in the output:
However, when I read the mesh in parallel, there are missing cells, as can be seen in the photos below:
I have tried reading in simple, smaller meshes in parallel using the above code, and it works fine. However, I have been reading in multiple different meshes for my geometry in parallel and have continued to see missing cells.
I appreciate any help, thanks!