meshio icon indicating copy to clipboard operation
meshio copied to clipboard

[BUG] create abaqus mesh with correct element type not possible anymore

Open reox opened this issue 3 years ago • 1 comments

Describe the bug In previous versions of meshio (for example 4.4.6), it was possible to create a mesh like this:

mesh = meshio.Mesh(nodes, {'C3D8': elements})
meshio.write("some_mesh.inp", mesh, translate_cell_names=False)

However, this fails with a recent version (5.3.2):

KeyError: 'C3D8'

Traceback (most recent calls last):
  File "/home/reox/hexMesh.py", line 150, in myplugin
    mesh = meshio.Mesh(nodes, {'C3D8': elements})
  File "/home/reox/.local/lib/python3.7/site-packages/meshio/_mesh.py", line 150, in __init__
    data if cell_type.startswith("polyhedron") else np.asarray(data),
  File "/home/reox/.local/lib/python3.7/site-packages/meshio/_mesh.py", line 99, in __init__
    self.dim = topological_dimension[cell_type]

The issue is, that if I create the mesh using hexahedron as a type, I can not be sure to get C3D8 elements in the resulting mesh. The reason is this lookup: https://github.com/nschloe/meshio/blob/3510aeaed8cf233363c5b365c8446d88b76b7a16/src/meshio/abaqus/_abaqus.py#L416-L418 The dictionary that is used will choose basically the last type that references to hexahedron, which in my case is C3D8RH: https://github.com/nschloe/meshio/blob/3510aeaed8cf233363c5b365c8446d88b76b7a16/src/meshio/abaqus/_abaqus.py#L100

One possibility would be to add a custom mapping to the abaqus.write function, i.e. something like: mapping=dict(hexahedron='C3D8') or to restore the old state were the element type could be written.

A workaround is to simply overwrite the dict before writing the mesh:

meshio.abaqus._abaqus.meshio_to_abaqus_type['hexahedron'] = 'C3D8'

To Reproduce

import numpy as np
import meshio
nodes = np.array([[0,0,0],[0,0,1],[0,1,1],[0,1,0],[1,0,0],[1,0,1],[1,1,1],[1,1,0],])
elements = np.array([[0,1,2,3,4,5,6,7]])
# old behaviour, working in 4.4.6, now crashing:
mesh = meshio.Mesh(nodes, {'C3D8': elements})  # crash
meshio.write("some_mesh.inp", mesh, translate_cell_names=False)

# new behaviour, not writing C3D8 (or the element type you want) but always C3D8RH
mesh = meshio.Mesh(nodes, {'hexahedron': elements}) 
meshio.write("some_mesh.inp", mesh)  # writes "wrong" element type

# Workaround:
mesh = meshio.Mesh(nodes, {'hexahedron': elements})
meshio.abaqus._abaqus.meshio_to_abaqus_type['hexahedron'] = 'C3D8'
meshio.write("some_mesh.inp", mesh)  # writes C3D8 again

Diagnose I may ask you to cut and paste the output of the following command.

$ pip freeze | grep meshio
meshio==5.3.2

Did I help?

If I was able to resolve your problem, consider sponsoring my work on meshio, or buy me a coffee to say thanks.

reox avatar Mar 04 '22 07:03 reox

Also noticed that for triangle elements: They are written as R3D3 by default, while I would like to have S3 as the default type.

I think there are two parts to this issue here:

  1. choosing a sensible default type
  2. being able to define the element type yourself.

Currently, the mapping of default types is a bit odd:

line            ->  B31H           
line3           ->  B33H           
quad            ->  CAX4P          
quad8           ->  S8R5           
quad9           ->  S9R5           
triangle        ->  R3D3           
triangle6       ->  CPE6           
hexahedron      ->  C3D8RH         
hexahedron20    ->  C3D20RH        
tetra           ->  C3D4           
tetra4          ->  C3D4H          
tetra10         ->  C3D10MH        
wedge           ->  C3D6           
wedge15         ->  C3D15          

For many element types, special element types are used, for example C3D20RH. The abaqus_to_meshio_type dict is already in somewhat a sensible order, with the most simplest element type first. Thus, a quick fix for the first issue would be to change https://github.com/nschloe/meshio/blob/3510aeaed8cf233363c5b365c8446d88b76b7a16/src/meshio/abaqus/_abaqus.py#L100 to:

meshio_to_abaqus_type = {k: v for k, v in reversed(abaqus_to_meshio_type.items())}

reox avatar Jul 25 '23 11:07 reox