Cannot read facetag when converting .bdf file genarted from hypermesh to xdmf

Hello all,

I’m attempting to replicate the example provided in this tutorial by creating a hexahedral mesh domain using Hypermesh. I have also identified and tagged the facets separately to apply boundary conditions.

However, when I try to read my face tags (quad), my code terminates with a PETSC error message. In contrast, reading the cell tags for the hexahedral elements works without any issues.

I have visualized both XDMF files in ParaView, and the tags are correctly identified.

Could someone help me identify the issue?

MWE (verion-0.7.3)


import meshio
import dolfinx
import mpi4py
from mpi4py import MPI
import numpy as np


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("nastran:ref", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells},
                           cell_data={"Region": [cell_data]})    
    return out_mesh

#'3d problem replication of [1](https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html) by creating mesh using hypermesh and then importing to fenicsx using meshio'

filename='./Schreibtisch/nasdran_files/linear_elasticity_test.bdf'

msh = meshio.read(filename)

#print(msh.cell_data_dict)

quad_mesh = create_mesh(msh, "quad", True)
hex_mesh=create_mesh(msh, "hexahedron", True)

meshio.write("hexa.xdmf", hex_mesh)

meshio.write("quad.xdmf", quad_mesh)

filename_domain = 'hexa.xdmf'

filename_facets = 'quad.xdmf'  #quad meshes to apply fixed boundary condition.

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename_domain, "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    cell_tags = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)


with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename_facets, "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    face_tags = xdmf.read_meshtags(mesh, name="Grid")

Error messaage

[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

Mesh data in .bdf format created from Hypermesh
Mesh_file

No, solution. But I located the error in this smaller script

import meshio
import dolfinx
from mpi4py import MPI


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("nastran:ref", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells},
                           cell_data={"Region": [cell_data]})
    return out_mesh


filename = 'linear_elasticity_test.bdf'

msh = meshio.read(filename)

quad_mesh = create_mesh(msh, "quad", True)

meshio.write("quad.xdmf", quad_mesh)

filename_facets = 'quad.xdmf'  # quad meshes to apply fixed boundary condition.

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename_facets, "r") as xdmf:
    fctmesh = xdmf.read_mesh(name="Grid")
    facet_tags = xdmf.read_meshtags(fctmesh, name="Grid")

to become (on my machine)

Invalid rank, error stack:
internal_Issend(60788): MPI_Issend(buf=0x11d7f46d1, count=1, MPI_BYTE, 1, 1, comm=0x84000001, request=0x11d7ecad4) failed
internal_Issend(60749): Invalid rank has value 1 but must be nonnegative and less than 1
Abort(665453830) on node 0 (rank 0 in comm 464): application called MPI_Abort(comm=0x84000001, 665453830) - process 0

which is obviously raised by the MPI …

Thanks.

But do you have any idea why this error is being raised ? Is there any issue with my mesh generation ?

@dokken @francesco-ballarin, could you please provide some insights into this?

The mesh file can be found in the link mesh file.

I have gone through many posts 1 which deal with almost the same issue but their mesh files are generated from gmesh but in my case, it was from hypermesh.

I have no experience with hypermesh, so it’s difficult for me to comment.

A quick google search shows this post How to save a mesh with *.msh format ? - HyperMesh - Altair Products - Altair Community : can you try what they say (i.e., re-read your mesh in gmsh, and save it in gmsh format)?
The rationale is that we have much better support for reading in gmsh files, and using meshio should now be discouraged.

@francesco-ballarin, Thanks for your reply.

I had tried that already. Although my code works without any error but the node numbers for my quad as well as my hex elements are completely different in comparision to my .bdf file. So the boundary conditions are applied incorrectly.

Below is the result of cells = mesh.get_cells_type(cell_type) using msh file.

[[ 104   76   75 ...  192  191  190]
 [ 110   75   74 ...  191  194  193]
 [ 106  105  111 ...  198  197  196]
 ...
 [ 869  874  886 ... 1014 1026 1028]
 [ 880  888  885 ... 1028 1025 1015]
 [ 888  886  884 ... 1026 1024 1025]]
[[ 4  7  8  3]
 [ 3  8  9  2]
 [ 2  9 10  1]
 [ 0 11 12  6]
 [ 6 12 13  5]
 [ 5 13  7  4]
 [ 7 14 15  8]
 [13 20 14  7]
 [ 8 15 16  9]
 [11 18 19 12]
 [ 9 16 17 10]
 [12 19 20 13]
 [14 21 22 15]
 [20 27 21 14]
 [15 22 23 16]
 [18 25 26 19]
 [16 23 24 17]
 [19 26 27 20]
 [21 28 29 22]
 [27 34 28 21]
 [22 29 30 23]
 [25 32 33 26]
 [23 30 31 24]
 [26 33 34 27]
 [28 35 36 29]
 [34 41 35 28]
 [29 36 37 30]
 [32 39 40 33]
 [30 37 38 31]
 [33 40 41 34]
 [35 42 43 36]
 [41 48 42 35]
 [36 43 44 37]
 [39 46 47 40]
 [37 44 45 38]
 [40 47 48 41]]

and using bdf file.

[[  62   28   27 ...  150  149  148]
 [  68   27   26 ...  149  152  151]
 [  64   63   69 ...  156  155  154]
 ...
 [ 855  860  879 ... 1007 1026 1028]
 [ 870  881  877 ... 1028 1024 1009]
 [ 881  879  875 ... 1026 1022 1024]]
[[  49  273  276   48]
 [  48  276  278   47]
 [  47  278  280   46]
 [   0  285  288   51]
 [  51  288  290   50]
 [  50  290  273   49]
 [ 273  420  423  276]
 [ 290  437  420  273]
 [ 276  423  425  278]
 [ 285  432  435  288]
 [ 278  425  427  280]
 [ 288  435  437  290]
 [ 420  567  570  423]
 [ 437  584  567  420]
 [ 423  570  572  425]
 [ 432  579  582  435]
 [ 425  572  574  427]
 [ 435  582  584  437]
 [ 567  714  717  570]
 [ 584  731  714  567]
 [ 570  717  719  572]
 [ 579  726  729  582]
 [ 572  719  721  574]
 [ 582  729  731  584]
 [ 714  861  864  717]
 [ 731  878  861  714]
 [ 717  864  866  719]
 [ 726  873  876  729]
 [ 719  866  868  721]
 [ 729  876  878  731]
 [ 861 1008 1011  864]
 [ 878 1025 1008  861]
 [ 864 1011 1013  866]
 [ 873 1020 1023  876]
 [ 866 1013 1015  868]
 [ 876 1023 1025  878]]
``

It’s possible that nodes are renumbered in the conversion. It shouldn’t be a problem, provided that the facet markers are modified accordingly, and you only read the mesh generated by gmsh (and not a mixture of the bdf one for the boundary markers and gmsh for the volume)

@francesco-ballarin , You are correct. Even though the nodes are different, the visualizations look exactly the same. I have attached images visualized from data in both MSH and BDF formats.

However, it seems that everything related to face tags is recognized correctly. The issue is that my boundary conditions are not being applied to the quad facet, resulting in a completely different deformed image

The BC is applied to the opposite faces.


import meshio
import dolfinx
import mpi4py
from mpi4py import MPI
import numpy as np
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import petsc4py
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from dolfinx import io


#The function with creates mesh for different element
def create_mesh(mesh, cell_type, prune_z=False):

    #cells --> elements (it creates the each element based on the supplied cell type)
    cells = mesh.get_cells_type(cell_type)
    #if (cell_type=="hexahedron"):
     #   print(cells)

    #cell_data --> This is the function which gets the mesh tags associated with each cells
    cell_data = mesh.get_cell_data("gmsh:geometrical", cell_type)

    #creates the entire mesh 
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells},
                           cell_data={"Region": [cell_data]})
    
    return out_mesh


'3d problem'

filename='./Schreibtisch/nasdran_files/linear_elasticity_test.msh'

msh = meshio.read(filename)
#print(msh.cell_data_dict)

quad_mesh = create_mesh(msh, "quad", True)
hex_mesh=create_mesh(msh, "hexahedron", True)

meshio.write("hexa_from_msh.xdmf", hex_mesh)
meshio.write("quad_from_msh.xdmf", quad_mesh)

filename_domain = 'hexa_from_msh.xdmf'

filename_facets = 'quad_from_msh.xdmf'



# Read it directly in the fenicsx code and start with further analysis of your work


with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename_domain, "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")
    cell_tags = xdmf.read_meshtags(domain, name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)


with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename_facets, "r") as xdmf:
    mesh1 = xdmf.read_mesh(name="Grid")
    face_tags = xdmf.read_meshtags(mesh1, name="Grid")

# Scaled variable

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

#filename='./Schreibtisch/nasdran_files/linear_elasticity_test.msh'
#domain,markers,facets=io.gmshio.read_from_msh(filename,MPI.COMM_WORLD)  


V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
fdim = domain.topology.dim - 1

# clamped boundary
facetClamped = face_tags.indices[face_tags.values==2] # physical group "bottom" with tag 2 (line)
#print(facetClamped)

dofsClamped = dolfinx.fem.locate_dofs_topological(V = V, entity_dim = fdim,  entities=facetClamped)
uClamped = dolfinx.fem.Constant(domain, petsc4py.PETSc.ScalarType((0.0,0.0,0.0)))
bc = dolfinx.fem.dirichletbc(uClamped, dofsClamped, V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)#facetClamped = dolfinx.mesh.locate_entities_boundary(mesh=domain, dim=facetDim, marker=lambda x: np.isclose(x[1], 0))

f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)

Stop using meshio, and use dolfinx.io.gmshio — DOLFINx 0.8.0 documentation instead on the msh file.

@francesco-ballarin , i have tried that too earlier. But it throws a error message.

Info    : Reading './Schreibtisch/nasdran_files/linear_elasticity_test.msh'...
Info    : 2 entities
Info    : 1029 nodes
Info    : 756 elements
Info    : Done reading './Schreibtisch/nasdran_files/linear_elasticity_test.msh'
Traceback (most recent call last):
  File "/home/peri_pr/Schreibtisch/testing.py", line 42, in <module>
    msh = dolfinx.io.gmshio.read_from_msh(filename,MPI.COMM_WORLD)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peri_pr/miniforge3/envs/fenicsx0.7.3-env/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 332, in read_from_msh
    msh = model_to_mesh(gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peri_pr/miniforge3/envs/fenicsx0.7.3-env/lib/python3.12/site-packages/dolfinx/io/gmshio.py", line 232, in model_to_mesh
    cell_id = cell_information[perm_sort[-1]]["id"]
                               ~~~~~~~~~^^^^
IndexError: index -1 is out of bounds for axis 0 with size 0
``

If you don’t have subdomain markers (cell tags in dolfinx), make sure to add them in your hyermesh generator (even if there isn’t an actual subdomain, and all cells get the same marker).

The domain and the subdomain are already provided with marker as specified below.

print(mesh.cell_data_dict)

The .msh file also produces the same set of labels when i read it using meshio

{'nastran:ref': {'hexahedron': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]), 'quad': array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])}}

There’s not much more I can do unless you prepare a link to the msh file.

@francesco-ballarin , Please find link to download mesh files in both bdf as well as in ```msh`` format.

Many thanks for your support.

Add a print(topologies) at the end of extract_topology_and_markers in dolfinx/io/gmshio.py. You’ll see that it prints the empty dictionary. This confirms my suspicion (as stated a few posts earlier) that subdomain and/or boundary markers are not present in the msh file. Please check if this is due to hypermesh, or due to having imported the file in gmsh, possibly on the corresponding support forums. I don’t have hypermesh, so it’s unlikely I can help much further.

1 Like

@francesco-ballarin, thanks for your suggestion.

I’m still wondering if extract_topology_and_markers in dolfinx/io/gmshio.py is creating an empty dictionary. How was I able to visualize the XDMF file in ParaView using threshold filters? I believe these filters are my mesh_tags.

Is this thinking correct?

That visualization isn’t helpful for debugging the issue, since it is based on xdmf files, and not the msh one. Hence, you must understand why the msh file is missing the tags (or physical groups as gmsh would call them).

Great post, learned a lot.

The mesh file has expired, so I cannot access it and try to reproduce the result.

@dokken , Please find the link to download both the mesh files (msh and bdf).