Writing CDB file with SOLID187 elements

[Telos4]

🔍 Before submitting the issue

• [X] I have searched among the existing issues
• [X] I am using a Python virtual environment

🐞 Description of the bug

I've encountered an issue when trying to write a CDB file containing SOLID187 elements using the ansys.mapdl.reader.archive.save_as_archive() function. Apparently, when writing a pyvista mesh of these elements, the mid-nodes get discarded.

I'm creating a tetrahedral mesh consisting of SOLID187 elements using pymapdl:

The elements are 10-node elements with mid-nodes in between the corner elements:

APDL element data (one-indexed):
 mapdl.mesh.elem[0]=array([  1,   1,   1,   1,   0,   0,  14,   0,   1,   0, 169, 167, 101,
       171, 157, 172, 173, 174, 175, 176], dtype=int32)

The last 10 entries of this vector are the node numbers that make up a single SOLID187 element. When saving this mesh as CDB file from APDL using the CDWRITE command the resulting element entries in the file also have 10 unique nodes:

...
EBLOCK,19,SOLID,       100,       100
(19i10)
         1         1         1         1         0         0         0         0        10         0         1       169       167       101       171       157       172       173       174
       175       176
...

Next, I'm converting converting the APDL mesh to a pyvista / VTK mesh (in my case I want to do some modifications / label certain parts of the mesh using python). In pyvista, each element still is made up of 10 unique nodes (4 nodes at the corners, 6 mid-nodes) as confirmed by the inspection of the node data:

VTK element type:
 g.cell[0].type=<CellType.TETRA: 10>
VTK element data (zero-indexed):
 g.cells[1:11]=array([168, 166, 100, 170, 156, 171, 172, 173, 174, 175])

Here's a visualization:

However, if I now save the pyvista mesh as a CDB file using ansys.mapdl.reader.archive.save_as_archive() the resulting elements only contain 4 unique nodes while the rest are duplicates of other nodes (similiar to if they were SOLID185 elements in the tedrahedral configuration):

...
EBLOCK,19,SOLID,       100,       100
(19i8)
       1       1       1       1       0       0       0       0      10       0       1     169     167     101     101     171     171     171     171
...

Is there a way to output a CDB file for SOLID187 elements with the correct nodes (including the mid-nodes)?

📝 Steps to reproduce

Full code to reproduce the above results:

import pyvista as pv

from ansys.mapdl.core import launch_mapdl
from ansys.mapdl import reader as pymapdl_reader

mapdl = launch_mapdl()
print(mapdl.directory)
mapdl.prep7()

# create a tetrahedral mesh with SOLID187 elements
L = 1
mapdl.block(x1=0, x2=L, y1=0, y2=L, z1=0, z2=L)
mapdl.et(1, "SOLID187")
mapdl.vmesh('all')
print(mapdl.mesh)
mapdl.eplot(savefig='ansys_mesh.png', off_screen=True)

# -> elements with 10 unique nodes per element (including mid-nodes)
print(f"APDL element data (one-indexed):\n {mapdl.mesh.elem[0]=}")

# write CDB file
mapdl.cdwrite('DB', 'apdl_archive.cdb')

# export APDL mesh to VTK
grid = mapdl.mesh.grid

# plot mesh including mid-nodes
mb = pv.MultiBlock([grid, pv.PolyData(grid.points)])
mb.plot(off_screen=True, show_edges=True, screenshot='vtk_mesh.png')

# data of first element
# VTK element type = 10 -> tetrahedral elements with 10 nodes (including mid-nodes)
print(f"VTK element type:\n {grid.cell[0].type=}")
print(f"VTK element data (zero-indexed):\n {grid.cells[1:11]=}")

# output VTK data as CDB archive
pymapdl_reader.save_as_archive('vtk_archive.cdb', grid)

mapdl.exit()

💻 Which operating system are you using?

Linux

📀 Which ANSYS version are you using?

2023 R1

🐍 Which Python version are you using?

3.10

📦 Installed packages

ansys-api-mapdl==0.5.1
ansys-api-platform-instancemanagement==1.0.0b3
ansys-mapdl-core==0.63.4
ansys-mapdl-reader==0.52.13
ansys-platform-instancemanagement==1.1.1
appdirs==1.4.4
certifi==2023.5.7
charset-normalizer==3.1.0
contourpy==1.0.7
cycler==0.11.0
fonttools==4.39.4
geomdl==5.3.1
googleapis-common-protos==1.59.0
grpcio==1.54.2
idna==3.4
imageio==2.29.0
importlib-metadata==6.6.0
kiwisolver==1.4.4
matplotlib==3.7.1
numpy==1.24.3
packaging==23.1
pexpect==4.8.0
Pillow==9.5.0
platformdirs==3.5.1
pooch==1.7.0
protobuf==3.20.3
protoc-gen-swagger==0.1.0
ptyprocess==0.7.0
pyiges==0.2.1
pyparsing==3.0.9
python-dateutil==2.8.2
pyvista==0.38.6
requests==2.31.0
scipy==1.10.1
scooby==0.7.2
six==1.16.0
tqdm==4.65.0
urllib3==2.0.2
vtk==9.2.6
zipp==3.15.0

[akaszynski]

Hi @Telos4, The root of the issue is actually within https://github.com/ansys/pymapdl, where we force all elements to be linear in:

https://github.com/ansys/pymapdl/blob/c1eb1e556b9e16923c422d39149ee934fc40c140/src/ansys/mapdl/core/mesh_grpc.py#L720

where we force all cell types to be their linear counterparts with force_linear=True, which I think was done for plotting. The side effect of this is that when writing the output as an archive file it will be written incorrectly.

In the meantime, you can override this behavior with:

import pyvista as pv

...

# export APDL mesh to VTK
grid = mapdl.mesh.grid

# convert any linear tetrahedron to quadratic tetrahedron
grid.celltypes[grid.celltypes == pv.CellType.TETRA] = pv.CellType.QUADRATIC_TETRA
...

@germa89, I recommend we set force_linear=False (default). This may affect plotting, so we might consider always storing the right representation as grid and then using linear_copy during plotting.

[germa89]

Thank you @akaszynski !!

[germa89]

I will transfer this issue to PyMAPDL then.

[germa89]

I think we are going to attack this issue later in the year, when we start to work on the MAPDL backend.