I’m trying to obtain the elements (cells) and nodes (points) from a gmsh file to manipulate them and then create a dolfinx mesh. However, I’m having trouble with the dolfinx.mesh.create_mesh function, which I’ve seen has been subjected to changes (I’m using dolfinx version 0.8).
This is my simple code so far:
import numpy as np
from mpi4py import MPI
import ufl
import basix.ufl
import dolfinx
from dolfinx.io import XDMFFile
import meshio
def create_custom_mesh_3d(cells, x):
# Ensure cells and x are numpy arrays with correct types
cells = np.asarray(cells, dtype=np.int64, order='C')
x = np.asarray(x, dtype=np.float64, order='C')
# Debugging statements
print(f"cells shape: {cells.shape}")
print(f"x shape: {x.shape}")
# Ensure x has the correct shape (N, 3)
if x.ndim != 2 or x.shape[1] != 3:
raise ValueError("The array x should be a 2D array with shape (N, 3).")
# Ensure cells has the correct shape (num_cells, 4)
if cells.ndim != 2 or cells.shape[1] != 4:
raise ValueError("The array cells should be a 2D array with shape (num_cells, 4).")
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(x.shape[1],)))
#c_el = dolfinx.fem.coordinate_element(dolfinx.mesh.CellType.tetrahedron, 1)
custom_mesh = dolfinx.mesh.create_mesh(MPI.COMM_SELF, cells, c_el, x)
return custom_mesh
def write_mesh_to_xdmf(mesh, filename):
# Write the mesh to an XDMF file
with XDMFFile(mesh.comm, filename, "w") as file:
file.write_mesh(mesh)
mesh_file_name = "sphere"
# Read the .vtk file using meshio
meshio_mesh = meshio.read(mesh_file_name + ".msh", file_format="gmsh")
# Extract data
cells = [cell.data for cell in meshio_mesh.cells if cell.type == "tetra"]
points = meshio_mesh.points
# TODO: MANIPULATE CELLS AND POINTS HERE WHEN CODE WORKS
# Convert cells to the required shape
if len(cells) > 0:
cells = np.concatenate(cells)
# Create the mesh
custom_mesh = create_custom_mesh_3d(cells, points)
print(custom_mesh)
write_mesh_to_xdmf(custom_mesh, "sphere_custom.xdmf")
When I run it, I’m getting this error:
cells shape: (683, 4)
x shape: (204, 3)
Traceback (most recent call last):
File "/home/cborau/fenicsx_projects/test7.py", line 50, in <module>
custom_mesh = create_custom_mesh_3d(cells, points)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/cborau/fenicsx_projects/test7.py", line 29, in create_custom_mesh_3d
custom_mesh = dolfinx.mesh.create_mesh(MPI.COMM_SELF, cells,c_el, x )
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/cborau/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/mesh.py", line 401, in create_mesh
gdim = x.shape[1]
~~~~~~~^^^
IndexError: tuple index out of range
And this is the code I use to create the sphere mesh:
from mpi4py import MPI
import gmsh
def gmsh_sphere(model: gmsh.model, name: str) -> gmsh.model:
"""Create a Gmsh model of a sphere.
Args:
model: Gmsh model to add the mesh to.
name: Name (identifier) of the mesh to add.
Returns:
Gmsh model with a sphere mesh added.
"""
model.add(name)
model.setCurrent(name)
sphere = model.occ.addSphere(0, 0, 0, 1, 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=3, tags=[sphere])
# Generate the mesh
model.mesh.generate(dim=3)
return model
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model = gmsh.model()
model = gmsh_sphere(model, "Sphere")
model.setCurrent("Sphere")
gmsh.write('sphere.msh')
Thanks in advance!