Skip to content

Instantly share code, notes, and snippets.

@QiangF
Created July 18, 2026 14:00
Show Gist options
  • Select an option

  • Save QiangF/68d1cd30c0a9aafa2d2d85b1512eceff to your computer and use it in GitHub Desktop.

Select an option

Save QiangF/68d1cd30c0a9aafa2d2d85b1512eceff to your computer and use it in GitHub Desktop.

Based on my thorough study of the extracted source code from both animtod3plot and openradioss-to-vtkhdf, I now have a comprehensive understanding of the Radioss animation file format. Let me summarize what I've learned and provide you with a complete PyVista reader implementation.

Radioss Animation File Format Summary

The Radioss animation file (.anim) is a big-endian binary format containing:

File Header (Flags)

  • magic - file format version
  • time - simulation time for this state
  • Various text fields (time_text, mod_anim_text, radioss_run_text)
  • Flag A0-A9: control presence of mass data, numbering, 3D geometry, 1D geometry, hierarchy, time history, skew tensors, SPH, etc.

2D Geometry (Shells)

  • nbNodes - total nodes
  • nbFacets - 4-node shell elements
  • nbParts - number of parts
  • nbFunc - nodal scalar values
  • nbEFunc - element scalar values
  • nbVect - vectors
  • nbTens - tensors
  • Arrays: coorA (3×nbNodes), connectA (4×nbFacets), delEltA, funcA, eFuncA, vectValA, tensValA

3D Geometry (Solids)

  • nbElts3D - 8-node solid elements
  • connect3DA (8×nbElts3D), eFunc3DA, tensVal3DA (6 components)

1D Geometry (Beams)

  • nbElts1D - 2-node beam elements
  • connect1DA (2×nbElts1D), torsVal1DA (9 components: 3 forces + 3 moments + 3 others)

SPH Elements

  • Point-like particles with connectivity

PyVista Reader Implementation

Here's a complete, standalone Radioss animation file reader for PyVista:

"""
Radioss Animation File Reader for PyVista
==========================================

A standalone reader for Altair Radioss animation files (.anim)
that converts the data into PyVista UnstructuredGrid and PolyData objects.

Based on the animtod3plot RadiossReader and openradioss-to-vtkhdf projects.

Usage:
    import pyvista as pv
    from radioss_reader import read_radioss_anim
    
    # Read a single animation file
    mesh = read_radioss_anim("file.anim")
    mesh.plot()
    
    # Or register as a PyVista reader
    import radioss_reader
    radioss_reader.register_pyvista_reader()
    
    mesh = pv.read("file.anim")
"""

from __future__ import annotations

import struct
import warnings
from pathlib import Path
from typing import BinaryIO, Dict, List, Optional, Tuple, Union

import numpy as np
import pyvista as pv


class RadiossAnimReader:
    """
    Reader for Radioss Animation Files (.anim)
    
    The Radioss animation format is a big-endian binary format containing
    finite element mesh data (nodes, shells, solids, beams, SPH) with
    associated scalar/vector/tensor results at a single time state.
    """
    
    # Data types - Radioss uses big-endian format
    ITYPE: np.dtype = np.dtype(">i4")   # 32-bit big-endian integer
    ITYPE_S: np.dtype = np.dtype(">i2") # 16-bit big-endian integer
    FTYPE: np.dtype = np.dtype(">f4")  # 32-bit big-endian float
    BTYPE: np.dtype = np.dtype("bool")
    
    WORDSIZE: int = 4
    WORDSIZE_S: int = 2
    WORDSIZE_B: int = 1
    
    def __init__(self, filepath: Union[str, Path, bytes, BinaryIO]):
        """
        Initialize reader and parse the animation file.
        
        Parameters
        ----------
        filepath : str, Path, bytes, or file-like
            Path to the .anim file or raw binary data
        """
        self.filepath = filepath
        self.raw_header: Dict = {}
        self.raw_arrays: Dict = {}
        self.arrays: Dict = {}
        
        # Parse the file
        self._load_file(filepath)
        self._unpack_arrays()
    
    def _read_buffer(self, source: Union[str, Path, bytes, BinaryIO]) -> bytes:
        """Read binary data from various sources."""
        if isinstance(source, (str, Path)):
            with open(source, 'rb') as f:
                return f.read()
        elif isinstance(source, bytes):
            return source
        elif hasattr(source, 'read'):
            return source.read()
        else:
            raise TypeError(f"Cannot read from type {type(source)}")
    
    def _read_single(self, data: bytes, pos: int, dtype: np.dtype) -> Tuple[Union[int, float], int]:
        """Read a single value and advance position."""
        size = dtype.itemsize
        value = np.frombuffer(data[pos:pos+size], dtype=dtype, count=1)[0]
        return value, pos + size
    
    def _read_string(self, data: bytes, pos: int, length: int) -> Tuple[str, int]:
        """Read a fixed-length string and advance position."""
        raw = data[pos:pos+length]
        # Decode and strip null bytes
        text = raw.decode('latin-1', errors='replace').replace('\x00', '').replace('\x01', '').strip()
        return text, pos + length
    
    def _read_array(self, data: bytes, pos: int, n_items: int, width: int, 
                    dtype: np.dtype) -> Tuple[np.ndarray, int]:
        """Read a numpy array and advance position."""
        if n_items == 0:
            return np.array([]), pos
        
        total_items = n_items * width
        size = total_items * dtype.itemsize
        arr = np.frombuffer(data[pos:pos+size], dtype=dtype, count=total_items)
        
        if width > 1:
            arr = arr.reshape((n_items, width))
        
        return arr, pos + size
    
    def _read_string_array(self, data: bytes, pos: int, n_items: int, 
                         str_len: int) -> Tuple[List[str], int]:
        """Read an array of fixed-length strings."""
        if n_items == 0:
            return [], pos
        
        strings = []
        for _ in range(n_items):
            s, pos = self._read_string(data, pos, str_len)
            strings.append(s)
        
        return strings, pos
    
    def _load_file(self, filepath: Union[str, Path, bytes, BinaryIO]) -> None:
        """Parse the binary animation file."""
        data = self._read_buffer(filepath)
        pos = 0
        
        # ==================== HEADER ====================
        # Magic number (file version)
        self.raw_header["magic"], pos = self._read_single(data, pos, self.ITYPE)
        
        # Time value
        self.raw_header["time"], pos = self._read_single(data, pos, self.FTYPE)
        
        # Text fields (81 chars each)
        self.raw_header["time_text"], pos = self._read_string(data, pos, 81)
        self.raw_header["mod_anim_text"], pos = self._read_string(data, pos, 81)
        self.raw_header["radioss_run_text"], pos = self._read_string(data, pos, 81)
        
        # Flags A0-A9
        for i in range(10):
            self.raw_header[f"flag_a_{i}"], pos = self._read_single(data, pos, self.ITYPE)
        
        # ==================== 2D GEOMETRY (SHELLS) ====================
        self.raw_header["nbNodes"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbFacets"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbParts"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbFunc"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbEFunc"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbVect"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbTens"], pos = self._read_single(data, pos, self.ITYPE)
        self.raw_header["nbSkew"], pos = self._read_single(data, pos, self.ITYPE)
        
        # Skew values (if any)
        self.raw_arrays["skewValA"], pos = self._read_array(
            data, pos, self.raw_header["nbSkew"], 6, self.ITYPE_S)
        
        # Node coordinates (3 × nbNodes)
        self.raw_arrays["coorA"], pos = self._read_array(
            data, pos, self.raw_header["nbNodes"], 3, self.FTYPE)
        
        # Shell connectivity (4 × nbFacets) - 1-based indexing in file
        self.raw_arrays["connectA"], pos = self._read_array(
            data, pos, self.raw_header["nbFacets"], 4, self.ITYPE)
        
        # Deleted element flags
        self.raw_arrays["delEltA"], pos = self._read_array(
            data, pos, self.raw_header["nbFacets"], 1, self.BTYPE)
        
        # Part definitions (cumulative counts)
        self.raw_arrays["defPartA"], pos = self._read_array(
            data, pos, self.raw_header["nbParts"], 1, self.ITYPE)
        
        # Part names
        self.raw_arrays["pTextA"], pos = self._read_string_array(
            data, pos, self.raw_header["nbParts"], 50)
        
        # Normal vectors (as int16, scaled)
        self.raw_arrays["normFloatA"], pos = self._read_array(
            data, pos, self.raw_header["nbNodes"], 3, self.ITYPE_S)
        
        # Scalar function names (nbFunc + nbEFunc)
        n_fnames = self.raw_header["nbFunc"] + self.raw_header["nbEFunc"]
        self.raw_arrays["fTextA"], pos = self._read_string_array(
            data, pos, n_fnames, 81)
        
        # Nodal scalar values (nbFunc × nbNodes)
        self.raw_arrays["funcA"], pos = self._read_array(
            data, pos, self.raw_header["nbFunc"] * self.raw_header["nbNodes"], 1, self.FTYPE)
        
        # Element scalar values (nbEFunc × nbFacets)
        self.raw_arrays["eFuncA"], pos = self._read_array(
            data, pos, self.raw_header["nbEFunc"] * self.raw_header["nbFacets"], 1, self.FTYPE)
        
        # Vector names
        self.raw_arrays["vTextA"], pos = self._read_string_array(
            data, pos, self.raw_header["nbVect"], 81)
        
        # Vector values (nbVect × nbNodes × 3)
        self.raw_arrays["vectValA"], pos = self._read_array(
            data, pos, self.raw_header["nbVect"] * self.raw_header["nbNodes"], 3, self.FTYPE)
        
        # Tensor names
        self.raw_arrays["tTextA"], pos = self._read_string_array(
            data, pos, self.raw_header["nbTens"], 81)
        
        # Tensor values (nbTens × nbFacets × 3) - note: 3 components in 2D
        self.raw_arrays["tensValA"], pos = self._read_array(
            data, pos, self.raw_header["nbTens"] * self.raw_header["nbFacets"], 3, self.FTYPE)
        
        # Mass data (if flag_a_0 == 1)
        if self.raw_header["flag_a_0"] == 1:
            self.raw_arrays["eMassA"], pos = self._read_array(
                data, pos, self.raw_header["nbFacets"], 1, self.FTYPE)
            self.raw_arrays["nMassA"], pos = self._read_array(
                data, pos, self.raw_header["nbNodes"], 1, self.FTYPE)
        
        # Node/element numbering (if flag_a_1)
        if self.raw_header["flag_a_1"]:
            self.raw_arrays["nodNumA"], pos = self._read_array(
                data, pos, self.raw_header["nbNodes"], 1, self.ITYPE)
            self.raw_arrays["elNumA"], pos = self._read_array(
                data, pos, self.raw_header["nbFacets"], 1, self.ITYPE)
        
        # Hierarchy data (if flag_a_4)
        if self.raw_header["flag_a_4"]:
            self.raw_arrays["part2subset2DA"], pos = self._read_array(
                data, pos, self.raw_header["nbParts"], 1, self.ITYPE)
            self.raw_arrays["partMaterial2DA"], pos = self._read_array(
                data, pos, self.raw_header["nbParts"], 1, self.ITYPE)
            self.raw_arrays["partProperties2DA"], pos = self._read_array(
                data, pos, self.raw_header["nbParts"], 1, self.ITYPE)
        
        # ==================== 3D GEOMETRY (SOLIDS) ====================
        if self.raw_header["flag_a_2"]:
            self.raw_header["nbElts3D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbParts3D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbEFunc3D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbTens3D"], pos = self._read_single(data, pos, self.ITYPE)
            
            if self.raw_header["nbElts3D"] > 0:
                self.raw_arrays["connect3DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts3D"], 8, self.ITYPE)
                self.raw_arrays["delElt3DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts3D"], 1, self.BTYPE)
                self.raw_arrays["defPart3DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbParts3D"], 1, self.ITYPE)
                self.raw_arrays["pText3DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbParts3D"], 50)
                self.raw_arrays["fText3DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbEFunc3D"], 81)
                self.raw_arrays["eFunc3DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbEFunc3D"] * self.raw_header["nbElts3D"], 1, self.FTYPE)
                self.raw_arrays["tText3DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbTens3D"], 81)
                self.raw_arrays["tensVal3DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts3D"] * self.raw_header["nbTens3D"], 6, self.FTYPE)
                
                if self.raw_header["flag_a_0"] == 1:
                    self.raw_arrays["eMass3DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbElts3D"], 1, self.FTYPE)
                if self.raw_header["flag_a_1"] == 1:
                    self.raw_arrays["elNum3DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbElts3D"], 1, self.ITYPE)
                if self.raw_header["flag_a_4"]:
                    self.raw_arrays["part2subset3DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts3D"], 1, self.ITYPE)
                    self.raw_arrays["partMaterial3DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts3D"], 1, self.ITYPE)
                    self.raw_arrays["partProperties3DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts3D"], 1, self.ITYPE)
        else:
            self.raw_header["nbElts3D"] = 0
            self.raw_header["nbParts3D"] = 0
            self.raw_header["nbEFunc3D"] = 0
            self.raw_header["nbTens3D"] = 0
        
        # ==================== 1D GEOMETRY (BEAMS) ====================
        if self.raw_header["flag_a_3"]:
            self.raw_header["nbElts1D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbParts1D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbEFunc1D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbTors1D"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["isSkew1D"], pos = self._read_single(data, pos, self.ITYPE)
            
            if self.raw_header["nbElts1D"] > 0:
                self.raw_arrays["connect1DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts1D"], 2, self.ITYPE)
                self.raw_arrays["delElt1DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts1D"], 1, self.BTYPE)
                self.raw_arrays["defPart1DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbParts1D"], 1, self.ITYPE)
                self.raw_arrays["pText1DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbParts1D"], 50)
                self.raw_arrays["fText1DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbEFunc1D"], 81)
                self.raw_arrays["eFunc1DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbEFunc1D"] * self.raw_header["nbElts1D"], 1, self.FTYPE)
                self.raw_arrays["tText1DA"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbTors1D"], 81)
                self.raw_arrays["torsVal1DA"], pos = self._read_array(
                    data, pos, self.raw_header["nbElts1D"] * self.raw_header["nbTors1D"], 9, self.FTYPE)
                
                if self.raw_header["isSkew1D"]:
                    self.raw_arrays["elt2Skew1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbElts1D"], 1, self.ITYPE)
                if self.raw_header["flag_a_0"] == 1:
                    self.raw_arrays["eMass1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbElts1D"], 1, self.FTYPE)
                if self.raw_header["flag_a_1"] == 1:
                    self.raw_arrays["elNum1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbElts1D"], 1, self.ITYPE)
                if self.raw_header["flag_a_4"]:
                    self.raw_arrays["part2subset1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts1D"], 1, self.ITYPE)
                    self.raw_arrays["partMaterial1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts1D"], 1, self.ITYPE)
                    self.raw_arrays["partProperties1DA"], pos = self._read_array(
                        data, pos, self.raw_header["nbParts1D"], 1, self.ITYPE)
        else:
            self.raw_header["nbElts1D"] = 0
            self.raw_header["nbParts1D"] = 0
            self.raw_header["nbEFunc1D"] = 0
            self.raw_header["nbTors1D"] = 0
            self.raw_header["isSkew1D"] = 0
        
        # ==================== HIERARCHY ====================
        if self.raw_header["flag_a_4"]:
            self.raw_header["nbSubsets"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_arrays["nbSubsets"] = {}
            
            for _ in range(self.raw_header["nbSubsets"]):
                subset_name, pos = self._read_string(data, pos, 50)
                self.raw_arrays["nbSubsets"][subset_name] = {}
                
                num_parent, pos = self._read_single(data, pos, self.ITYPE)
                self.raw_arrays["nbSubsets"][subset_name]["numParent"] = num_parent
                
                nb_sons, pos = self._read_single(data, pos, self.ITYPE)
                self.raw_arrays["nbSubsets"][subset_name]["nbSubsetSon"] = nb_sons
                
                if nb_sons > 0:
                    sons, pos = self._read_array(data, pos, nb_sons, 1, self.ITYPE)
                    self.raw_arrays["nbSubsets"][subset_name]["subsetSonA"] = sons
                
                nb_subpart2d, pos = self._read_single(data, pos, self.ITYPE)
                self.raw_arrays["nbSubsets"][subset_name]["nbSubPart2D"] = nb_subpart2d
                if nb_subpart2d > 0:
                    parts, pos = self._read_array(data, pos, nb_subpart2d, 1, self.ITYPE)
                    self.raw_arrays["nbSubsets"][subset_name]["subPart2DA"] = parts
                
                nb_subpart3d, pos = self._read_single(data, pos, self.ITYPE)
                self.raw_arrays["nbSubsets"][subset_name]["nbSubPart3D"] = nb_subpart3d
                if nb_subpart3d > 0:
                    parts, pos = self._read_array(data, pos, nb_subpart3d, 1, self.ITYPE)
                    self.raw_arrays["nbSubsets"][subset_name]["subPart3DA"] = parts
                
                nb_subpart1d, pos = self._read_single(data, pos, self.ITYPE)
                self.raw_arrays["nbSubsets"][subset_name]["nbSubPart1D"] = nb_subpart1d
                if nb_subpart1d > 0:
                    parts, pos = self._read_array(data, pos, nb_subpart1d, 1, self.ITYPE)
                    self.raw_arrays["nbSubsets"][subset_name]["subPart1DA"] = parts
            
            self.raw_header["nbMaterials"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbProperties"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_arrays["materialTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbMaterials"], 50)
            self.raw_arrays["materialTypeA"], pos = self._read_array(
                data, pos, self.raw_header["nbMaterials"], 1, self.ITYPE)
            self.raw_arrays["propertiesTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbProperties"], 50)
            self.raw_arrays["propertiesTypeA"], pos = self._read_array(
                data, pos, self.raw_header["nbProperties"], 1, self.ITYPE)
        else:
            self.raw_header["nbSubsets"] = 0
            self.raw_header["nbMaterials"] = 0
            self.raw_header["nbProperties"] = 0
        
        # ==================== TIME HISTORY NODES/ELEMENTS ====================
        if self.raw_header["flag_a_5"]:
            self.raw_header["nbNodesTH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbElts2DTH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbElts3DTH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbElts1DTH"], pos = self._read_single(data, pos, self.ITYPE)
            
            self.raw_arrays["nodes2THA"], pos = self._read_array(
                data, pos, self.raw_header["nbNodesTH"], 1, self.ITYPE)
            self.raw_arrays["n2thTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbNodesTH"], 50)
            self.raw_arrays["elt2DTHA"], pos = self._read_array(
                data, pos, self.raw_header["nbElts2DTH"], 1, self.ITYPE)
            self.raw_arrays["elt2DthTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbElts2DTH"], 50)
            self.raw_arrays["elt3DTHA"], pos = self._read_array(
                data, pos, self.raw_header["nbElts3DTH"], 1, self.ITYPE)
            self.raw_arrays["elt3DthTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbElts3DTH"], 50)
            self.raw_arrays["elt1DTHA"], pos = self._read_array(
                data, pos, self.raw_header["nbElts1DTH"], 1, self.ITYPE)
            self.raw_arrays["elt1DthTextA"], pos = self._read_string_array(
                data, pos, self.raw_header["nbElts1DTH"], 50)
        else:
            self.raw_header["nbNodesTH"] = 0
            self.raw_header["nbElts2DTH"] = 0
            self.raw_header["nbElts3DTH"] = 0
            self.raw_header["nbElts1DTH"] = 0
        
        # ==================== SPH ELEMENTS ====================
        if self.raw_header["flag_a_7"]:
            self.raw_header["nbEltsSPH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbPartsSPH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbEFuncSPH"], pos = self._read_single(data, pos, self.ITYPE)
            self.raw_header["nbTensSPH"], pos = self._read_single(data, pos, self.ITYPE)
            
            if self.raw_header["nbEltsSPH"] > 0:
                self.raw_arrays["connecSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"], 1, self.ITYPE)
                self.raw_arrays["delEltSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"], 1, self.BTYPE)
            
            if self.raw_header["nbPartsSPH"] > 0:
                self.raw_arrays["defPartSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbPartsSPH"], 1, self.ITYPE)
                self.raw_arrays["pTextSPH"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbPartsSPH"], 50)
            
            if self.raw_header["nbEFuncSPH"] > 0:
                self.raw_arrays["scalTextSPH"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbEFuncSPH"], 81)
                self.raw_arrays["eFuncSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"] * self.raw_header["nbEFuncSPH"], 1, self.ITYPE)
            
            if self.raw_header["nbTensSPH"] > 0:
                self.raw_arrays["tensTextSPH"], pos = self._read_string_array(
                    data, pos, self.raw_header["nbTensSPH"], 81)
                self.raw_arrays["tensValSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"] * self.raw_header["nbTensSPH"], 6, self.FTYPE)
            
            if self.raw_header["flag_a_0"] == 1:
                self.raw_arrays["eMassSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"], 1, self.FTYPE)
            if self.raw_header["flag_a_1"] == 1:
                self.raw_arrays["nodNumSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbEltsSPH"], 1, self.ITYPE)
            if self.raw_header["flag_a_4"]:
                self.raw_arrays["numParentSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbPartsSPH"], 1, self.ITYPE)
                self.raw_arrays["matPartSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbPartsSPH"], 1, self.ITYPE)
                self.raw_arrays["propPartSPH"], pos = self._read_array(
                    data, pos, self.raw_header["nbPartsSPH"], 1, self.ITYPE)
        else:
            self.raw_header["nbEltsSPH"] = 0
            self.raw_header["nbPartsSPH"] = 0
            self.raw_header["nbEFuncSPH"] = 0
            self.raw_header["nbTensSPH"] = 0
    
    def _unpack_arrays(self) -> None:
        """Unpack raw arrays into user-friendly format."""
        self.arrays = {}
        self.arrays["timesteps"] = self.raw_header["time"]
        
        # Node coordinates
        if "coorA" in self.raw_arrays and len(self.raw_arrays["coorA"]) > 0:
            self.arrays["node_coordinates"] = self.raw_arrays["coorA"]
        
        # Element connectivities (convert from 1-based to 0-based indexing)
        if "connect1DA" in self.raw_arrays and len(self.raw_arrays["connect1DA"]) > 0:
            self.arrays["element_beam_node_indexes"] = self.raw_arrays["connect1DA"] - 1
        
        if "connectA" in self.raw_arrays and len(self.raw_arrays["connectA"]) > 0:
            self.arrays["element_shell_node_indexes"] = self.raw_arrays["connectA"] - 1
        
        if "connect3DA" in self.raw_arrays and len(self.raw_arrays["connect3DA"]) > 0:
            self.arrays["element_solid_node_indexes"] = self.raw_arrays["connect3DA"] - 1
        
        if "connecSPH" in self.raw_arrays and len(self.raw_arrays["connecSPH"]) > 0:
            self.arrays["sph_node_indexes"] = self.raw_arrays["connecSPH"] - 1
        
        # Node/element IDs
        if "nodNumA" in self.raw_arrays:
            self.arrays["node_ids"] = self.raw_arrays["nodNumA"]
        if "elNumA" in self.raw_arrays:
            self.arrays["element_shell_ids"] = self.raw_arrays["elNumA"]
        if "elNum3DA" in self.raw_arrays:
            self.arrays["element_solid_ids"] = self.raw_arrays["elNum3DA"]
        if "elNum1DA" in self.raw_arrays:
            self.arrays["element_beam_ids"] = self.raw_arrays["elNum1DA"]
        if "nodNumSPH" in self.raw_arrays:
            self.arrays["element_sph_ids"] = self.raw_arrays["nodNumSPH"]
        
        # Unpack nodal scalars
        if "fTextA" in self.raw_arrays and self.raw_header["nbFunc"] > 0:
            for ifun in range(self.raw_header["nbFunc"]):
                name = f"node_{self.raw_arrays['fTextA'][ifun].lower().replace(' ', '_').strip()}"
                start = ifun * self.raw_header["nbNodes"]
                end = (ifun + 1) * self.raw_header["nbNodes"]
                self.arrays[name] = self.raw_arrays["funcA"][start:end]
        
        # Unpack nodal vectors
        if "vTextA" in self.raw_arrays and self.raw_header["nbVect"] > 0:
            for ivect in range(self.raw_header["nbVect"]):
                name = f"node_{self.raw_arrays['vTextA'][ivect].lower().replace(' ', '_').strip()}"
                start = ivect * self.raw_header["nbNodes"]
                end = (ivect + 1) * self.raw_header["nbNodes"]
                self.arrays[name] = self.raw_arrays["vectValA"][start:end]
        
        # Unpack element scalars (2D)
        if "fTextA" in self.raw_arrays and self.raw_header["nbEFunc"] > 0:
            for ifun in range(self.raw_header["nbEFunc"]):
                name = f"element_shell_{self.raw_arrays['fTextA'][self.raw_header['nbFunc'] + ifun].lower().replace(' ', '_').strip()}"
                start = ifun * self.raw_header["nbFacets"]
                end = (ifun + 1) * self.raw_header["nbFacets"]
                self.arrays[name] = self.raw_arrays["eFuncA"][start:end]
        
        # Unpack element tensors (2D)
        if "tTextA" in self.raw_arrays and self.raw_header["nbTens"] > 0:
            for itens in range(self.raw_header["nbTens"]):
                name = f"element_shell_{self.raw_arrays['tTextA'][itens].lower().replace(' ', '_').strip()}"
                start = itens * self.raw_header["nbFacets"]
                end = (itens + 1) * self.raw_header["nbFacets"]
                self.arrays[name] = self.raw_arrays["tensValA"][start:end]
        
        # Unpack 1D element part indexes
        if "defPart1DA" in self.raw_arrays and len(self.raw_arrays["defPart1DA"]) > 0:
            self._unpack_part_indexes("1D", "beam")
        
        # Unpack 2D element part indexes
        if self.raw_header["nbParts"] > 0 and "defPartA" in self.raw_arrays and len(self.raw_arrays["defPartA"]) > 0:
            self._unpack_part_indexes("2D", "shell")
        
        # Unpack 3D element part indexes
        if self.raw_header["nbParts3D"] > 0 and "defPart3DA" in self.raw_arrays and len(self.raw_arrays["defPart3DA"]) > 0:
            self._unpack_part_indexes("3D", "solid")
    
    def _unpack_part_indexes(self, dim: str, elem_type: str) -> None:
        """Unpack cumulative part definitions into per-element arrays."""
        def_part = self.raw_arrays[f"defPart{dim}A"]
        p_text = self.raw_arrays[f"pText{dim}A"]
        
        if len(def_part) == 0:
            return
        
        n_elements = def_part[-1] if len(def_part) > 0 else 0
        if n_elements == 0:
            return
        
        tmp_list_i = np.empty(n_elements, dtype=int)
        tmp_list_n = np.empty(n_elements, dtype=object)
        
        start = 0
        for ipart, end in enumerate(def_part):
            end_val = int(end)
            if start < n_elements:
                actual_end = min(end_val, n_elements)
                tmp_list_i[start:actual_end] = ipart
                name = p_text[ipart].strip() if ipart < len(p_text) else f"part_{ipart}"
                tmp_list_n[start:actual_end] = name
                start = actual_end
        
        self.arrays[f"element_{elem_type}_part_indexes"] = tmp_list_i
        self.arrays[f"element_{elem_type}_part_ids"] = tmp_list_n
    
    def get_time(self) -> float:
        """Get the simulation time for this animation state."""
        return float(self.raw_header.get("time", 0.0))
    
    def get_node_coordinates(self) -> Optional[np.ndarray]:
        """Get node coordinates array (n_nodes × 3)."""
        return self.arrays.get("node_coordinates")
    
    def get_shell_connectivity(self) -> Optional[np.ndarray]:
        """Get shell element connectivity (n_shells × 4), 0-based."""
        return self.arrays.get("element_shell_node_indexes")
    
    def get_solid_connectivity(self) -> Optional[np.ndarray]:
        """Get solid element connectivity (n_solids × 8), 0-based."""
        return self.arrays.get("element_solid_node_indexes")
    
    def get_beam_connectivity(self) -> Optional[np.ndarray]:
        """Get beam element connectivity (n_beams × 2), 0-based."""
        return self.arrays.get("element_beam_node_indexes")
    
    def get_sph_connectivity(self) -> Optional[np.ndarray]:
        """Get SPH element connectivity (n_sph × 1), 0-based."""
        return self.arrays.get("sph_node_indexes")
    
    def get_nodal_scalars(self) -> Dict[str, np.ndarray]:
        """Get all nodal scalar arrays."""
        return {k: v for k, v in self.arrays.items() if k.startswith("node_") and v.ndim == 1}
    
    def get_nodal_vectors(self) -> Dict[str, np.ndarray]:
        """Get all nodal vector arrays."""
        return {k: v for k, v in self.arrays.items() if k.startswith("node_") and v.ndim == 2 and v.shape[1] == 3}
    
    def get_element_scalars(self) -> Dict[str, np.ndarray]:
        """Get all element scalar arrays."""
        return {k: v for k, v in self.arrays.items() if k.startswith("element_") and v.ndim == 1}
    
    def to_pyvista(self) -> pv.MultiBlock:
        """
        Convert the animation data to PyVista datasets.
        
        Returns
        -------
        pv.MultiBlock
            MultiBlock dataset containing:
            - 'shells': PolyData with quad shell elements
            - 'solids': UnstructuredGrid with hexahedral elements
            - 'beams': PolyData with line elements
            - 'sph': PolyData with vertex elements
        """
        blocks = pv.MultiBlock()
        coords = self.get_node_coordinates()
        
        if coords is None or len(coords) == 0:
            warnings.warn("No node coordinates found in animation file")
            return blocks
        
        # ========== SHELL ELEMENTS (2D) ==========
        shell_conn = self.get_shell_connectivity()
        if shell_conn is not None and len(shell_conn) > 0:
            # Convert quads to triangles for better VTK compatibility
            # or keep as quads using PolyData
            n_shells = len(shell_conn)
            
            # Create faces array for PolyData: [n_pts, i0, i1, i2, i3, ...]
            faces = np.column_stack([
                np.full(n_shells, 4, dtype=np.int64),
                shell_conn.astype(np.int64)
            ]).flatten()
            
            shells = pv.PolyData(coords, faces)
            shells.field_data["time"] = [self.get_time()]
            shells.field_data["element_type"] = ["shell"]
            
            # Add node IDs
            if "node_ids" in self.arrays:
                shells.point_data["node_ids"] = self.arrays["node_ids"]
            
            # Add nodal scalars
            for name, arr in self.get_nodal_scalars().items():
                if len(arr) == shells.n_points:
                    shells.point_data[name.replace("node_", "")] = arr
            
            # Add nodal vectors
            for name, arr in self.get_nodal_vectors().items():
                if len(arr) == shells.n_points:
                    shells.point_data[name.replace("node_", "")] = arr
            
            # Add element scalars
            for name, arr in self.get_element_scalars().items():
                if "shell" in name and len(arr) == shells.n_cells:
                    clean_name = name.replace("element_shell_", "")
                    shells.cell_data[clean_name] = arr
            
            # Add part IDs
            if "element_shell_part_indexes" in self.arrays:
                shells.cell_data["part_index"] = self.arrays["element_shell_part_indexes"]
            if "element_shell_part_ids" in self.arrays:
                # Convert to categorical/coded array
                part_ids = self.arrays["element_shell_part_ids"]
                unique_parts = np.unique(part_ids)
                part_map = {p: i for i, p in enumerate(unique_parts)}
                shells.cell_data["part_id"] = [part_map.get(p, -1) for p in part_ids]
                shells.field_data["part_names"] = list(unique_parts)
            
            # Add element IDs
            if "element_shell_ids" in self.arrays:
                shells.cell_data["element_id"] = self.arrays["element_shell_ids"]
            
            # Add deletion flag
            if "delEltA" in self.raw_arrays and len(self.raw_arrays["delEltA"]) == shells.n_cells:
                shells.cell_data["deleted"] = self.raw_arrays["delEltA"].astype(int)
            
            blocks["shells"] = shells
        
        # ========== SOLID ELEMENTS (3D) ==========
        solid_conn = self.get_solid_connectivity()
        if solid_conn is not None and len(solid_conn) > 0:
            n_solids = len(solid_conn)
            
            # VTK_HEXAHEDRON = 12
            cell_types = np.full(n_solids, 12, dtype=np.uint8)
            cells = np.column_stack([
                np.full(n_solids, 8, dtype=np.int64),
                solid_conn.astype(np.int64)
            ]).flatten()
            
            solids = pv.UnstructuredGrid(cells, cell_types, coords)
            solids.field_data["time"] = [self.get_time()]
            solids.field_data["element_type"] = ["solid"]
            
            # Add node IDs
            if "node_ids" in self.arrays:
                solids.point_data["node_ids"] = self.arrays["node_ids"]
            
            # Add nodal scalars
            for name, arr in self.get_nodal_scalars().items():
                if len(arr) == solids.n_points:
                    solids.point_data[name.replace("node_", "")] = arr
            
            # Add nodal vectors
            for name, arr in self.get_nodal_vectors().items():
                if len(arr) == solids.n_points:
                    solids.point_data[name.replace("node_", "")] = arr
            
            # Add element scalars (3D)
            if "eFunc3DA" in self.raw_arrays and self.raw_header["nbEFunc3D"] > 0:
                for ifun in range(self.raw_header["nbEFunc3D"]):
                    name = self.raw_arrays["fText3DA"][ifun].strip()
                    start = ifun * self.raw_header["nbElts3D"]
                    end = (ifun + 1) * self.raw_header["nbElts3D"]
                    solids.cell_data[name] = self.raw_arrays["eFunc3DA"][start:end]
            
            # Add element tensors (3D)
            if "tensVal3DA" in self.raw_arrays and self.raw_header["nbTens3D"] > 0:
                for itens in range(self.raw_header["nbTens3D"]):
                    name = self.raw_arrays["tText3DA"][itens].strip()
                    start = itens * self.raw_header["nbElts3D"]
                    end = (itens + 1) * self.raw_header["nbElts3D"]
                    solids.cell_data[name] = self.raw_arrays["tensVal3DA"][start:end]
            
            # Add part IDs
            if "element_solid_part_indexes" in self.arrays:
                solids.cell_data["part_index"] = self.arrays["element_solid_part_indexes"]
            if "element_solid_part_ids" in self.arrays:
                part_ids = self.arrays["element_solid_part_ids"]
                unique_parts = np.unique(part_ids)
                part_map = {p: i for i, p in enumerate(unique_parts)}
                solids.cell_data["part_id"] = [part_map.get(p, -1) for p in part_ids]
                solids.field_data["part_names"] = list(unique_parts)
            
            # Add element IDs
            if "element_solid_ids" in self.arrays:
                solids.cell_data["element_id"] = self.arrays["element_solid_ids"]
            
            # Add deletion flag
            if "delElt3DA" in self.raw_arrays and len(self.raw_arrays["delElt3DA"]) == solids.n_cells:
                solids.cell_data["deleted"] = self.raw_arrays["delElt3DA"].astype(int)
            
            blocks["solids"] = solids
        
        # ========== BEAM ELEMENTS (1D) ==========
        beam_conn = self.get_beam_connectivity()
        if beam_conn is not None and len(beam_conn) > 0:
            n_beams = len(beam_conn)
            
            # Create lines
            lines = np.column_stack([
                np.full(n_beams, 2, dtype=np.int64),
                beam_conn.astype(np.int64)
            ]).flatten()
            
            beams = pv.PolyData(coords, lines=lines)
            beams.field_data["time"] = [self.get_time()]
            beams.field_data["element_type"] = ["beam"]
            
            # Add node IDs
            if "node_ids" in self.arrays:
                beams.point_data["node_ids"] = self.arrays["node_ids"]
            
            # Add nodal scalars
            for name, arr in self.get_nodal_scalars().items():
                if len(arr) == beams.n_points:
                    beams.point_data[name.replace("node_", "")] = arr
            
            # Add nodal vectors
            for name, arr in self.get_nodal_vectors().items():
                if len(arr) == beams.n_points:
                    beams.point_data[name.replace("node_", "")] = arr
            
            # Add beam element data
            if "eFunc1DA" in self.raw_arrays and self.raw_header["nbEFunc1D"] > 0:
                for ifun in range(self.raw_header["nbEFunc1D"]):
                    name = self.raw_arrays["fText1DA"][ifun].strip()
                    start = ifun * self.raw_header["nbElts1D"]
                    end = (ifun + 1) * self.raw_header["nbElts1D"]
                    beams.cell_data[name] = self.raw_arrays["eFunc1DA"][start:end]
            
            # Add torsor data (forces/moments)
            if "torsVal1DA" in self.raw_arrays and self.raw_header["nbTors1D"] > 0:
                for itors in range(self.raw_header["nbTors1D"]):
                    name = self.raw_arrays["tText1DA"][itors].strip()
                    start = itors * self.raw_header["nbElts1D"]
                    end = (itors + 1) * self.raw_header["nbElts1D"]
                    # 9 components: [Fx, Fy, Fz, Mx, My, Mz, ...]
                    tors_data = self.raw_arrays["torsVal1DA"][start:end]
                    for icomp, comp_name in enumerate(["Fx", "Fy", "Fz", "Mx", "My", "Mz", "Sx", "Sy", "Sz"]):
                        beams.cell_data[f"{name}_{comp_name}"] = tors_data[:, icomp]
            
            # Add part IDs
            if "element_beam_part_indexes" in self.arrays:
                beams.cell_data["part_index"] = self.arrays["element_beam_part_indexes"]
            if "element_beam_part_ids" in self.arrays:
                part_ids = self.arrays["element_beam_part_ids"]
                unique_parts = np.unique(part_ids)
                part_map = {p: i for i, p in enumerate(unique_parts)}
                beams.cell_data["part_id"] = [part_map.get(p, -1) for p in part_ids]
                beams.field_data["part_names"] = list(unique_parts)
            
            # Add element IDs
            if "element_beam_ids" in self.arrays:
                beams.cell_data["element_id"] = self.arrays["element_beam_ids"]
            
            # Add deletion flag
            if "delElt1DA" in self.raw_arrays and len(self.raw_arrays["delElt1DA"]) == beams.n_cells:
                beams.cell_data["deleted"] = self.raw_arrays["delElt1DA"].astype(int)
            
            blocks["beams"] = beams
        
        # ========== SPH ELEMENTS ==========
        sph_conn = self.get_sph_connectivity()
        if sph_conn is not None and len(sph_conn) > 0:
            n_sph = len(sph_conn)
            
            # SPH particles as vertices
            verts = np.column_stack([
                np.full(n_sph, 1, dtype=np.int64),
                sph_conn.astype(np.int64)
            ]).flatten()
            
            sph = pv.PolyData(coords, verts)
            sph.field_data["time"] = [self.get_time()]
            sph.field_data["element_type"] = ["sph"]
            
            # Add node IDs
            if "node_ids" in self.arrays:
                sph.point_data["node_ids"] = self.arrays["node_ids"]
            
            # Add SPH-specific data
            if "eFuncSPH" in self.raw_arrays and self.raw_header["nbEFuncSPH"] > 0:
                for ifun in range(self.raw_header["nbEFuncSPH"]):
                    name = self.raw_arrays["scalTextSPH"][ifun].strip()
                    start = ifun * self.raw_header["nbEltsSPH"]
                    end = (ifun + 1) * self.raw_header["nbEltsSPH"]
                    sph.cell_data[name] = self.raw_arrays["eFuncSPH"][start:end]
            
            # Add SPH tensors
            if "tensValSPH" in self.raw_arrays and self.raw_header["nbTensSPH"] > 0:
                for itens in range(self.raw_header["nbTensSPH"]):
                    name = self.raw_arrays["tensTextSPH"][itens].strip()
                    start = itens * self.raw_header["nbEltsSPH"]
                    end = (itens + 1) * self.raw_header["nbEltsSPH"]
                    sph.cell_data[name] = self.raw_arrays["tensValSPH"][start:end]
            
            # Add element IDs
            if "element_sph_ids" in self.arrays:
                sph.cell_data["element_id"] = self.arrays["element_sph_ids"]
            
            # Add deletion flag
            if "delEltSPH" in self.raw_arrays and len(self.raw_arrays["delEltSPH"]) == sph.n_cells:
                sph.cell_data["deleted"] = self.raw_arrays["delEltSPH"].astype(int)
            
            blocks["sph"] = sph
        
        return blocks


def read_radioss_anim(filepath: Union[str, Path]) -> pv.MultiBlock:
    """
    Read a Radioss animation file and return a PyVista MultiBlock dataset.
    
    Parameters
    ----------
    filepath : str or Path
        Path to the .anim file
        
    Returns
    -------
    pv.MultiBlock
        MultiBlock containing shells, solids, beams, and SPH as separate blocks
        
    Examples
    --------
    >>> import pyvista as pv
    >>> from radioss_reader import read_radioss_anim
    >>> mesh = read_radioss_anim("animation.anim")
    >>> mesh["shells"].plot()
    """
    reader = RadiossAnimReader(filepath)
    return reader.to_pyvista()


def register_pyvista_reader() -> None:
    """
    Register the Radioss animation reader with PyVista's reader registry.
    
    After calling this function, you can use:
        pv.read("file.anim")
    """
    try:
        pv.register_reader(".anim", read_radioss_anim)
        print("Registered Radioss .anim reader with PyVista")
    except Exception as e:
        warnings.warn(f"Could not register reader: {e}")


# Convenience function for entry points
def _pyvista_reader_entry(path: str, **kwargs) -> pv.MultiBlock:
    """Entry point for PyVista reader registration."""
    return read_radioss_anim(path)


if __name__ == "__main__":
    # Example usage
    import sys
    
    if len(sys.argv) > 1:
        filepath = sys.argv[1]
        print(f"Reading: {filepath}")
        
        reader = RadiossAnimReader(filepath)
        print(f"Time: {reader.get_time()}")
        print(f"Nodes: {reader.raw_header.get('nbNodes', 0)}")
        print(f"Shells: {reader.raw_header.get('nbFacets', 0)}")
        print(f"Solids: {reader.raw_header.get('nbElts3D', 0)}")
        print(f"Beams: {reader.raw_header.get('nbElts1D', 0)}")
        print(f"SPH: {reader.raw_header.get('nbEltsSPH', 0)}")
        
        mesh = reader.to_pyvista()
        print(f"\nBlocks: {list(mesh.keys())}")
        
        for name, block in mesh.items():
            print(f"\n{name}:")
            print(f"  Points: {block.n_points}")
            print(f"  Cells: {block.n_cells}")
            print(f"  Point data: {list(block.point_data.keys())}")
            print(f"  Cell data: {list(block.cell_data.keys())}")
    else:
        print("Usage: python radioss_reader.py <file.anim>")

Setup / Installation

Create a pyproject.toml for your package:

[project]
name = "pyvista-radioss"
version = "0.1.0"
description = "Radioss animation file reader for PyVista"
dependencies = [
    "numpy",
    "pyvista",
]

[project.entry-points."pyvista.readers"]
".anim" = "radioss_reader:_pyvista_reader_entry"

Usage Examples

import pyvista as pv

# Method 1: Direct import
from radioss_reader import read_radioss_anim
mesh = read_radioss_anim("model.anim")

# Method 2: Register with PyVista
from radioss_reader import register_pyvista_reader
register_pyvista_reader()
mesh = pv.read("model.anim")

# Method 3: After pip install (entry points auto-register)
mesh = pv.read("model.anim")

# Access individual element types
shells = mesh["shells"]
solids = mesh["solids"]
beams = mesh["beams"]

# Plot with results
shells.plot(scalars="von_mises", cmap="jet")
solids.plot(scalars="plastic_strain", cmap="viridis")

# Animate multiple time states
import glob
files = sorted(glob.glob("animation*.anim"))
plotter = pv.Plotter()
for f in files:
    state = read_radioss_anim(f)["shells"]
    plotter.add_mesh(state, scalars="stress", show_edges=True)
    plotter.write_frame()

This implementation provides a complete, standalone Radioss animation file reader that integrates seamlessly with PyVista's ecosystem, supporting all major element types (shells, solids, beams, SPH) with their associated scalar, vector, and tensor data.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment