I’m trying to load a mesh from an XDMF and H5 file (created in Gmsh) using the following code:
from mpi4py import MPI
# Create the MPI communicator
comm = MPI.COMM_WORLD
import dolfinx # using dolfinx v0.7.2
import numpy as np
import h5py
# Read mesh from the XDMF file
with dolfinx.io.XDMFFile(comm, "collimator_tetra.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
print("mesh read")
cell_tags = xdmf.read_meshtags(mesh, name="Grid")
This produces the following error:
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
I’ve seen this previous question, where the user resolved the problem by adding all surfaces to a physical group. I’ve added all entities (including surfaces) to physical groups, but I’m still experiencing the same issue.
I’ve also tried to locate the source of the problem by starting with the simplest shapes and gradually adding shapes, physical groups, and fields until I get an error. The error only occurs after I define the field.
Here is the mesh that is causing the issue:
import gmsh # For creating mesh
import numpy as np # For adding magnets
import meshio # For converting to xdmf format for FEniCSx
from collections import defaultdict # For grouping unsorted entities
\# VARIABLES ------------------------------------------
run_gui = False # Show the mesh in GMSH GUI
\# Thruster / collimator variables:
\# Note that these are arbitrary values
r_c = 0.0125 # m, = 12.5mm, cathode radius
r_coll = 0.05 # m, = 3cm, collimator radius
r_a = 2\*r_c # m, = 25mm, anode radius
r_trig = 0.001 # m = 1mm, trigger radius
kappa = 0.03 # m, = 3cm, cathode-collimator distance
l_c = r_c # m = 12.5mm, length of cathode
l_coll = 0.08 # m, = 10cm, length of collimator
l_a = l_c+kappa # m, = 42.5mm, length of anode
l_trig = l_c \* 3/2 # m, = 18.75mm, length of trigger
n_mag = 3 # number of pairs of magnets
h_mag = 0.0025 # m, = 2.5mm, height of magnet
l_mag = l_coll # m, = 10cm, length of magnet
w_mag = 2\*h_mag # m, = 5mm, width of magnet
\# Plasma variables:
debye_length = 0.005
\# Model variables:
size_min = 0.001
size_max = debye_length
\# File variables:
meshname = "collimator.msh" # Mesh will be saved to this file
\# FUNCTIONS ------------------------------------------
def add_magnet_pair(model:gmsh.model, x:float, y:float, z:float, dx:float, dy:float, dz:float, angle:float, last_tag:int, magnet_surface_list:list, magnet_volume_list:list)->gmsh.model:
'''
Adds a pair of magnets at a given point on the collimator-thruster mesh.
The magnet is a rectangle defined by a point (x\', y\', z\') and the extents along the x-, y- and z-axes, and is rotated by angle.
The second magnet is added diagonally opposite the first.
Args:
model: the model to be updated
x: x-coordinate of one corner of box
y: y-coordinate of one corner of box
z: z-coordinate of one corner of box
dx: x-coordinate of other corner of box
dy: y-coordinate of other corner of box
dz: z-coordinate of other corner of box
angle: (radians) the rectangle is rotated about the x-axis by this angle
last_tag: the tag of the entity to add the magnets to
magnet_surface_list: list of dimtags of all surfaces of the magnets
magnet_volume_list: list of all the tags of the magnet volumes. For defining the magnets as a physical group later.
Returns:
model: updated gmsh model
last_tag: the tag of the updated model
magnet_surface_list: list of magnet surfaces by \[(dimension=2, tag)\] (for defining field later)
'''
\# Add first magnet
mag_1_vol = model.occ.addBox(x,-y+1e-6,-z, dx, dy, dz, tag=last_tag+1)
model.occ.rotate(\[(3, mag_1_vol)\], x+dx/2,0,0, 1,0,0, angle)
\# Add second magnet opposite the first
mag_2_vol = model.occ.addBox(x,-y+1e-6,-z, dx, dy, dz, tag=last_tag+2)
model.occ.rotate(\[(3, mag_2_vol)\], x+dx/2,0,0, 1,0,0, angle+np.pi)
gmsh.model.occ.synchronize() # Cut both at the same time - reduces synchronise calls.
mag_1_surf = model.getBoundary(\[(3, mag_1_vol)\], combined=False, oriented=False)
mag_2_surf = model.getBoundary(\[(3, mag_2_vol)\], combined=False, oriented=False)
magnet_surface_list.append(mag_1_surf)
magnet_surface_list.append(mag_2_surf)
magnet_volume_list.append(mag_1_vol)
magnet_volume_list.append(mag_2_vol)
gmsh.model.occ.synchronize()
\# Update tag
last_tag = last_tag+2
return model, last_tag, magnet_surface_list, magnet_volume_list
def add_physical_boundary_surface(model:gmsh.model, vol_tag:int|list, phys_tag=int, surf_name=str)->int:
'''
Adds a physical group containing the boundary of the volume vol_tag.
The resulting physical group has the tag phys_tag, and the name surf_name.
Args:
- model: gmsh model to add the physical group to.
- vol_tag: (positive integer) the tag of the volume whose surface will be added.
- phys_tag: (positive integer) the tag of the resulting physical group. -1 to let gmsh assign a tag automatically.
- surf_name: (string) the name of the physical group.
Returns: the tag of the physical group, phys_tag.
'''
if type(vol_tag) == int:
boundary_entities = model.getBoundary(\[(3, vol_tag)\])
surface_tags = \[tag for dim, tag in boundary_entities if dim == 2\]
return model.addPhysicalGroup(2, surface_tags, tag=phys_tag, name=surf_name)
elif type(vol_tag) == list:
surface_tags = \[\]
for vol in vol_tag:
boundary_entities = model.getBoundary(\[(3, vol)\])
for dim, tag in boundary_entities:
if dim == 2:
surface_tags.append(tag)
return model.addPhysicalGroup(2, surface_tags, tag=phys_tag, name=surf_name)
def create_mesh(mesh, cell_type, prune_z=False):
'''
Modified from code written by J.S. Dokken, December 2021, available here: (https://fenicsproject.discourse.group/t/importing-mesh-gives-error/7282/5)
This reads in a mesh of type gmsh.mesh, and returns a mesh consisting of only the cells of type cell_type.
Args:
- mesh: gmsh.mesh containing mixed cell types
- cell_type: either "tetra" or "triangle" - the cell type to be included in the output mesh.
Returns:
- mesh containing only cells of cell_type.
'''
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type) # For encoding facet tags (required to preserve physical groups)
points = mesh.points\[:, :2\] 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
\# DEFINING A MODEL WITH GMSH -------------------------
\# ---- Set parameters ----
gmsh.initialize()
gmsh.model.add("collimator") # Add model with name collimator
\# ---- Define geometry ----
vol_coll = gmsh.model.occ.addCylinder(l_a,0,0, l_coll,0,0, r_coll, 1) # Collimator cylinder, tag = 1
vol_a = gmsh.model.occ.addCylinder(0,0,0, l_a,0,0, r_a,2) # Anode cylinder, tag = 2
\# Then remove the cathode from the anode cylinder
vol_c = gmsh.model.occ.addCylinder(0,0,0, l_c,0,0, r_c, 4) # tag = 4
vol_trig = gmsh.model.occ.addCylinder(0,0,0, l_trig,0,0, r_trig, tag=5) # tag = 5
gmsh.model.occ.synchronize()
\# Add magnets:
if n_mag == 0:
pass
elif n_mag > 0:
last_tag = 5 # Tag of most-recently created entity
angle_array = np.linspace(0, np.pi, n_mag+1)\[0:n_mag\] # Generate array of equally spaced magnet pairs around spherical cathode surface
magnet_surface_list = \[\]
magnet_volume_list = \[\]
for loop_idx, angle in enumerate(angle_array):
gmsh.model, last_tag, magnet_surface_list, magnet_volume_list = add_magnet_pair(gmsh.model, l_a, r_coll, w_mag/2, l_mag, h_mag, w_mag, angle, last_tag, magnet_surface_list, magnet_volume_list)
gmsh.model.occ.synchronize()
flat_magnet_surface_tags = \[tag for surf in magnet_surface_list for (dim, tag) in surf\]
gmsh.model.mesh.embed(2, flat_magnet_surface_tags, 3, vol_coll)
\# DEFINING PHYSICAL GROUPS ----------------------------------
\# ---- Define physical groups -----
print(f"magnet surfaces: {magnet_surface_list}")
print(f"magnet volume: {magnet_volume_list}")
phys_col = gmsh.model.addPhysicalGroup(3, \[vol_coll\], 49, "collimator")
phys_a = gmsh.model.addPhysicalGroup(3, \[vol_a\], 50, "anode")
phys_c = gmsh.model.addPhysicalGroup(3, \[vol_c\], 51, "cathode")
phys_trig = gmsh.model.addPhysicalGroup(3, \[vol_trig\], 52, "trigger")
phys_mag = gmsh.model.addPhysicalGroup(3, magnet_volume_list, 53, "magnets")
phys_surf_col = add_physical_boundary_surface(gmsh.model, vol_coll, 54, "collimator_surface")
phys_surf_a = add_physical_boundary_surface(gmsh.model, vol_a, 55, "anode_surface")
phys_surf_c = add_physical_boundary_surface(gmsh.model, vol_c, 56, "cathode_surface")
phys_surf_trig = add_physical_boundary_surface(gmsh.model, vol_trig, 57, "trigger_surface")
phys_surf_mag = gmsh.model.addPhysicalGroup(2, flat_magnet_surface_tags, 58, "magnets_surface")
\# DEFINING FIELDS -------------------------------------------
# NOTE: when this section is commented out, the resulting xdmf / h5 file is able to be read without error by fenics.
\# ---- Define mesh field ----
print(f"magnet surfaces: {magnet_surface_list}")
\# Get magnet boundary
magnet_array = np.zeros((n_mag\*2\*6))
for mag_idx, magnet in enumerate(magnet_surface_list):
for surf_idx, surface in enumerate(magnet):
magnet_array\[6\*mag_idx+surf_idx\] = surface\[1\]
distance_field = gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "FacesList", magnet_array)
gmsh.model.mesh.field.setNumber(1, "Sampling", 1000)
threshold_field = gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", distance_field) # Create a mesh field that defines mesh size based on proximity to magnets
gmsh.model.mesh.field.setNumber(2, "SizeMin", size_min)
gmsh.model.mesh.field.setNumber(2, "SizeMax", size_max)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.05) # 0.5mm
gmsh.model.mesh.field.setAsBackgroundMesh(2)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
gmsh.model.occ.synchronize()
\# UNGROUPED ENTITIES PHYSICAL GROUPS ---------------------------------
entities = gmsh.model.get_entities(-1)
print("Entities")
for e in entities:
print(f"dim:{e\[0\]}, tag:{e\[1\]}")
\# Find entities not already part of a physical group and add them to a physical group by dimension
ungrouped_entities = entities
misctags = defaultdict(list)
physicalgroups = gmsh.model.get_physical_groups(-1)
print("Physical groups")
for p in physicalgroups:
print(f"{gmsh.model.get_physical_name(p\[0\], p\[1\])}, dim: {p\[0\]}, tag = {p\[1\]}")
grouped_entities = gmsh.model.get_entities_for_physical_group(p\[0\], p\[1\])
print(f"entities in {gmsh.model.get_physical_name(p\[0\], p\[1\])}: {grouped_entities}")
for g in grouped_entities:
dt = (p\[0\], int(g))
if dt in ungrouped_entities:
ungrouped_entities.remove(dt)
\# Sort ungrouped entities by dimension
for d, e in ungrouped_entities:
misctags\[d\].append(e)
phys_misc = defaultdict(int)
phys_names = \["points", "lines", "surfaces", "volumes"\]
for ii in range(0, 3):
if ii in misctags:
phys_misc\[ii\] = gmsh.model.addPhysicalGroup(ii, misctags\[ii\], tag=-1, name=phys_names\[ii\])
print(f"Physical group created for ungrouped {phys_names\[ii\]}, contains these tags: {misctags\[ii\]}")
for t in misctags\[ii\]:
dt = (ii, t)
if dt in ungrouped_entities:
ungrouped_entities.remove(dt)
\# Check again for ungrouped entities: (expect this to print an empty list)
print(f"Ungrouped entities after grouping: {ungrouped_entities}")
gmsh.model.occ.synchronize()
\# GENERATE MESH ------------------------------------------
\# ---- Parallelize ----
gmsh.option.setNumber("General.NumThreads", 5) # Use multithreading to generate mesh faster
gmsh.option.setNumber("Mesh.Algorithm3D", 4) # Frontal-Delaunay (default. If this causes issues consider changing to 10 - HXT which is a version of Delaunay that is optimised for parallel meshing).
\# ---- Generate Mesh ----
gmsh.model.mesh.generate() #(3) # Generate in 3d
gmsh.write(meshname) # Save mesh as meshname
\# ---- Output to GUI ----
if run_gui: # Off by default for remote server.
gmsh.fltk.run()
gmsh.finalize()
print("Done")
\# CONVERT TO XDMF FORMAT -------------------------------------
\# This step reads in the msh file so that the cell can be run without the use of gmsh.
in_mesh = meshio.read(meshname)
collimator_tetra = create_mesh(in_mesh, "tetra")
collimator_triangle = create_mesh(in_mesh, "triangle")
meshio.write("collimator_tetra.xdmf", collimator_tetra, compression=None)
meshio.write("collimator_tetra_facets.xdmf", collimator_triangle, compression=None)
For reference, when the field section is commented out, the resulting XDMF/H5 files can be read by FEniCSx without any problems.
Any help would be very much appreciated! 