Kernel Crashes When Loading 1D Mesh with XDMFFile

General
1

Hello, FEniCSx developers and users,

I am a structural engineering researcher who finds this program extremely useful.

In our field, it is common to model rebar as 1D elements, so I have been creating 1D elements for rebar within concrete and performing FEM analysis.

However, when I try to load a file containing 1D elements as shown below, I frequently encounter the message “restarting kernel…” and Spyder shuts down. Could there be an issue with how I am saving and loading the mesh?

I have attached the mesh-related code I wrote below.

    def create_mesh(self, 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[:, :2] if prune_z else mesh.points
        out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                               "name_to_read": [cell_data]})
        return out_mesh
        
    def XDMF_generator(self):
        # read msh file
        self.mesh = meshio.read("mesh3D.msh")
        # Save 3D mesh
        if "3D" in self.dimension_list:
            meshio.write("mesh3D.xdmf" , self.create_mesh(self.mesh,"tetra"))
            # meshio.write("mesh3D.xdmf" , meshio.Mesh(points=self.mesh.points,
            #                                            cells={"tetra": self.mesh.cells_dict["tetra"]},
            #                                            cell_data={"tetra_group": [self.mesh.cell_data_dict["gmsh:physical"]["tetra"]]}))
            
        # Save 1D mesh
        if "1D" in self.dimension_list:
            meshio.write("mesh1D.xdmf", self.create_mesh(self.mesh,"line"))
        MWV.comm.barrier()

with XDMFFile(MPI.COMM_WORLD, "mesh3D.xdmf", "r") as xdmf:
    domain_1 = xdmf.read_mesh(name="Grid")
    #Save cell tag
    cell_marker_1 = xdmf.read_meshtags(domain_1, name="Grid")
xdmf.close()

with XDMFFile(MPI.COMM_WORLD, "mesh1D.xdmf", "r") as xdmf:
    domain_2 = xdmf.read_mesh(name="Grid")
    #Save cell tag
    edge_marker_2 = xdmf.read_meshtags(domain_2, name="Grid")
xdmf.close()

The xdmf file of mesh 1D is

<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="1804 3" Format="HDF" Precision="8">mesh1D.h5:/data0</DataItem></Geometry><Topology TopologyType="Polyline" NumberOfElements="24" NodesPerElement="2"><DataItem DataType="Int" Dimensions="24 2" Format="HDF" Precision="8">mesh1D.h5:/data1</DataItem></Topology><Attribute Name="name_to_read" AttributeType="Scalar" Center="Cell"><DataItem DataType="Int" Dimensions="24" Format="HDF" Precision="8">mesh1D.h5:/data2</DataItem></Attribute></Grid></Domain></Xdmf>
2
  • the xdmf file by itself is not enough for anyone to reproduce, since the actual mesh data is stored in mesh1D.h5
  • prepare a standalone (but complete) script, save it in a python file, and run it in a terminal. That will give you more useful information than running in a notebook and seeing a kernel crash.
3

Thank you kindly for your assistance.
As per your advice, I have organized a functional version of the code.
I am truly grateful for your tremendous help, and I would also like to express my deep respect and appreciation for providing such a useful software.

import gmsh
import math
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_mesh
import numpy as np
import ufl
from mpi4py import MPI
import meshio
from dolfinx.io import gmshio
import pyvista as pv

gmsh.initialize()
gmsh.model.add('DFG 3D')

#define dimension
length_member = 500
width_member = 100
depth_member = 100
diameter_longitudinal_rebar = 8
cover_spacing_longitudinal_rebar = 20
embedded_length_longitudinal_rebar = 500
diameter_transverse_rebar = 8
cover_spacing_transverse_rebar = 80
spacing_transverse_rebar = 200
additional_diameter_of_unbonded_region = 1
clear_spacing_between_longitudinal_and_transverse_rebar = 5
radius_of_transverse_rebar = 3 * diameter_transverse_rebar

#####tag list#####
#volume 1X0
tag_volume_longitudinal_rebar = 110
tag_volume_transverse_rebar = 120
tag_volume_concrete = 130
tag_volume_unbonded_region_1 = 140
tag_volume_unbonded_region_2 = 141

#Surface 1X00
tag_surface_longitudinal_rebar = 1100
tag_surface_transverse_rebar = 1200

#Line 1X000
tag_line_transverse_rebar = 12000

#Point 1X0000
tag_point_count_transverse_rebar = 120000

#volume 2X0
marker_volume_longitudinal_rebar = 210
marker_volume_transverse_rebar = 220
marker_volume_concrete = 230

#Surface 2X00
marker_surface_longitudinal_rebar = 2100
marker_surface_transverse_rebar = 2200

##### model concrete #####
domain_concrete = gmsh.model.occ.addBox(0, 0, 0, length_member, width_member, depth_member, tag=tag_volume_concrete)
##### generate concrete #####
gmsh.model.occ.synchronize()

def Calculate_transverse_rebar_positions(length_member, spacing_transverse_rebar):
    if length_member % spacing_transverse_rebar == 0:
        The_number_of_transverse_rebar = math.floor(length_member / spacing_transverse_rebar) - 1
    else:
        The_number_of_transverse_rebar = math.floor(length_member / spacing_transverse_rebar)
    remaining_length = length_member - The_number_of_transverse_rebar * spacing_transverse_rebar
    cover_distance = remaining_length / 2
    positions = [cover_distance + i * spacing_transverse_rebar for i in range(The_number_of_transverse_rebar + 1)]
    return positions

position_transverse_rebar = Calculate_transverse_rebar_positions(length_member, spacing_transverse_rebar)

##### define points for straight line #####
   
point_straight_line_transverse_rebar=[[(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), (cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [width_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), (cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [width_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), depth_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), depth_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)]
                ] # 모따기 아닌부분의 점좌표

##### define start points #####
corner_transverse_rebar = ['Bottom left', 'Bottom right', 'Top right', 'Top left']
list_point_position_individual_transverse_rebar = []
list_point_position_total_transverse_rebar = []


for position in position_transverse_rebar:
    for corner_index in range(len(corner_transverse_rebar)):
        ##### model straight point #####
        gmsh.model.occ.addPoint(position, point_straight_line_transverse_rebar[corner_index][0], point_straight_line_transverse_rebar[corner_index][1], tag=tag_point_count_transverse_rebar)
        gmsh.model.occ.synchronize()
        list_point_position_individual_transverse_rebar.append(tag_point_count_transverse_rebar)
        tag_point_count_transverse_rebar = tag_point_count_transverse_rebar + 1
    list_point_position_total_transverse_rebar.append(list_point_position_individual_transverse_rebar)
    list_point_position_individual_transverse_rebar = []
list_transverse_rebar = []
list_loop_transverse_rebar = []


for position_index in range(len(position_transverse_rebar)):
    ##### model straight line #####
    line1 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][0], list_point_position_total_transverse_rebar[position_index][1], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line1)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line2 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][1], list_point_position_total_transverse_rebar[position_index][2], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line2)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line3 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][2], list_point_position_total_transverse_rebar[position_index][3], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line3)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line4 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][3], list_point_position_total_transverse_rebar[position_index][0], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line4)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    
##### embed straight line #####
for transverse_rebar in list_transverse_rebar:
    gmsh.model.mesh.embed(1, [transverse_rebar[1]], 3, tag_volume_concrete)
    gmsh.model.occ.synchronize()


gmsh.model.occ.removeAllDuplicates()
volume_3D = gmsh.model.getEntities(dim=3)
##### allocate marker and name #####
gmsh.model.addPhysicalGroup(3, [tag_volume_concrete], marker_volume_concrete)
gmsh.model.setPhysicalName(3, marker_volume_concrete, 'volume_concrete')

##### allocate marker and name #####
dimension_list_transverse_rebar = [i[0] for i in list_transverse_rebar]
tag_list_transverse_rebar = [i[1] for i in list_transverse_rebar]
gmsh.model.addPhysicalGroup(1, tag_list_transverse_rebar, marker_volume_transverse_rebar)
gmsh.model.setPhysicalName(1, marker_volume_transverse_rebar, 'line_transverse_rebar')
gmsh.model.occ.synchronize()
   
resolution = min(diameter_longitudinal_rebar, diameter_transverse_rebar) / 2
gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 6)
gmsh.option.setNumber('Mesh.MeshSizeMin', 1 * resolution)
gmsh.option.setNumber('Mesh.MeshSizeMax', 10 * resolution)
gmsh.option.setNumber("Mesh.MshFileVersion", 4.1)  # MSH 파일 버전 (4.1, 2.2 등)
gmsh.option.setNumber("Mesh.Binary", 0)  # ASCII 형식 (Binary 사용 시 1)
gmsh.model.occ.synchronize()
#Generate mesh 3 : Cell, 2 : Facet, 1 : Edge
gmsh.model.mesh.generate(3)
#save model
gmsh_model=gmsh.model
#save msh file
gmsh.write('mesh3D.msh')
gmsh.finalize()
        
mesh_target_for_plot = pv.read('mesh3D.msh')


# gmsh.model.occ.synchronize()
# gmsh.fltk.run()


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[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh
        
# read msh file
mesh1 = meshio.read("mesh3D.msh")

mesh2=create_mesh(mesh1,"tetra")

# Save 3D mesh
meshio.write("mesh3D.xdmf" , mesh2)
   
# Save 1D mesh
mesh3=create_mesh(mesh1,"line")
meshio.write("mesh1D.xdmf", mesh3)

    
with XDMFFile(MPI.COMM_WORLD, "mesh3D.xdmf", "r") as xdmf:
    domain_1 = xdmf.read_mesh(name="Grid")
    #Save cell tag
    cell_marker_1 = xdmf.read_meshtags(domain_1, name="Grid")
xdmf.close()

with XDMFFile(MPI.COMM_WORLD, "mesh1D.xdmf", "r") as xdmf:
    domain_2 = xdmf.read_mesh(name="Grid")
    #Save cell tag
    edge_marker_2 = xdmf.read_meshtags(domain_2, name="Grid")
xdmf.close()
4

What happens if you use the builtin gmshio.model_to_mesh?

5

Thank you, Francesco.
The following error was encountered:

ValueError: Topological dimension cannot be larger than geometric dimension

I noticed that the code runs successfully when the following line is removed:

python

edge_marker_2 = xdmf.read_meshtags(domain_2, name="Grid")

Is there a better way to read the mesh tags for lines?

6

@Henrik_Nicolay_Finsb added the capability to get 1D mesh tags out of a 3D mesh in Extract edge tags and vertex tags from gmsh and return a NamedTuple by finsberg · Pull Request #3547 · FEniCS/dolfinx · GitHub

This commit is not in any release: if you build dolfinx from source and can update to the main branch, please do. If you can’t update or do not build dolfinx from source, a quick-and-dirty workaround that I have verified works on dolfinx 0.9.0 is to:

  1. copy the code from dolfinx/python/dolfinx/io/gmshio.py at 2f4a497e7a3ea03f4ed638eccdafbc31ac1336fe · finsberg/dolfinx · GitHub
  2. save it in a file, e.g. called updated_gmshio.py
  3. remove AdjacencyList from lines 22 and and 455
import gmsh
from mpi4py import MPI
#from dolfinx.io import gmshio  # if you can update to dolfinx main
import updated_gmshio as gmshio  # if you can't update to dolfinx main
import math

gmsh.initialize()
gmsh.model.add('DFG 3D')

#define dimension
length_member = 500
width_member = 100
depth_member = 100
diameter_longitudinal_rebar = 8
cover_spacing_longitudinal_rebar = 20
embedded_length_longitudinal_rebar = 500
diameter_transverse_rebar = 8
cover_spacing_transverse_rebar = 80
spacing_transverse_rebar = 200
additional_diameter_of_unbonded_region = 1
clear_spacing_between_longitudinal_and_transverse_rebar = 5
radius_of_transverse_rebar = 3 * diameter_transverse_rebar

#####tag list#####
#volume 1X0
tag_volume_longitudinal_rebar = 110
tag_volume_transverse_rebar = 120
tag_volume_concrete = 130
tag_volume_unbonded_region_1 = 140
tag_volume_unbonded_region_2 = 141

#Surface 1X00
tag_surface_longitudinal_rebar = 1100
tag_surface_transverse_rebar = 1200

#Line 1X000
tag_line_transverse_rebar = 12000

#Point 1X0000
tag_point_count_transverse_rebar = 120000

#volume 2X0
marker_volume_longitudinal_rebar = 210
marker_volume_transverse_rebar = 220
marker_volume_concrete = 230

#Surface 2X00
marker_surface_longitudinal_rebar = 2100
marker_surface_transverse_rebar = 2200

##### model concrete #####
domain_concrete = gmsh.model.occ.addBox(0, 0, 0, length_member, width_member, depth_member, tag=tag_volume_concrete)
##### generate concrete #####
gmsh.model.occ.synchronize()

def Calculate_transverse_rebar_positions(length_member, spacing_transverse_rebar):
    if length_member % spacing_transverse_rebar == 0:
        The_number_of_transverse_rebar = math.floor(length_member / spacing_transverse_rebar) - 1
    else:
        The_number_of_transverse_rebar = math.floor(length_member / spacing_transverse_rebar)
    remaining_length = length_member - The_number_of_transverse_rebar * spacing_transverse_rebar
    cover_distance = remaining_length / 2
    positions = [cover_distance + i * spacing_transverse_rebar for i in range(The_number_of_transverse_rebar + 1)]
    return positions

position_transverse_rebar = Calculate_transverse_rebar_positions(length_member, spacing_transverse_rebar)

##### define points for straight line #####
   
point_straight_line_transverse_rebar=[[(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), (cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [width_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), (cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [width_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), depth_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)],
                [(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar), depth_member-(cover_spacing_transverse_rebar-diameter_transverse_rebar/2-clear_spacing_between_longitudinal_and_transverse_rebar)]
                ] # 모따기 아닌부분의 점좌표

##### define start points #####
corner_transverse_rebar = ['Bottom left', 'Bottom right', 'Top right', 'Top left']
list_point_position_individual_transverse_rebar = []
list_point_position_total_transverse_rebar = []


for position in position_transverse_rebar:
    for corner_index in range(len(corner_transverse_rebar)):
        ##### model straight point #####
        gmsh.model.occ.addPoint(position, point_straight_line_transverse_rebar[corner_index][0], point_straight_line_transverse_rebar[corner_index][1], tag=tag_point_count_transverse_rebar)
        gmsh.model.occ.synchronize()
        list_point_position_individual_transverse_rebar.append(tag_point_count_transverse_rebar)
        tag_point_count_transverse_rebar = tag_point_count_transverse_rebar + 1
    list_point_position_total_transverse_rebar.append(list_point_position_individual_transverse_rebar)
    list_point_position_individual_transverse_rebar = []
list_transverse_rebar = []
list_loop_transverse_rebar = []


for position_index in range(len(position_transverse_rebar)):
    ##### model straight line #####
    line1 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][0], list_point_position_total_transverse_rebar[position_index][1], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line1)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line2 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][1], list_point_position_total_transverse_rebar[position_index][2], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line2)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line3 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][2], list_point_position_total_transverse_rebar[position_index][3], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line3)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    line4 = gmsh.model.occ.addLine(list_point_position_total_transverse_rebar[position_index][3], list_point_position_total_transverse_rebar[position_index][0], tag=tag_line_transverse_rebar)
    gmsh.model.occ.synchronize()
    transverse_rebar = (1, line4)
    list_transverse_rebar.append(transverse_rebar)
    tag_line_transverse_rebar = tag_line_transverse_rebar + 1
    
##### embed straight line #####
for transverse_rebar in list_transverse_rebar:
    gmsh.model.mesh.embed(1, [transverse_rebar[1]], 3, tag_volume_concrete)
    gmsh.model.occ.synchronize()


gmsh.model.occ.removeAllDuplicates()
volume_3D = gmsh.model.getEntities(dim=3)
##### allocate marker and name #####
gmsh.model.addPhysicalGroup(3, [tag_volume_concrete], marker_volume_concrete)
gmsh.model.setPhysicalName(3, marker_volume_concrete, 'volume_concrete')

gmsh.model.addPhysicalGroup(2, [1], 1)

##### allocate marker and name #####
dimension_list_transverse_rebar = [i[0] for i in list_transverse_rebar]
tag_list_transverse_rebar = [i[1] for i in list_transverse_rebar]
gmsh.model.addPhysicalGroup(1, tag_list_transverse_rebar, marker_volume_transverse_rebar)
gmsh.model.setPhysicalName(1, marker_volume_transverse_rebar, 'line_transverse_rebar')
gmsh.model.occ.synchronize()
   
resolution = min(diameter_longitudinal_rebar, diameter_transverse_rebar) / 2
gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 6)
gmsh.option.setNumber('Mesh.MeshSizeMin', 1 * resolution)
gmsh.option.setNumber('Mesh.MeshSizeMax', 10 * resolution)
gmsh.option.setNumber("Mesh.MshFileVersion", 4.1)  # MSH 파일 버전 (4.1, 2.2 등)
gmsh.option.setNumber("Mesh.Binary", 0)  # ASCII 형식 (Binary 사용 시 1)
gmsh.model.occ.synchronize()
#Generate mesh 3 : Cell, 2 : Facet, 1 : Edge
gmsh.model.mesh.generate(3)
mesh, *tags, names = gmshio.model_to_mesh(
    gmsh.model, comm=MPI.COMM_WORLD, rank=0, gdim=3)
for i in range(len(tags)):
    print(3 - i, tags[i])

Note that I had to add a useless (for you) facet marker with

gmsh.model.addPhysicalGroup(2, [1], 1)

because of [BUG]: #3547 does not handle missing tags for intermediate dimensions · Issue #3594 · FEniCS/dolfinx · GitHub . Subscribe to the discussion there for more updates, or even better contribute there with a minimal example or a test. Nonetheless, the two workarounds reported here should be enough for you to correctly get your cell and line tags in the meantime.

7

Thank you so much for your kind assistance.
I believe that FEM will be one of the most advanced technologies alongside the development of quantum computing and artificial intelligence. In this context, having FEniCSx as an open-source tool is truly a blessing for modern researchers. Once again, I sincerely appreciate your help.