API Reference

Mechanics

SymbolicHandler

A class that handles symbolic computations for Continuum Mechanics Hyperelastic Frameworks using SymPy.

Attributes:

Name	Type	Description
c_tensor	Matrix	A 3x3 matrix of symbols.

Source code in hyper_surrogate/mechanics/symbolic.py

class SymbolicHandler:
    """
    A class that handles symbolic computations for Continuum Mechanics Hyperelastic Frameworks using SymPy.

    Attributes:
        c_tensor (Matrix): A 3x3 matrix of symbols.
    """

    def __init__(self) -> None:
        self.c_tensor = self._c_tensor()
        self.f_tensor = self._f_tensor()

    def f_symbols(self) -> list[Symbol]:
        """
        Return the f_tensor flattened symbols.

        Returns:
            list: A list of f_tensor symbols.
        """
        return [self.f_tensor[i, j] for i in range(3) for j in range(3)]

    def _f_tensor(self) -> ImmutableMatrix:
        """
        Create a 3x3 matrix of symbols.

        Returns:
            Matrix: A 3x3 matrix of symbols.
        """
        return ImmutableMatrix(3, 3, lambda i, j: Symbol(f"F_{i + 1}{j + 1}"))

    def c_symbols(self) -> list[Symbol]:
        """
        Return the c_tensor flattened symbols.

        Returns:
            list: A list of c_tensor symbols.
        """
        return [self.c_tensor[i, j] for i in range(3) for j in range(3)]

    def _c_tensor(self) -> ImmutableMatrix:
        """
        Create a 3x3 matrix of symbols.

        Returns:
            Matrix: A 3x3 matrix of symbols.
        """
        return ImmutableMatrix(3, 3, lambda i, j: Symbol(f"C_{i + 1}{j + 1}"))

    # multiply c_tensor by itself
    def _c_tensor_squared(self) -> ImmutableMatrix:
        """
        Compute the square of the c_tensor.

        Returns:
            Matrix: The square of the c_tensor.
        """
        # matrix product
        return self.c_tensor**2

    @property
    def isochoric_invariant1(self) -> Expr:
        """
        Compute the first isochoric invariant of the c_tensor.

        Returns:
            Expr: The first isochoric invariant of the c_tensor.
        """
        return self.c_tensor.trace() * (self.invariant3 ** (-Rational(1, 3)))

    @property
    def isochoric_invariant2(self) -> Expr:
        """
        Compute the second isochoric invariant of the c_tensor.
        I2_bar = 1/2 * (I1^2 - trace(C^2)) * J^{-4/3}

        Returns:
            Expr: The second isochoric invariant of the c_tensor.
        """
        i1 = self.c_tensor.trace()
        return Rational(1, 2) * (i1**2 - self._c_tensor_squared().trace()) * (self.invariant3 ** (-Rational(2, 3)))

    @property
    def invariant3(self) -> Expr:
        """
        Compute the third invariant of the c_tensor.

        Returns:
            Expr: The third invariant of the c_tensor.
        """
        return self.c_tensor.det()

    def pk2(self, sef: Expr) -> Matrix:
        """
        Compute the pk2 tensor.

        Args:
            sef (Expr): The strain energy function.

        Returns:
            Matrix: The pk2 tensor.
        """
        pk2 = Matrix([[diff(sef, self.c_tensor[i, j]) for j in range(3)] for i in range(3)])
        # force symmetry
        return pk2 + pk2.T

    def cmat(self, pk2: Matrix) -> ImmutableDenseNDimArray:
        """
        Compute the cmat tensor.

        Args:
            pk2 (Matrix): The pk2 tensor.

        Returns:
            ImmutableDenseNDimArray: The stiffness tensor (3x3x3x3) with minor symmetry.
        """
        return ImmutableDenseNDimArray([
            [
                [
                    [(diff(pk2[i, j], self.c_tensor[k, ll]) + diff(pk2[i, j], self.c_tensor[ll, k])) for ll in range(3)]
                    for k in range(3)
                ]
                for j in range(3)
            ]
            for i in range(3)
        ])

    def cauchy(self, sef: Expr, f: Matrix) -> Matrix:
        """
        Compute the Cauchy stress tensor.

        Args:
            sef (Expr): The strain energy function.
            f (Matrix): The deformation gradient tensor.

        Returns:
            Matrix: The Cauchy stress tensor.
        """
        return self.pushforward_2nd_order(self.pk2(sef), f)

    def spatial_tangent(self, pk2: Matrix, f: Matrix) -> ImmutableDenseNDimArray:
        """
        Compute the spatial tangent stiffness tensor.

        Args:
            pk2 (Matrix): The pk2 tensor.
            f (Matrix): The deformation gradient tensor.

        Returns:
            ImmutableDenseNDimArray: The spatial tangent stiffness tensor.
        """
        return self.pushforward_4th_order(self.cmat(pk2), f)

    def substitute(
        self,
        symbolic_tensor: MutableDenseMatrix,
        numerical_c_tensor: np.ndarray,
        *args: dict,
    ) -> Matrix:
        """
        Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

        Args:
            symbolic_tensor (Matrix): A symbolic tensor to substitute numerical values into.
            numerical_c_tensor (np.ndarray): A 3x3 numerical matrix to substitute into c_tensor.
            args (dict): Additional substitution dictionaries.

        Returns:
            Matrix: The symbolic_tensor with numerical values substituted.

        Raises:
            ValueError: If numerical_tensor is not a 3x3 matrix.
        """
        if not isinstance(numerical_c_tensor, np.ndarray) or numerical_c_tensor.shape != (3, 3):
            raise ValueError("c_tensor.shape")

        # Start with substitutions for c_tensor elements
        substitutions = {self.c_tensor[i, j]: numerical_c_tensor[i, j] for i in range(3) for j in range(3)}
        # Merge additional substitution dictionaries from *args
        substitutions.update(*args)
        return symbolic_tensor.subs(substitutions)

    def substitute_iterator(
        self,
        symbolic_tensor: MutableDenseMatrix,
        numerical_c_tensors: np.ndarray,
        *args: dict,
    ) -> Generator[Matrix, None, None]:
        """
        Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

        Args:
            symbolic_tensor (Matrix): A symbolic tensor to substitute numerical values into.
            numerical_c_tensors (np.ndarray): N 3x3 numerical matrices to substitute into c_tensor.
            args (dict): Additional substitution dictionaries.

        Returns:
            Generator[Matrix, None, None]: The symbolic_tensor with numerical values substituted.

        Raises:
            ValueError: If numerical_tensor is not a 3x3 matrix.
        """
        for numerical_c_tensor in numerical_c_tensors:
            yield self.substitute(symbolic_tensor, numerical_c_tensor, *args)

    def lambdify(self, symbolic_tensor: Matrix, *args: Iterable[Any]) -> Any:
        """
        Create a lambdified function from a symbolic tensor that can be used for numerical evaluation.

        Args:
            symbolic_tensor (Expr or Matrix): The symbolic tensor to be lambdified.
            args (dict): Additional substitution lists of symbols.
        Returns:
            function: A function that can be used to numerically evaluate the tensor with specific values.
        """

        return lambdify((self.c_symbols(), *args), symbolic_tensor.tolist(), modules="numpy")

    def _evaluate(self, lambdified_tensor: Any, *args: Any, **kwargs: Any) -> Any:
        """
        Evaluate a lambdified tensor with specific values.

        Args:
            lambdified_tensor (function): A lambdified tensor function.
            args (dict): Additional substitution lists of symbols.
            kwargs (dict): Additional keyword arguments.

        Returns:
            Generator[Any, None, None]: The evaluated tensor.
        """
        return lambdified_tensor(*args, **kwargs)

    def evaluate_iterator(
        self, lambdified_tensor: Any, numerical_c_tensors: np.ndarray, *args: Any, **kwargs: Any
    ) -> Generator[Any, None, None]:
        """
        Evaluate a lambdified tensor with specific values.

        Args:
            lambdified_tensor (function): A lambdified tensor function.
            args (dict): Additional substitution lists of symbols.
            kwargs (dict): Additional keyword arguments.

        Returns:
            Generator[Any, None, None]: The evaluated tensor.
        """
        for numerical_c_tensor in numerical_c_tensors:
            yield self._evaluate(lambdified_tensor, numerical_c_tensor.flatten(), *args, **kwargs)

    @staticmethod
    def to_voigt_2(tensor: Matrix) -> Matrix:
        """
        Convert a 3x3 matrix to 6x1 matrix using Voigt notation

        Args:
            tensor (sp.Matrix): A 3x3 symmetric matrix.

        Returns:
            sp.Matrix: A 6x1 matrix.
        """
        # Validate the input tensor dimensions
        if tensor.shape != (3, 3):
            raise ValueError("Wrong.shape.")
        # Voigt notation conversion: xx, yy, zz, xy, xz, yz
        return Matrix([
            tensor[0, 0],  # xx
            tensor[1, 1],  # yy
            tensor[2, 2],  # zz
            tensor[0, 1],  # xy
            tensor[0, 2],  # xz
            tensor[1, 2],  # yz
        ])

    @staticmethod
    def to_voigt_4(tensor: ImmutableDenseNDimArray) -> Matrix:
        """
        Convert a 3x3x3x3 matrix to 6x6 matrix using Voigt notation

        Args:
            tensor (ImmutableDenseNDimArray): A 3x3x3x3 matrix.

        Returns:
            Matrix: A 6x6 matrix.
        """
        # Validate the input tensor dimensions
        if tensor.shape != (3, 3, 3, 3):
            raise ValueError("Wrong.shape.")

        # Voigt notation index mapping
        ii1 = [0, 1, 2, 0, 0, 1]
        ii2 = [0, 1, 2, 1, 2, 2]

        # Initialize a 6x6 matrix
        voigt_matrix = Matrix.zeros(6, 6)

        # Fill the Voigt matrix by averaging symmetric components of the 4th-order tensor
        for i in range(6):
            for j in range(6):
                # Access the relevant tensor components using ii1 and ii2 mappings
                pp1 = tensor[ii1[i], ii2[i], ii1[j], ii2[j]]
                pp2 = tensor[ii1[i], ii2[i], ii2[j], ii1[j]]
                # Average the two permutations to ensure symmetry
                voigt_matrix[i, j] = 0.5 * (pp1 + pp2)
        return voigt_matrix

    @staticmethod
    def pushforward_2nd_order(tensor2: Matrix, f: Matrix) -> Matrix:
        """
        Push forward a 2nd order tensor in material configuration.

        args:
        tensor2: Any - The 2nd order tensor
        f: Any - The deformation gradient tensor

        returns:
        Any - The pushforwarded 2nd order tensor
        """
        return simplify(f * tensor2 * f.T)

    @staticmethod
    def pushforward_4th_order(tensor4: ImmutableDenseNDimArray, f: Matrix) -> ImmutableDenseNDimArray:
        """
        Push forward a 4th order tensor in material configuration.

        Args:
            tensor4 (MutableDenseNDimArray): The 4th order tensor.
            f (Matrix): The deformation gradient tensor (2nd order tensor).

        Returns:
            ImmutableDenseNDimArray: The pushforwarded 4th order tensor.
        """
        # Calculate the pushforwarded tensor using comprehensions and broadcasting
        return ImmutableDenseNDimArray([
            [
                [
                    [
                        sum(
                            f[i, ii] * f[j, jj] * f[k, kk] * f[l0, ll] * tensor4[ii, jj, kk, ll]
                            for ii in range(3)
                            for jj in range(3)
                            for kk in range(3)
                            for ll in range(3)
                        )
                        for l0 in range(3)
                    ]
                    for k in range(3)
                ]
                for j in range(3)
            ]
            for i in range(3)
        ])

    @staticmethod
    def jr(sigma: Matrix) -> ImmutableDenseNDimArray:
        """
        Compute the Jaumann rate contribution for the spatial elasticity tensor.

        Args:
            sigma (Matrix): The Cauchy stress tensor (2nd order tensor).

        Returns:
            ImmutableDenseNDimArray: The Jaumann rate contribution (4th order tensor).
        """
        # Ensure sigma is a 3x3 matrix
        if sigma.shape != (3, 3):
            raise ValueError("Wrongshape")

        # Use list comprehension to build the 4th-order tensor directly
        return ImmutableDenseNDimArray([
            [
                [
                    [
                        (1 / 2)
                        * (
                            int(i == k) * sigma[j, ll]
                            + sigma[i, k] * int(j == ll)
                            + int(i == ll) * sigma[j, k]
                            + sigma[i, ll] * int(j == k)
                        )
                        for ll in range(3)
                    ]
                    for k in range(3)
                ]
                for j in range(3)
            ]
            for i in range(3)
        ])

    # Alias for jr
    jaumann_correction = jr

invariant3 property

Compute the third invariant of the c_tensor.

Returns:

Name	Type	Description
Expr	Expr	The third invariant of the c_tensor.

isochoric_invariant1 property

Compute the first isochoric invariant of the c_tensor.

Returns:

Name	Type	Description
Expr	Expr	The first isochoric invariant of the c_tensor.

isochoric_invariant2 property

Compute the second isochoric invariant of the c_tensor. I2_bar = ½ * (I1^2 - trace(C^2)) * J^{-4/3}

Returns:

Name	Type	Description
Expr	Expr	The second isochoric invariant of the c_tensor.

c_symbols()

Return the c_tensor flattened symbols.

Returns:

Name	Type	Description
list	list[Symbol]	A list of c_tensor symbols.

Source code in hyper_surrogate/mechanics/symbolic.py

def c_symbols(self) -> list[Symbol]:
    """
    Return the c_tensor flattened symbols.

    Returns:
        list: A list of c_tensor symbols.
    """
    return [self.c_tensor[i, j] for i in range(3) for j in range(3)]

cauchy(sef, f)

Compute the Cauchy stress tensor.

Parameters:

Name	Type	Description	Default
sef	Expr	The strain energy function.	required
f	Matrix	The deformation gradient tensor.	required

Returns:

Name	Type	Description
Matrix	Matrix	The Cauchy stress tensor.

Source code in hyper_surrogate/mechanics/symbolic.py

def cauchy(self, sef: Expr, f: Matrix) -> Matrix:
    """
    Compute the Cauchy stress tensor.

    Args:
        sef (Expr): The strain energy function.
        f (Matrix): The deformation gradient tensor.

    Returns:
        Matrix: The Cauchy stress tensor.
    """
    return self.pushforward_2nd_order(self.pk2(sef), f)

cmat(pk2)

Compute the cmat tensor.

Parameters:

Name	Type	Description	Default
pk2	Matrix	The pk2 tensor.	required

Returns:

Name	Type	Description
ImmutableDenseNDimArray	ImmutableDenseNDimArray	The stiffness tensor (3x3x3x3) with minor symmetry.

Source code in hyper_surrogate/mechanics/symbolic.py

def cmat(self, pk2: Matrix) -> ImmutableDenseNDimArray:
    """
    Compute the cmat tensor.

    Args:
        pk2 (Matrix): The pk2 tensor.

    Returns:
        ImmutableDenseNDimArray: The stiffness tensor (3x3x3x3) with minor symmetry.
    """
    return ImmutableDenseNDimArray([
        [
            [
                [(diff(pk2[i, j], self.c_tensor[k, ll]) + diff(pk2[i, j], self.c_tensor[ll, k])) for ll in range(3)]
                for k in range(3)
            ]
            for j in range(3)
        ]
        for i in range(3)
    ])

evaluate_iterator(lambdified_tensor, numerical_c_tensors, *args, **kwargs)

Evaluate a lambdified tensor with specific values.

Parameters:

Name	Type	Description	Default
lambdified_tensor	function	A lambdified tensor function.	required
args	dict	Additional substitution lists of symbols.	()
kwargs	dict	Additional keyword arguments.	{}

Returns:

Type	Description
None	Generator[Any, None, None]: The evaluated tensor.

Source code in hyper_surrogate/mechanics/symbolic.py

def evaluate_iterator(
    self, lambdified_tensor: Any, numerical_c_tensors: np.ndarray, *args: Any, **kwargs: Any
) -> Generator[Any, None, None]:
    """
    Evaluate a lambdified tensor with specific values.

    Args:
        lambdified_tensor (function): A lambdified tensor function.
        args (dict): Additional substitution lists of symbols.
        kwargs (dict): Additional keyword arguments.

    Returns:
        Generator[Any, None, None]: The evaluated tensor.
    """
    for numerical_c_tensor in numerical_c_tensors:
        yield self._evaluate(lambdified_tensor, numerical_c_tensor.flatten(), *args, **kwargs)

f_symbols()

Return the f_tensor flattened symbols.

Returns:

Name	Type	Description
list	list[Symbol]	A list of f_tensor symbols.

Source code in hyper_surrogate/mechanics/symbolic.py

def f_symbols(self) -> list[Symbol]:
    """
    Return the f_tensor flattened symbols.

    Returns:
        list: A list of f_tensor symbols.
    """
    return [self.f_tensor[i, j] for i in range(3) for j in range(3)]

jr(sigma) staticmethod

Compute the Jaumann rate contribution for the spatial elasticity tensor.

Parameters:

Name	Type	Description	Default
sigma	Matrix	The Cauchy stress tensor (2nd order tensor).	required

Returns:

Name	Type	Description
ImmutableDenseNDimArray	ImmutableDenseNDimArray	The Jaumann rate contribution (4th order tensor).

Source code in hyper_surrogate/mechanics/symbolic.py

@staticmethod
def jr(sigma: Matrix) -> ImmutableDenseNDimArray:
    """
    Compute the Jaumann rate contribution for the spatial elasticity tensor.

    Args:
        sigma (Matrix): The Cauchy stress tensor (2nd order tensor).

    Returns:
        ImmutableDenseNDimArray: The Jaumann rate contribution (4th order tensor).
    """
    # Ensure sigma is a 3x3 matrix
    if sigma.shape != (3, 3):
        raise ValueError("Wrongshape")

    # Use list comprehension to build the 4th-order tensor directly
    return ImmutableDenseNDimArray([
        [
            [
                [
                    (1 / 2)
                    * (
                        int(i == k) * sigma[j, ll]
                        + sigma[i, k] * int(j == ll)
                        + int(i == ll) * sigma[j, k]
                        + sigma[i, ll] * int(j == k)
                    )
                    for ll in range(3)
                ]
                for k in range(3)
            ]
            for j in range(3)
        ]
        for i in range(3)
    ])

lambdify(symbolic_tensor, *args)

Create a lambdified function from a symbolic tensor that can be used for numerical evaluation.

Parameters:

Name	Type	Description	Default
symbolic_tensor	Expr or Matrix	The symbolic tensor to be lambdified.	required
args	dict	Additional substitution lists of symbols.	()

Returns: function: A function that can be used to numerically evaluate the tensor with specific values.

Source code in hyper_surrogate/mechanics/symbolic.py

def lambdify(self, symbolic_tensor: Matrix, *args: Iterable[Any]) -> Any:
    """
    Create a lambdified function from a symbolic tensor that can be used for numerical evaluation.

    Args:
        symbolic_tensor (Expr or Matrix): The symbolic tensor to be lambdified.
        args (dict): Additional substitution lists of symbols.
    Returns:
        function: A function that can be used to numerically evaluate the tensor with specific values.
    """

    return lambdify((self.c_symbols(), *args), symbolic_tensor.tolist(), modules="numpy")

pk2(sef)

Compute the pk2 tensor.

Parameters:

Name	Type	Description	Default
sef	Expr	The strain energy function.	required

Returns:

Name	Type	Description
Matrix	Matrix	The pk2 tensor.

Source code in hyper_surrogate/mechanics/symbolic.py

def pk2(self, sef: Expr) -> Matrix:
    """
    Compute the pk2 tensor.

    Args:
        sef (Expr): The strain energy function.

    Returns:
        Matrix: The pk2 tensor.
    """
    pk2 = Matrix([[diff(sef, self.c_tensor[i, j]) for j in range(3)] for i in range(3)])
    # force symmetry
    return pk2 + pk2.T

pushforward_2nd_order(tensor2, f) staticmethod

Push forward a 2^nd order tensor in material configuration.

args: tensor2: Any - The 2^nd order tensor f: Any - The deformation gradient tensor

returns: Any - The pushforwarded 2^nd order tensor

Source code in hyper_surrogate/mechanics/symbolic.py

@staticmethod
def pushforward_2nd_order(tensor2: Matrix, f: Matrix) -> Matrix:
    """
    Push forward a 2nd order tensor in material configuration.

    args:
    tensor2: Any - The 2nd order tensor
    f: Any - The deformation gradient tensor

    returns:
    Any - The pushforwarded 2nd order tensor
    """
    return simplify(f * tensor2 * f.T)

pushforward_4th_order(tensor4, f) staticmethod

Push forward a 4^th order tensor in material configuration.

Parameters:

Name	Type	Description	Default
tensor4	MutableDenseNDimArray	The 4th order tensor.	required
f	Matrix	The deformation gradient tensor (2nd order tensor).	required

Returns:

Name	Type	Description
ImmutableDenseNDimArray	ImmutableDenseNDimArray	The pushforwarded 4th order tensor.

Source code in hyper_surrogate/mechanics/symbolic.py

@staticmethod
def pushforward_4th_order(tensor4: ImmutableDenseNDimArray, f: Matrix) -> ImmutableDenseNDimArray:
    """
    Push forward a 4th order tensor in material configuration.

    Args:
        tensor4 (MutableDenseNDimArray): The 4th order tensor.
        f (Matrix): The deformation gradient tensor (2nd order tensor).

    Returns:
        ImmutableDenseNDimArray: The pushforwarded 4th order tensor.
    """
    # Calculate the pushforwarded tensor using comprehensions and broadcasting
    return ImmutableDenseNDimArray([
        [
            [
                [
                    sum(
                        f[i, ii] * f[j, jj] * f[k, kk] * f[l0, ll] * tensor4[ii, jj, kk, ll]
                        for ii in range(3)
                        for jj in range(3)
                        for kk in range(3)
                        for ll in range(3)
                    )
                    for l0 in range(3)
                ]
                for k in range(3)
            ]
            for j in range(3)
        ]
        for i in range(3)
    ])

spatial_tangent(pk2, f)

Compute the spatial tangent stiffness tensor.

Parameters:

Name	Type	Description	Default
pk2	Matrix	The pk2 tensor.	required
f	Matrix	The deformation gradient tensor.	required

Returns:

Name	Type	Description
ImmutableDenseNDimArray	ImmutableDenseNDimArray	The spatial tangent stiffness tensor.

Source code in hyper_surrogate/mechanics/symbolic.py

def spatial_tangent(self, pk2: Matrix, f: Matrix) -> ImmutableDenseNDimArray:
    """
    Compute the spatial tangent stiffness tensor.

    Args:
        pk2 (Matrix): The pk2 tensor.
        f (Matrix): The deformation gradient tensor.

    Returns:
        ImmutableDenseNDimArray: The spatial tangent stiffness tensor.
    """
    return self.pushforward_4th_order(self.cmat(pk2), f)

substitute(symbolic_tensor, numerical_c_tensor, *args)

Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

Parameters:

Name	Type	Description	Default
symbolic_tensor	Matrix	A symbolic tensor to substitute numerical values into.	required
numerical_c_tensor	ndarray	A 3x3 numerical matrix to substitute into c_tensor.	required
args	dict	Additional substitution dictionaries.	()

Returns:

Name	Type	Description
Matrix	Matrix	The symbolic_tensor with numerical values substituted.

Raises:

Type	Description
ValueError	If numerical_tensor is not a 3x3 matrix.

Source code in hyper_surrogate/mechanics/symbolic.py

def substitute(
    self,
    symbolic_tensor: MutableDenseMatrix,
    numerical_c_tensor: np.ndarray,
    *args: dict,
) -> Matrix:
    """
    Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

    Args:
        symbolic_tensor (Matrix): A symbolic tensor to substitute numerical values into.
        numerical_c_tensor (np.ndarray): A 3x3 numerical matrix to substitute into c_tensor.
        args (dict): Additional substitution dictionaries.

    Returns:
        Matrix: The symbolic_tensor with numerical values substituted.

    Raises:
        ValueError: If numerical_tensor is not a 3x3 matrix.
    """
    if not isinstance(numerical_c_tensor, np.ndarray) or numerical_c_tensor.shape != (3, 3):
        raise ValueError("c_tensor.shape")

    # Start with substitutions for c_tensor elements
    substitutions = {self.c_tensor[i, j]: numerical_c_tensor[i, j] for i in range(3) for j in range(3)}
    # Merge additional substitution dictionaries from *args
    substitutions.update(*args)
    return symbolic_tensor.subs(substitutions)

substitute_iterator(symbolic_tensor, numerical_c_tensors, *args)

Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

Parameters:

Name	Type	Description	Default
symbolic_tensor	Matrix	A symbolic tensor to substitute numerical values into.	required
numerical_c_tensors	ndarray	N 3x3 numerical matrices to substitute into c_tensor.	required
args	dict	Additional substitution dictionaries.	()

Returns:

Type	Description
None	Generator[Matrix, None, None]: The symbolic_tensor with numerical values substituted.

Raises:

Type	Description
ValueError	If numerical_tensor is not a 3x3 matrix.

Source code in hyper_surrogate/mechanics/symbolic.py

def substitute_iterator(
    self,
    symbolic_tensor: MutableDenseMatrix,
    numerical_c_tensors: np.ndarray,
    *args: dict,
) -> Generator[Matrix, None, None]:
    """
    Automatically substitute numerical values from a given 3x3 numerical matrix into c_tensor.

    Args:
        symbolic_tensor (Matrix): A symbolic tensor to substitute numerical values into.
        numerical_c_tensors (np.ndarray): N 3x3 numerical matrices to substitute into c_tensor.
        args (dict): Additional substitution dictionaries.

    Returns:
        Generator[Matrix, None, None]: The symbolic_tensor with numerical values substituted.

    Raises:
        ValueError: If numerical_tensor is not a 3x3 matrix.
    """
    for numerical_c_tensor in numerical_c_tensors:
        yield self.substitute(symbolic_tensor, numerical_c_tensor, *args)

to_voigt_2(tensor) staticmethod

Convert a 3x3 matrix to 6x1 matrix using Voigt notation

Parameters:

Name	Type	Description	Default
tensor	Matrix	A 3x3 symmetric matrix.	required

Returns:

Type	Description
Matrix	sp.Matrix: A 6x1 matrix.

Source code in hyper_surrogate/mechanics/symbolic.py

@staticmethod
def to_voigt_2(tensor: Matrix) -> Matrix:
    """
    Convert a 3x3 matrix to 6x1 matrix using Voigt notation

    Args:
        tensor (sp.Matrix): A 3x3 symmetric matrix.

    Returns:
        sp.Matrix: A 6x1 matrix.
    """
    # Validate the input tensor dimensions
    if tensor.shape != (3, 3):
        raise ValueError("Wrong.shape.")
    # Voigt notation conversion: xx, yy, zz, xy, xz, yz
    return Matrix([
        tensor[0, 0],  # xx
        tensor[1, 1],  # yy
        tensor[2, 2],  # zz
        tensor[0, 1],  # xy
        tensor[0, 2],  # xz
        tensor[1, 2],  # yz
    ])

to_voigt_4(tensor) staticmethod

Convert a 3x3x3x3 matrix to 6x6 matrix using Voigt notation

Parameters:

Name	Type	Description	Default
tensor	ImmutableDenseNDimArray	A 3x3x3x3 matrix.	required

Returns:

Name	Type	Description
Matrix	Matrix	A 6x6 matrix.

Source code in hyper_surrogate/mechanics/symbolic.py

@staticmethod
def to_voigt_4(tensor: ImmutableDenseNDimArray) -> Matrix:
    """
    Convert a 3x3x3x3 matrix to 6x6 matrix using Voigt notation

    Args:
        tensor (ImmutableDenseNDimArray): A 3x3x3x3 matrix.

    Returns:
        Matrix: A 6x6 matrix.
    """
    # Validate the input tensor dimensions
    if tensor.shape != (3, 3, 3, 3):
        raise ValueError("Wrong.shape.")

    # Voigt notation index mapping
    ii1 = [0, 1, 2, 0, 0, 1]
    ii2 = [0, 1, 2, 1, 2, 2]

    # Initialize a 6x6 matrix
    voigt_matrix = Matrix.zeros(6, 6)

    # Fill the Voigt matrix by averaging symmetric components of the 4th-order tensor
    for i in range(6):
        for j in range(6):
            # Access the relevant tensor components using ii1 and ii2 mappings
            pp1 = tensor[ii1[i], ii2[i], ii1[j], ii2[j]]
            pp2 = tensor[ii1[i], ii2[i], ii2[j], ii1[j]]
            # Average the two permutations to ensure symmetry
            voigt_matrix[i, j] = 0.5 * (pp1 + pp2)
    return voigt_matrix

Kinematics

A class that provides various kinematic methods.

Attributes:

Name	Type	Description
None		This class does not have any attributes.

Methods:

Name	Description
jacobian	Compute the Jacobian of the deformation gradient.
trace_invariant	Calculate the first invariant (trace) of each tensor in the batch.
quadratic_invariant	Calculate the second invariant of the tensor.
det_invariant	Calculate the third invariant (determinant) of the tensor.
isochoric_invariant1	Calculate the isochoric first invariant.
isochoric_invariant2	Calculate the isochoric second invariant.
right_cauchy_green	Compute the right Cauchy-Green deformation tensor for a batch of deformation gradients.
left_cauchy_green	Compute the left Cauchy-Green deformation tensor for a batch of deformation gradients.
rotation_tensor	Compute the rotation tensors.
pushforward	Forward tensor configuration.
principal_stretches	Compute the principal stretches.
principal_directions	Compute the principal directions.

Source code in hyper_surrogate/mechanics/kinematics.py

class Kinematics:
    """
    A class that provides various kinematic methods.

    Attributes:
        None: This class does not have any attributes.

    Methods:
        jacobian: Compute the Jacobian of the deformation gradient.
        trace_invariant: Calculate the first invariant (trace) of each tensor in the batch.
        quadratic_invariant: Calculate the second invariant of the tensor.
        det_invariant: Calculate the third invariant (determinant) of the tensor.
        isochoric_invariant1: Calculate the isochoric first invariant.
        isochoric_invariant2: Calculate the isochoric second invariant.
        right_cauchy_green: Compute the right Cauchy-Green deformation tensor for a batch of deformation gradients.
        left_cauchy_green: Compute the left Cauchy-Green deformation tensor for a batch of deformation gradients.
        rotation_tensor: Compute the rotation tensors.
        pushforward: Forward tensor configuration.
        principal_stretches: Compute the principal stretches.
        principal_directions: Compute the principal directions.

    """

    @staticmethod
    def jacobian(f: np.ndarray) -> Any:
        """
        Compute the Jacobian of the deformation gradient.

        Args:
            f: 4D tensor of shape (N, 3, 3).

        Returns:
            np.ndarray: The Jacobian of the deformation gradient.
        """
        return np.linalg.det(f)

    @staticmethod
    def trace_invariant(c: np.ndarray) -> Any:
        """
        Calculate the first invariant (trace) of each tensor in the batch.

        Args:
            c: Tensor of shape (N, 3, 3).

        Returns:
            The first invariant of each tensor in the batch.
        """
        return np.einsum("nii->n", c)

    @staticmethod
    def quadratic_invariant(c: np.ndarray) -> Any:
        """
        Calculate the second invariant of the tensor.

        Args:
            c: Tensor of shape (N, 3, 3).

        Returns:
            The second invariant.
        """
        return 0.5 * (np.einsum("nii->n", c) ** 2 - np.einsum("nij,nji->n", c, c))

    @staticmethod
    def det_invariant(c: np.ndarray) -> Any:
        """
        Calculate the third invariant (determinant) of the tensor.

        Args:
            c: Tensor of shape (N, 3, 3).

        Returns:
            The third invariant.
        """
        return np.linalg.det(c)

    @staticmethod
    def isochoric_invariant1(c: np.ndarray) -> np.ndarray:
        """Isochoric first invariant: tr(C) * det(C)^(-1/3)."""
        return np.einsum("nii->n", c) * np.linalg.det(c) ** (-1.0 / 3.0)  # type: ignore[no-any-return]

    @staticmethod
    def isochoric_invariant2(c: np.ndarray) -> np.ndarray:
        """Isochoric second invariant: 0.5*(I1^2 - tr(C^2)) * det(C)^(-2/3)."""
        i1 = np.einsum("nii->n", c)
        i1_sq = i1**2
        tr_c2 = np.einsum("nij,nji->n", c, c)
        return 0.5 * (i1_sq - tr_c2) * np.linalg.det(c) ** (-2.0 / 3.0)  # type: ignore[no-any-return]

    @staticmethod
    def right_cauchy_green(f: np.ndarray) -> Any:
        """
        Compute the right Cauchy-Green deformation tensor for a batch of deformation gradients
        using a more efficient vectorized approach.
        $$C = F^T F$$

        Args:
            f (np.ndarray): The deformation gradient tensor with shape (N, 3, 3),
                            where N is the number of deformation gradients.

        Returns:
            np.ndarray: The batch of right Cauchy-Green deformation tensors, shape (N, 3, 3).
        """
        # Use np.einsum to perform batch matrix multiplication: f's transpose @ f
        # The einsum subscript 'nij,nkj->nik' denotes batched matrix multiplication
        # where 'n' iterates over each matrix in the batch,
        # 'ji' are the indices of the transposed matrix,
        # and 'jk' are the indices for the second matrix.
        # Note: The difference from the left Cauchy-Green tensor is in the order of multiplication.
        return np.einsum("nji,njk->nik", f, f)

    @staticmethod
    def left_cauchy_green(f: np.ndarray) -> Any:
        """
        Compute the left Cauchy-Green deformation tensor for a batch of deformation gradients
        using a more efficient vectorized approach.

        Args:
            f (np.ndarray): The deformation gradient tensor with shape (N, 3, 3),
                            where N is the number of deformation gradients.

        Returns:
            np.ndarray: The batch of left Cauchy-Green deformation tensors, shape (N, 3, 3).
        """
        # Use np.einsum to perform batch matrix multiplication: f @ f's transpose
        # The einsum subscript 'nij,njk->nik' denotes batched matrix multiplication
        # where 'n' iterates over each matrix in the batch,
        # 'ij' are the indices of the first matrix,
        # and 'jk' are the indices for the second matrix (transposed to 'kj' for multiplication).
        return np.einsum("nij,nkj->nik", f, f)

    @staticmethod
    def pushforward(f: np.ndarray, tensor2D: np.ndarray) -> Any:
        """
        Forward tensor configuration.
        F*tensor2D*F^T. This is the forward transformation of a 2D tensor.

        Args:
            f (np.ndarray): deformation gradient # (N, 3, 3)
            tensor2D (np.ndarray): The 2D tensor to be mapped # (N, 3, 3)

        Returns:
            np.ndarray: The transformed tensor.
        """
        return np.einsum("nik,njl,nkl->nij", f, f, tensor2D)

    def principal_stretches(self, f: np.ndarray) -> np.ndarray:
        """
        Compute the principal stretches.

        Args:
            f (np.ndarray): The deformation gradient.

        Returns:
            np.ndarray: The principal stretches.
        """
        return np.sqrt(np.linalg.eigvals(self.right_cauchy_green(f)))

    def principal_directions(self, f: np.ndarray) -> np.ndarray:
        """
        Compute the principal directions.

        Args:
            f (np.ndarray): The deformation gradient.

        Returns:
            np.ndarray: The principal directions.
        """
        return np.linalg.eig(self.right_cauchy_green(f))[1]

    def right_stretch_tensor(self, f: np.ndarray) -> Any:
        """
        Compute the right stretch tensor.

        Args:
            f (np.ndarray): The deformation gradient.

        Returns:
            np.ndarray: The right stretch tensor.
        """
        v, vv = np.linalg.eig(self.right_cauchy_green(f))
        return np.einsum("...ij,...j->...ij", vv, v)

    def left_stretch_tensor(self, f: np.ndarray) -> Any:
        """
        Compute the left stretch tensor.

        Args:
            f (np.ndarray): The deformation gradient.

        Returns:
            np.ndarray: The left stretch tensor.
        """
        v, vv = np.linalg.eig(self.left_cauchy_green(f))
        return np.einsum("...ij,...j->...ij", vv, v)

    @staticmethod
    def fiber_invariant4(c: np.ndarray, a0: np.ndarray) -> np.ndarray:
        """Fiber invariant I4 = a0 . C . a0.

        Args:
            c: Right Cauchy-Green tensor (N, 3, 3).
            a0: Fiber direction (3,) or (N, 3).

        Returns:
            I4 values (N,).
        """
        a0 = np.asarray(a0)
        if a0.ndim == 1:
            return np.einsum("i,nij,j->n", a0, c, a0)  # type: ignore[no-any-return]
        return np.einsum("ni,nij,nj->n", a0, c, a0)  # type: ignore[no-any-return]

    @staticmethod
    def fiber_invariant5(c: np.ndarray, a0: np.ndarray) -> np.ndarray:
        """Fiber invariant I5 = a0 . C^2 . a0.

        Args:
            c: Right Cauchy-Green tensor (N, 3, 3).
            a0: Fiber direction (3,) or (N, 3).

        Returns:
            I5 values (N,).
        """
        a0 = np.asarray(a0)
        c2 = np.einsum("nij,njk->nik", c, c)
        if a0.ndim == 1:
            return np.einsum("i,nij,j->n", a0, c2, a0)  # type: ignore[no-any-return]
        return np.einsum("ni,nij,nj->n", a0, c2, a0)  # type: ignore[no-any-return]

    @staticmethod
    def fiber_invariants_multi(c: np.ndarray, fiber_directions: list[np.ndarray]) -> np.ndarray:
        """Compute I4, I5 for multiple fiber families.

        Args:
            c: Right Cauchy-Green tensor (N, 3, 3).
            fiber_directions: List of fiber direction vectors, each (3,).

        Returns:
            Array of shape (N, 2*num_fibers) with columns [I4_1, I5_1, I4_2, I5_2, ...].
        """
        cols: list[np.ndarray] = []
        for a0 in fiber_directions:
            cols.append(Kinematics.fiber_invariant4(c, a0))
            cols.append(Kinematics.fiber_invariant5(c, a0))
        return np.column_stack(cols)

    def rotation_tensor(self, f: np.ndarray) -> Any:
        """
        Compute the rotation tensors.

        Args:
            f (np.ndarray): The deformation gradients. batched with shape (N, 3, 3).

        Returns:
            np.ndarray: The rotation tensors. batched with shape (N, 3, 3).
        """
        return np.einsum("nij,njk->nik", f, np.linalg.inv(self.right_stretch_tensor(f)))

det_invariant(c) staticmethod

Calculate the third invariant (determinant) of the tensor.

Parameters:

Name	Type	Description	Default
c	ndarray	Tensor of shape (N, 3, 3).	required

Returns:

Type	Description
Any	The third invariant.

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def det_invariant(c: np.ndarray) -> Any:
    """
    Calculate the third invariant (determinant) of the tensor.

    Args:
        c: Tensor of shape (N, 3, 3).

    Returns:
        The third invariant.
    """
    return np.linalg.det(c)

fiber_invariant4(c, a0) staticmethod

Fiber invariant I4 = a0 . C . a0.

Parameters:

Name	Type	Description	Default
c	ndarray	Right Cauchy-Green tensor (N, 3, 3).	required
a0	ndarray	Fiber direction (3,) or (N, 3).	required

Returns:

Type	Description
ndarray	I4 values (N,).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def fiber_invariant4(c: np.ndarray, a0: np.ndarray) -> np.ndarray:
    """Fiber invariant I4 = a0 . C . a0.

    Args:
        c: Right Cauchy-Green tensor (N, 3, 3).
        a0: Fiber direction (3,) or (N, 3).

    Returns:
        I4 values (N,).
    """
    a0 = np.asarray(a0)
    if a0.ndim == 1:
        return np.einsum("i,nij,j->n", a0, c, a0)  # type: ignore[no-any-return]
    return np.einsum("ni,nij,nj->n", a0, c, a0)  # type: ignore[no-any-return]

fiber_invariant5(c, a0) staticmethod

Fiber invariant I5 = a0 . C^2 . a0.

Parameters:

Name	Type	Description	Default
c	ndarray	Right Cauchy-Green tensor (N, 3, 3).	required
a0	ndarray	Fiber direction (3,) or (N, 3).	required

Returns:

Type	Description
ndarray	I5 values (N,).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def fiber_invariant5(c: np.ndarray, a0: np.ndarray) -> np.ndarray:
    """Fiber invariant I5 = a0 . C^2 . a0.

    Args:
        c: Right Cauchy-Green tensor (N, 3, 3).
        a0: Fiber direction (3,) or (N, 3).

    Returns:
        I5 values (N,).
    """
    a0 = np.asarray(a0)
    c2 = np.einsum("nij,njk->nik", c, c)
    if a0.ndim == 1:
        return np.einsum("i,nij,j->n", a0, c2, a0)  # type: ignore[no-any-return]
    return np.einsum("ni,nij,nj->n", a0, c2, a0)  # type: ignore[no-any-return]

fiber_invariants_multi(c, fiber_directions) staticmethod

Compute I4, I5 for multiple fiber families.

Parameters:

Name	Type	Description	Default
c	ndarray	Right Cauchy-Green tensor (N, 3, 3).	required
fiber_directions	list[ndarray]	List of fiber direction vectors, each (3,).	required

Returns:

Type	Description
ndarray	Array of shape (N, 2*num_fibers) with columns [I4_1, I5_1, I4_2, I5_2, ...].

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def fiber_invariants_multi(c: np.ndarray, fiber_directions: list[np.ndarray]) -> np.ndarray:
    """Compute I4, I5 for multiple fiber families.

    Args:
        c: Right Cauchy-Green tensor (N, 3, 3).
        fiber_directions: List of fiber direction vectors, each (3,).

    Returns:
        Array of shape (N, 2*num_fibers) with columns [I4_1, I5_1, I4_2, I5_2, ...].
    """
    cols: list[np.ndarray] = []
    for a0 in fiber_directions:
        cols.append(Kinematics.fiber_invariant4(c, a0))
        cols.append(Kinematics.fiber_invariant5(c, a0))
    return np.column_stack(cols)

isochoric_invariant1(c) staticmethod

Isochoric first invariant: tr(C) * det(C)^(-⅓).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def isochoric_invariant1(c: np.ndarray) -> np.ndarray:
    """Isochoric first invariant: tr(C) * det(C)^(-1/3)."""
    return np.einsum("nii->n", c) * np.linalg.det(c) ** (-1.0 / 3.0)  # type: ignore[no-any-return]

isochoric_invariant2(c) staticmethod

Isochoric second invariant: 0.5*(I1^2 - tr(C^2)) * det(C)^(-⅔).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def isochoric_invariant2(c: np.ndarray) -> np.ndarray:
    """Isochoric second invariant: 0.5*(I1^2 - tr(C^2)) * det(C)^(-2/3)."""
    i1 = np.einsum("nii->n", c)
    i1_sq = i1**2
    tr_c2 = np.einsum("nij,nji->n", c, c)
    return 0.5 * (i1_sq - tr_c2) * np.linalg.det(c) ** (-2.0 / 3.0)  # type: ignore[no-any-return]

jacobian(f) staticmethod

Compute the Jacobian of the deformation gradient.

Parameters:

Name	Type	Description	Default
f	ndarray	4D tensor of shape (N, 3, 3).	required

Returns:

Type	Description
Any	np.ndarray: The Jacobian of the deformation gradient.

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def jacobian(f: np.ndarray) -> Any:
    """
    Compute the Jacobian of the deformation gradient.

    Args:
        f: 4D tensor of shape (N, 3, 3).

    Returns:
        np.ndarray: The Jacobian of the deformation gradient.
    """
    return np.linalg.det(f)

left_cauchy_green(f) staticmethod

Compute the left Cauchy-Green deformation tensor for a batch of deformation gradients using a more efficient vectorized approach.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient tensor with shape (N, 3, 3), where N is the number of deformation gradients.	required

Returns:

Type	Description
Any	np.ndarray: The batch of left Cauchy-Green deformation tensors, shape (N, 3, 3).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def left_cauchy_green(f: np.ndarray) -> Any:
    """
    Compute the left Cauchy-Green deformation tensor for a batch of deformation gradients
    using a more efficient vectorized approach.

    Args:
        f (np.ndarray): The deformation gradient tensor with shape (N, 3, 3),
                        where N is the number of deformation gradients.

    Returns:
        np.ndarray: The batch of left Cauchy-Green deformation tensors, shape (N, 3, 3).
    """
    # Use np.einsum to perform batch matrix multiplication: f @ f's transpose
    # The einsum subscript 'nij,njk->nik' denotes batched matrix multiplication
    # where 'n' iterates over each matrix in the batch,
    # 'ij' are the indices of the first matrix,
    # and 'jk' are the indices for the second matrix (transposed to 'kj' for multiplication).
    return np.einsum("nij,nkj->nik", f, f)

left_stretch_tensor(f)

Compute the left stretch tensor.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient.	required

Returns:

Type	Description
Any	np.ndarray: The left stretch tensor.

Source code in hyper_surrogate/mechanics/kinematics.py

def left_stretch_tensor(self, f: np.ndarray) -> Any:
    """
    Compute the left stretch tensor.

    Args:
        f (np.ndarray): The deformation gradient.

    Returns:
        np.ndarray: The left stretch tensor.
    """
    v, vv = np.linalg.eig(self.left_cauchy_green(f))
    return np.einsum("...ij,...j->...ij", vv, v)

principal_directions(f)

Compute the principal directions.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient.	required

Returns:

Type	Description
ndarray	np.ndarray: The principal directions.

Source code in hyper_surrogate/mechanics/kinematics.py

def principal_directions(self, f: np.ndarray) -> np.ndarray:
    """
    Compute the principal directions.

    Args:
        f (np.ndarray): The deformation gradient.

    Returns:
        np.ndarray: The principal directions.
    """
    return np.linalg.eig(self.right_cauchy_green(f))[1]

principal_stretches(f)

Compute the principal stretches.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient.	required

Returns:

Type	Description
ndarray	np.ndarray: The principal stretches.

Source code in hyper_surrogate/mechanics/kinematics.py

def principal_stretches(self, f: np.ndarray) -> np.ndarray:
    """
    Compute the principal stretches.

    Args:
        f (np.ndarray): The deformation gradient.

    Returns:
        np.ndarray: The principal stretches.
    """
    return np.sqrt(np.linalg.eigvals(self.right_cauchy_green(f)))

pushforward(f, tensor2D) staticmethod

Forward tensor configuration. Ftensor2DF^T. This is the forward transformation of a 2D tensor.

Parameters:

Name	Type	Description	Default
f	ndarray	deformation gradient # (N, 3, 3)	required
tensor2D	ndarray	The 2D tensor to be mapped # (N, 3, 3)	required

Returns:

Type	Description
Any	np.ndarray: The transformed tensor.

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def pushforward(f: np.ndarray, tensor2D: np.ndarray) -> Any:
    """
    Forward tensor configuration.
    F*tensor2D*F^T. This is the forward transformation of a 2D tensor.

    Args:
        f (np.ndarray): deformation gradient # (N, 3, 3)
        tensor2D (np.ndarray): The 2D tensor to be mapped # (N, 3, 3)

    Returns:
        np.ndarray: The transformed tensor.
    """
    return np.einsum("nik,njl,nkl->nij", f, f, tensor2D)

quadratic_invariant(c) staticmethod

Calculate the second invariant of the tensor.

Parameters:

Name	Type	Description	Default
c	ndarray	Tensor of shape (N, 3, 3).	required

Returns:

Type	Description
Any	The second invariant.

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def quadratic_invariant(c: np.ndarray) -> Any:
    """
    Calculate the second invariant of the tensor.

    Args:
        c: Tensor of shape (N, 3, 3).

    Returns:
        The second invariant.
    """
    return 0.5 * (np.einsum("nii->n", c) ** 2 - np.einsum("nij,nji->n", c, c))

right_cauchy_green(f) staticmethod

Compute the right Cauchy-Green deformation tensor for a batch of deformation gradients using a more efficient vectorized approach. $\(C = F^T F\)$

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient tensor with shape (N, 3, 3), where N is the number of deformation gradients.	required

Returns:

Type	Description
Any	np.ndarray: The batch of right Cauchy-Green deformation tensors, shape (N, 3, 3).

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def right_cauchy_green(f: np.ndarray) -> Any:
    """
    Compute the right Cauchy-Green deformation tensor for a batch of deformation gradients
    using a more efficient vectorized approach.
    $$C = F^T F$$

    Args:
        f (np.ndarray): The deformation gradient tensor with shape (N, 3, 3),
                        where N is the number of deformation gradients.

    Returns:
        np.ndarray: The batch of right Cauchy-Green deformation tensors, shape (N, 3, 3).
    """
    # Use np.einsum to perform batch matrix multiplication: f's transpose @ f
    # The einsum subscript 'nij,nkj->nik' denotes batched matrix multiplication
    # where 'n' iterates over each matrix in the batch,
    # 'ji' are the indices of the transposed matrix,
    # and 'jk' are the indices for the second matrix.
    # Note: The difference from the left Cauchy-Green tensor is in the order of multiplication.
    return np.einsum("nji,njk->nik", f, f)

right_stretch_tensor(f)

Compute the right stretch tensor.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradient.	required

Returns:

Type	Description
Any	np.ndarray: The right stretch tensor.

Source code in hyper_surrogate/mechanics/kinematics.py

def right_stretch_tensor(self, f: np.ndarray) -> Any:
    """
    Compute the right stretch tensor.

    Args:
        f (np.ndarray): The deformation gradient.

    Returns:
        np.ndarray: The right stretch tensor.
    """
    v, vv = np.linalg.eig(self.right_cauchy_green(f))
    return np.einsum("...ij,...j->...ij", vv, v)

rotation_tensor(f)

Compute the rotation tensors.

Parameters:

Name	Type	Description	Default
f	ndarray	The deformation gradients. batched with shape (N, 3, 3).	required

Returns:

Type	Description
Any	np.ndarray: The rotation tensors. batched with shape (N, 3, 3).

Source code in hyper_surrogate/mechanics/kinematics.py

def rotation_tensor(self, f: np.ndarray) -> Any:
    """
    Compute the rotation tensors.

    Args:
        f (np.ndarray): The deformation gradients. batched with shape (N, 3, 3).

    Returns:
        np.ndarray: The rotation tensors. batched with shape (N, 3, 3).
    """
    return np.einsum("nij,njk->nik", f, np.linalg.inv(self.right_stretch_tensor(f)))

trace_invariant(c) staticmethod

Calculate the first invariant (trace) of each tensor in the batch.

Parameters:

Name	Type	Description	Default
c	ndarray	Tensor of shape (N, 3, 3).	required

Returns:

Type	Description
Any	The first invariant of each tensor in the batch.

Source code in hyper_surrogate/mechanics/kinematics.py

@staticmethod
def trace_invariant(c: np.ndarray) -> Any:
    """
    Calculate the first invariant (trace) of each tensor in the batch.

    Args:
        c: Tensor of shape (N, 3, 3).

    Returns:
        The first invariant of each tensor in the batch.
    """
    return np.einsum("nii->n", c)

Demiray

Bases: Material

Demiray exponential hyperelastic model.

W = (C1 / C2) * (exp(C2 * (I1_bar - 3)) - 1) + vol(J)

Source code in hyper_surrogate/mechanics/materials.py

class Demiray(Material):
    """Demiray exponential hyperelastic model.

    W = (C1 / C2) * (exp(C2 * (I1_bar - 3)) - 1) + vol(J)
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "C1": 0.05,
        "C2": 8.0,
        "KBULK": 1000.0,
    }

    def __init__(self, parameters: dict[str, float] | None = None) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        super().__init__(params, fiber_direction=None)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        h = self._handler
        C1 = self._symbols["C1"]
        C2 = self._symbols["C2"]
        KBULK = self._symbols["KBULK"]
        return (C1 / C2) * (sympy_exp(C2 * (h.isochoric_invariant1 - 3)) - 1) + Rational(1, 4) * KBULK * (
            h.invariant3 - 1 - 2 * log(h.invariant3 ** Rational(1, 2))
        )

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        C1 = self._symbols["C1"]
        C2 = self._symbols["C2"]
        return (C1 / C2) * (sympy_exp(C2 * (i1_bar - 3)) - 1) + self._volumetric(j)

Fung

Bases: Material

Fung-type exponential model with configurable Q quadratic form.

W = (c / 2) * (exp(Q) - 1) + vol(J)

Q = sum_ij b_ij * E_ij * E_ij (Green-Lagrange strain components)

Default: isotropic Q = b1(E11^2 + E22^2 + E33^2) + b2(E12^2 + E13^2 + E23^2)

Source code in hyper_surrogate/mechanics/materials.py

class Fung(Material):
    """Fung-type exponential model with configurable Q quadratic form.

    W = (c / 2) * (exp(Q) - 1) + vol(J)

    Q = sum_ij b_ij * E_ij * E_ij  (Green-Lagrange strain components)

    Default: isotropic Q = b1*(E11^2 + E22^2 + E33^2) + b2*(E12^2 + E13^2 + E23^2)
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "c": 1.0,
        "b1": 10.0,
        "b2": 5.0,
        "KBULK": 1000.0,
    }

    def __init__(self, parameters: dict[str, float] | None = None) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        super().__init__(params, fiber_direction=None)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "Fung.sef is expressed in Green strain; use evaluate_energy instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        msg = "Fung model is expressed in Green strain components, not invariants"
        raise NotImplementedError(msg)

    def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
        """Evaluate Fung energy via Green strain."""
        from hyper_surrogate.mechanics.kinematics import Kinematics

        E = 0.5 * (c_batch - np.eye(3))

        c_p = self._params["c"]
        b1 = self._params["b1"]
        b2 = self._params["b2"]

        Q = b1 * (E[:, 0, 0] ** 2 + E[:, 1, 1] ** 2 + E[:, 2, 2] ** 2) + b2 * 2 * (
            E[:, 0, 1] ** 2 + E[:, 0, 2] ** 2 + E[:, 1, 2] ** 2
        )

        W = 0.5 * c_p * (np.exp(Q) - 1)

        # Volumetric
        j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
        KBULK = self._params["KBULK"]
        j2 = j_vals**2
        W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

        return W

    def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
        """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

        Returns (N, 6) array. Not to be confused with invariant gradients.
        """
        n = len(c_batch)
        eps = 1e-7
        W0 = self.evaluate_energy(c_batch)
        grad = np.zeros((n, 6))
        idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
        for k, (i, j) in enumerate(idx_map):
            c_p = c_batch.copy()
            c_p[:, i, j] += eps
            c_p[:, j, i] += eps
            W_p = self.evaluate_energy(c_p)
            grad[:, k] = (W_p - W0) / eps
        return grad

    def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
        """Not available for Fung — uses C-component formulation, not invariants."""
        msg = "Fung uses C-component formulation; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
        raise NotImplementedError(msg)

evaluate_energy(c_batch)

Evaluate Fung energy via Green strain.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
    """Evaluate Fung energy via Green strain."""
    from hyper_surrogate.mechanics.kinematics import Kinematics

    E = 0.5 * (c_batch - np.eye(3))

    c_p = self._params["c"]
    b1 = self._params["b1"]
    b2 = self._params["b2"]

    Q = b1 * (E[:, 0, 0] ** 2 + E[:, 1, 1] ** 2 + E[:, 2, 2] ** 2) + b2 * 2 * (
        E[:, 0, 1] ** 2 + E[:, 0, 2] ** 2 + E[:, 1, 2] ** 2
    )

    W = 0.5 * c_p * (np.exp(Q) - 1)

    # Volumetric
    j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
    KBULK = self._params["KBULK"]
    j2 = j_vals**2
    W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

    return W

evaluate_energy_grad_invariants(c_batch)

Not available for Fung — uses C-component formulation, not invariants.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
    """Not available for Fung — uses C-component formulation, not invariants."""
    msg = "Fung uses C-component formulation; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
    raise NotImplementedError(msg)

evaluate_energy_grad_voigt(c_batch)

Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

Returns (N, 6) array. Not to be confused with invariant gradients.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
    """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

    Returns (N, 6) array. Not to be confused with invariant gradients.
    """
    n = len(c_batch)
    eps = 1e-7
    W0 = self.evaluate_energy(c_batch)
    grad = np.zeros((n, 6))
    idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
    for k, (i, j) in enumerate(idx_map):
        c_p = c_batch.copy()
        c_p[:, i, j] += eps
        c_p[:, j, i] += eps
        W_p = self.evaluate_energy(c_p)
        grad[:, k] = (W_p - W0) / eps
    return grad

GasserOgdenHolzapfel

Bases: Material

Gasser-Ogden-Holzapfel model with fiber dispersion (kappa).

W = (a/2b)(exp(b(I1_bar - 3)) - 1) + (af/2bf)(exp(bf * E_bar^2) - 1) + vol(J)

where E_bar = kappa(I1_bar - 3) + (1 - 3kappa)*(I4 - 1).

kappa in [0, ⅓]: kappa=0 recovers HolzapfelOgden, kappa=⅓ gives isotropic fiber.

Source code in hyper_surrogate/mechanics/materials.py

class GasserOgdenHolzapfel(Material):
    """Gasser-Ogden-Holzapfel model with fiber dispersion (kappa).

    W = (a/2b)(exp(b(I1_bar - 3)) - 1)
      + (af/2bf)(exp(bf * E_bar^2) - 1)
      + vol(J)

    where E_bar = kappa*(I1_bar - 3) + (1 - 3*kappa)*(I4 - 1).

    kappa in [0, 1/3]: kappa=0 recovers HolzapfelOgden, kappa=1/3 gives isotropic fiber.
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "a": 0.059,
        "b": 8.023,
        "af": 18.472,
        "bf": 16.026,
        "kappa": 0.226,
        "KBULK": 1000.0,
    }

    def __init__(
        self,
        parameters: dict[str, float] | None = None,
        fiber_direction: np.ndarray | None = None,
    ) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        if fiber_direction is None:
            fiber_direction = np.array([1.0, 0.0, 0.0])
        super().__init__(params, fiber_direction=fiber_direction)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "GasserOgdenHolzapfel.sef requires fiber invariants; use sef_from_invariants instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        a_s = self._symbols["a"]
        b_s = self._symbols["b"]
        af_s = self._symbols["af"]
        bf_s = self._symbols["bf"]
        kappa_s = self._symbols["kappa"]

        # Isotropic ground substance
        W_iso = (a_s / (2 * b_s)) * (sympy_exp(b_s * (i1_bar - 3)) - 1)

        # Fiber contribution with dispersion
        if i4 is not None:
            E_bar = kappa_s * (i1_bar - 3) + (1 - 3 * kappa_s) * (i4 - 1)
            W_fiber: Expr = (af_s / (2 * bf_s)) * (sympy_exp(bf_s * E_bar**2) - 1)
        else:
            W_fiber = Symbol("zero") * 0

        return W_iso + W_fiber + self._volumetric(j)

Guccione

Bases: Material

Guccione transversely isotropic cardiac model.

W = (C / 2) * (exp(Q) - 1) + vol(J)

Q = bf * E_ff^2 + bt * (E_ss^2 + E_nn^2 + E_sn^2 + E_ns^2) + bfs * (E_fs^2 + E_sf^2 + E_fn^2 + E_nf^2)

where E = (C - I) / 2 is the Green-Lagrange strain in the fiber frame.

Source code in hyper_surrogate/mechanics/materials.py

class Guccione(Material):
    """Guccione transversely isotropic cardiac model.

    W = (C / 2) * (exp(Q) - 1) + vol(J)

    Q = bf * E_ff^2 + bt * (E_ss^2 + E_nn^2 + E_sn^2 + E_ns^2)
      + bfs * (E_fs^2 + E_sf^2 + E_fn^2 + E_nf^2)

    where E = (C - I) / 2 is the Green-Lagrange strain in the fiber frame.
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "C": 0.876,
        "bf": 18.48,
        "bt": 3.58,
        "bfs": 1.627,
        "KBULK": 1000.0,
    }

    def __init__(
        self,
        parameters: dict[str, float] | None = None,
        fiber_direction: np.ndarray | None = None,
        sheet_direction: np.ndarray | None = None,
    ) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        if fiber_direction is None:
            fiber_direction = np.array([1.0, 0.0, 0.0])
        self._sheet_direction = (
            np.asarray(sheet_direction, dtype=float) if sheet_direction is not None else np.array([0.0, 1.0, 0.0])
        )
        super().__init__(params, fiber_direction=fiber_direction)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "Guccione.sef requires fiber-frame computation; use evaluate_energy instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        msg = "Guccione model is expressed in Green strain components, not invariants"
        raise NotImplementedError(msg)

    def _rotation_matrix(self) -> np.ndarray:
        """Build rotation from global to fiber (f, s, n) frame."""
        fiber_dir = self.fiber_direction
        if fiber_dir is None:  # pragma: no cover
            msg = "Guccione requires a fiber direction"
            raise ValueError(msg)
        f0 = fiber_dir / np.linalg.norm(fiber_dir)
        s0 = self._sheet_direction / np.linalg.norm(self._sheet_direction)
        n0 = np.cross(f0, s0)
        n0 = n0 / np.linalg.norm(n0)
        return np.array([f0, s0, n0])  # (3, 3) rows = local axes

    def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
        """Evaluate Guccione energy via Green strain in fiber frame."""
        from hyper_surrogate.mechanics.kinematics import Kinematics

        Q_mat = self._rotation_matrix()  # (3, 3)

        # Green-Lagrange strain: E = (C - I) / 2
        E_global = 0.5 * (c_batch - np.eye(3))

        # Rotate to fiber frame: E_local = Q * E_global * Q^T
        E_local = np.einsum("ij,njk,lk->nil", Q_mat, E_global, Q_mat)

        C_p = self._params["C"]
        bf_p = self._params["bf"]
        bt_p = self._params["bt"]
        bfs_p = self._params["bfs"]

        E_ff = E_local[:, 0, 0]
        E_ss = E_local[:, 1, 1]
        E_nn = E_local[:, 2, 2]
        E_fs = E_local[:, 0, 1]
        E_fn = E_local[:, 0, 2]
        E_sn = E_local[:, 1, 2]

        Q_val = bf_p * E_ff**2 + bt_p * (E_ss**2 + E_nn**2 + 2 * E_sn**2) + bfs_p * (2 * E_fs**2 + 2 * E_fn**2)

        W = 0.5 * C_p * (np.exp(Q_val) - 1)

        # Volumetric
        j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
        KBULK = self._params["KBULK"]
        j2 = j_vals**2
        W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

        return np.asarray(W)

    def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
        """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

        Returns (N, 6) array. Not to be confused with invariant gradients.
        """
        n = len(c_batch)
        eps = 1e-7
        W0 = self.evaluate_energy(c_batch)
        grad = np.zeros((n, 6))
        idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
        for k, (i, j) in enumerate(idx_map):
            c_p = c_batch.copy()
            c_p[:, i, j] += eps
            c_p[:, j, i] += eps
            W_p = self.evaluate_energy(c_p)
            grad[:, k] = (W_p - W0) / eps
        return grad

    def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
        """Not available for Guccione — uses C-component formulation, not invariants."""
        msg = "Guccione uses C-component formulation; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
        raise NotImplementedError(msg)

evaluate_energy(c_batch)

Evaluate Guccione energy via Green strain in fiber frame.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
    """Evaluate Guccione energy via Green strain in fiber frame."""
    from hyper_surrogate.mechanics.kinematics import Kinematics

    Q_mat = self._rotation_matrix()  # (3, 3)

    # Green-Lagrange strain: E = (C - I) / 2
    E_global = 0.5 * (c_batch - np.eye(3))

    # Rotate to fiber frame: E_local = Q * E_global * Q^T
    E_local = np.einsum("ij,njk,lk->nil", Q_mat, E_global, Q_mat)

    C_p = self._params["C"]
    bf_p = self._params["bf"]
    bt_p = self._params["bt"]
    bfs_p = self._params["bfs"]

    E_ff = E_local[:, 0, 0]
    E_ss = E_local[:, 1, 1]
    E_nn = E_local[:, 2, 2]
    E_fs = E_local[:, 0, 1]
    E_fn = E_local[:, 0, 2]
    E_sn = E_local[:, 1, 2]

    Q_val = bf_p * E_ff**2 + bt_p * (E_ss**2 + E_nn**2 + 2 * E_sn**2) + bfs_p * (2 * E_fs**2 + 2 * E_fn**2)

    W = 0.5 * C_p * (np.exp(Q_val) - 1)

    # Volumetric
    j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
    KBULK = self._params["KBULK"]
    j2 = j_vals**2
    W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

    return np.asarray(W)

evaluate_energy_grad_invariants(c_batch)

Not available for Guccione — uses C-component formulation, not invariants.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
    """Not available for Guccione — uses C-component formulation, not invariants."""
    msg = "Guccione uses C-component formulation; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
    raise NotImplementedError(msg)

evaluate_energy_grad_voigt(c_batch)

Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

Returns (N, 6) array. Not to be confused with invariant gradients.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
    """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

    Returns (N, 6) array. Not to be confused with invariant gradients.
    """
    n = len(c_batch)
    eps = 1e-7
    W0 = self.evaluate_energy(c_batch)
    grad = np.zeros((n, 6))
    idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
    for k, (i, j) in enumerate(idx_map):
        c_p = c_batch.copy()
        c_p[:, i, j] += eps
        c_p[:, j, i] += eps
        W_p = self.evaluate_energy(c_p)
        grad[:, k] = (W_p - W0) / eps
    return grad

HolzapfelOgden

Bases: Material

Holzapfel-Ogden transversely isotropic hyperelastic model.

W = (a/2b)(exp(b(I1_bar - 3)) - 1) + (af/2bf)(exp(bf*^2) - 1) + vol(J)

where = max(x, 0) is the Macaulay bracket (fiber only under tension).

Source code in hyper_surrogate/mechanics/materials.py

class HolzapfelOgden(Material):
    """Holzapfel-Ogden transversely isotropic hyperelastic model.

    W = (a/2b)(exp(b(I1_bar - 3)) - 1) + (af/2bf)(exp(bf*<I4-1>^2) - 1) + vol(J)

    where <x> = max(x, 0) is the Macaulay bracket (fiber only under tension).
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "a": 0.059,
        "b": 8.023,
        "af": 18.472,
        "bf": 16.026,
        "KBULK": 1000.0,
    }

    def __init__(
        self,
        parameters: dict[str, float] | None = None,
        fiber_direction: np.ndarray | None = None,
    ) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        if fiber_direction is None:
            fiber_direction = np.array([1.0, 0.0, 0.0])
        super().__init__(params, fiber_direction=fiber_direction)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "HolzapfelOgden.sef requires fiber invariants; use sef_from_invariants instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        a_s = self._symbols["a"]
        b_s = self._symbols["b"]
        af_s = self._symbols["af"]
        bf_s = self._symbols["bf"]

        # Isotropic ground substance
        W_iso = (a_s / (2 * b_s)) * (sympy_exp(b_s * (i1_bar - 3)) - 1)

        # Fiber contribution (I4 term only; I5 omitted for simplicity)
        # Macaulay bracket: symbolic form; numerically clamped during evaluation
        W_fiber = (af_s / (2 * bf_s)) * (sympy_exp(bf_s * (i4 - 1) ** 2) - 1) if i4 is not None else Symbol("zero") * 0

        return W_iso + W_fiber + self._volumetric(j)

HolzapfelOgdenBiaxial

Bases: Material

Two-fiber Holzapfel-Ogden model for arterial wall tissue.

W = (a/2b)(exp(b(I1_bar - 3)) - 1) + sum_{k=1}^{2} (af/2bf)(exp(bf*(I4_k - 1)^2) - 1) + vol(J)

Each fiber family contributes an I4 term; both share the same parameters.

Source code in hyper_surrogate/mechanics/materials.py

class HolzapfelOgdenBiaxial(Material):
    """Two-fiber Holzapfel-Ogden model for arterial wall tissue.

    W = (a/2b)(exp(b(I1_bar - 3)) - 1)
      + sum_{k=1}^{2} (af/2bf)(exp(bf*(I4_k - 1)^2) - 1)
      + vol(J)

    Each fiber family contributes an I4 term; both share the same parameters.
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "a": 0.059,
        "b": 8.023,
        "af": 18.472,
        "bf": 16.026,
        "KBULK": 1000.0,
    }

    def __init__(
        self,
        parameters: dict[str, float] | None = None,
        fiber_directions: list[np.ndarray] | None = None,
    ) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        if fiber_directions is None:
            theta = np.radians(39.0)
            fiber_directions = [
                np.array([np.cos(theta), np.sin(theta), 0.0]),
                np.array([np.cos(theta), -np.sin(theta), 0.0]),
            ]
        if len(fiber_directions) != 2:
            msg = f"HolzapfelOgdenBiaxial requires exactly 2 fiber directions, got {len(fiber_directions)}"
            raise ValueError(msg)
        super().__init__(params, fiber_directions=fiber_directions)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "HolzapfelOgdenBiaxial.sef requires fiber invariants; use sef_from_all_invariants instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        msg = "HolzapfelOgdenBiaxial has 2 fibers; use sef_from_all_invariants instead"
        raise NotImplementedError(msg)

    def sef_from_all_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        fiber_invariants: list[tuple[Symbol, Symbol]] | None = None,
    ) -> Expr:
        a_s = self._symbols["a"]
        b_s = self._symbols["b"]
        af_s = self._symbols["af"]
        bf_s = self._symbols["bf"]

        # Isotropic ground substance
        W: Expr = (a_s / (2 * b_s)) * (sympy_exp(b_s * (i1_bar - 3)) - 1)

        # Fiber contributions
        if fiber_invariants is not None:
            for i4_k, _i5_k in fiber_invariants:
                W = W + (af_s / (2 * bf_s)) * (sympy_exp(bf_s * (i4_k - 1) ** 2) - 1)

        return W + self._volumetric(j)

    def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
        """Evaluate strain energy for (N,3,3) C tensors via invariants."""
        from sympy import lambdify as sym_lambdify

        from hyper_surrogate.mechanics.kinematics import Kinematics

        i1s, i2s, js = Symbol("I1b"), Symbol("I2b"), Symbol("J")
        fiber_inv_pairs: list[tuple[Symbol, Symbol]] = []
        inv_syms: list[Symbol] = [i1s, i2s, js]
        for k in range(self.num_fiber_families):
            i4k, i5k = Symbol(f"I4_{k}"), Symbol(f"I5_{k}")
            inv_syms.extend([i4k, i5k])
            fiber_inv_pairs.append((i4k, i5k))

        W = self.sef_from_all_invariants(i1s, i2s, js, fiber_inv_pairs)
        param_syms = list(self._symbols.values())
        fn = sym_lambdify((*inv_syms, *param_syms), W, modules="numpy")

        i1 = Kinematics.isochoric_invariant1(c_batch)
        i2 = Kinematics.isochoric_invariant2(c_batch)
        j = np.sqrt(Kinematics.det_invariant(c_batch))

        inv_vals: list[np.ndarray] = [i1, i2, j]
        for a0 in self._fiber_directions:
            inv_vals.append(Kinematics.fiber_invariant4(c_batch, a0))
            inv_vals.append(Kinematics.fiber_invariant5(c_batch, a0))

        param_vals = list(self._params.values())
        results = fn(*inv_vals, *param_vals)
        return np.asarray(results, dtype=float)

evaluate_energy(c_batch)

Evaluate strain energy for (N,3,3) C tensors via invariants.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
    """Evaluate strain energy for (N,3,3) C tensors via invariants."""
    from sympy import lambdify as sym_lambdify

    from hyper_surrogate.mechanics.kinematics import Kinematics

    i1s, i2s, js = Symbol("I1b"), Symbol("I2b"), Symbol("J")
    fiber_inv_pairs: list[tuple[Symbol, Symbol]] = []
    inv_syms: list[Symbol] = [i1s, i2s, js]
    for k in range(self.num_fiber_families):
        i4k, i5k = Symbol(f"I4_{k}"), Symbol(f"I5_{k}")
        inv_syms.extend([i4k, i5k])
        fiber_inv_pairs.append((i4k, i5k))

    W = self.sef_from_all_invariants(i1s, i2s, js, fiber_inv_pairs)
    param_syms = list(self._symbols.values())
    fn = sym_lambdify((*inv_syms, *param_syms), W, modules="numpy")

    i1 = Kinematics.isochoric_invariant1(c_batch)
    i2 = Kinematics.isochoric_invariant2(c_batch)
    j = np.sqrt(Kinematics.det_invariant(c_batch))

    inv_vals: list[np.ndarray] = [i1, i2, j]
    for a0 in self._fiber_directions:
        inv_vals.append(Kinematics.fiber_invariant4(c_batch, a0))
        inv_vals.append(Kinematics.fiber_invariant5(c_batch, a0))

    param_vals = list(self._params.values())
    results = fn(*inv_vals, *param_vals)
    return np.asarray(results, dtype=float)

Material

Base class for constitutive material models using composition.

Source code in hyper_surrogate/mechanics/materials.py

class Material:
    """Base class for constitutive material models using composition."""

    def __init__(
        self,
        parameters: dict[str, float],
        fiber_direction: np.ndarray | None = None,
        fiber_directions: list[np.ndarray] | None = None,
    ) -> None:
        self._handler = SymbolicHandler()
        self._params = parameters
        self._symbols = {k: Symbol(k) for k in parameters}
        # Normalize fiber storage: always a list
        if fiber_directions is not None:
            self._fiber_directions = [np.asarray(d, dtype=float) for d in fiber_directions]
        elif fiber_direction is not None:
            self._fiber_directions = [np.asarray(fiber_direction, dtype=float)]
        else:
            self._fiber_directions = []

    @property
    def is_anisotropic(self) -> bool:
        """Whether this material depends on fiber invariants (I4, I5)."""
        return len(self._fiber_directions) > 0

    @property
    def fiber_direction(self) -> np.ndarray | None:
        """First fiber direction (backward compat). None if isotropic."""
        return self._fiber_directions[0] if self._fiber_directions else None

    @property
    def fiber_directions(self) -> list[np.ndarray]:
        """All fiber directions. Empty list if isotropic."""
        return self._fiber_directions

    @property
    def num_fiber_families(self) -> int:
        """Number of fiber families."""
        return len(self._fiber_directions)

    @property
    def input_dim(self) -> int:
        """Invariant input dimension: 3 (isotropic) + 2 per fiber family."""
        return 3 + 2 * self.num_fiber_families

    @property
    def handler(self) -> SymbolicHandler:
        return self._handler

    @property
    def sef(self) -> Expr:
        raise NotImplementedError

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        """Return SEF as a function of invariant symbols (I1_bar, I2_bar, J[, I4, I5]).

        Override in subclasses to enable energy gradient computation w.r.t. invariants.
        """
        raise NotImplementedError

    def sef_from_all_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        fiber_invariants: list[tuple[Symbol, Symbol]] | None = None,
    ) -> Expr:
        """Return SEF as a function of all invariant symbols including multi-fiber.

        ``fiber_invariants`` is a list of (I4_k, I5_k) tuples, one per fiber family.
        Default implementation delegates to :meth:`sef_from_invariants` for 0 or 1 fibers.
        Override in multi-fiber subclasses.
        """
        if not fiber_invariants:
            return self.sef_from_invariants(i1_bar, i2_bar, j)
        if len(fiber_invariants) == 1:
            i4, i5 = fiber_invariants[0]
            return self.sef_from_invariants(i1_bar, i2_bar, j, i4, i5)
        msg = f"{type(self).__name__} does not support {len(fiber_invariants)} fiber families"
        raise NotImplementedError(msg)

    @cached_property
    def pk2_expr(self) -> Matrix:
        return self._handler.pk2(self.sef)

    @cached_property
    def cmat_expr(self) -> Any:
        return self._handler.cmat(self.pk2_expr)

    @cached_property
    def pk2_func(self) -> Callable:
        return self._handler.lambdify(self.pk2_expr, *self._symbols.values())  # type: ignore[no-any-return]

    @cached_property
    def cmat_func(self) -> Callable:
        return self._handler.lambdify(self.cmat_expr, *self._symbols.values())  # type: ignore[no-any-return]

    def evaluate_pk2(self, c_batch: np.ndarray) -> np.ndarray:
        """Vectorized PK2 evaluation over (N,3,3) C tensors."""
        param_values = list(self._params.values())
        results = []
        for c in c_batch:
            result = self.pk2_func(c.flatten(), *param_values)
            results.append(np.array(result, dtype=float))
        return np.array(results)

    def evaluate_cmat(self, c_batch: np.ndarray) -> np.ndarray:
        """Vectorized CMAT evaluation over (N,3,3) C tensors."""
        param_values = list(self._params.values())
        results = []
        for c in c_batch:
            result = self.cmat_func(c.flatten(), *param_values)
            results.append(np.array(result, dtype=float))
        return np.array(results)

    def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
        """Compute dW/d(invariants) for (N,3,3) C tensors.

        Returns (N, input_dim) where input_dim = 3 + 2*num_fiber_families.
        """
        from sympy import diff
        from sympy import lambdify as sym_lambdify

        from hyper_surrogate.mechanics.kinematics import Kinematics

        i1s, i2s, js = Symbol("I1b"), Symbol("I2b"), Symbol("J")

        inv_syms: list[Symbol] = [i1s, i2s, js]
        fiber_inv_pairs: list[tuple[Symbol, Symbol]] = []

        for k in range(self.num_fiber_families):
            i4k = Symbol(f"I4_{k}")
            i5k = Symbol(f"I5_{k}")
            inv_syms.extend([i4k, i5k])
            fiber_inv_pairs.append((i4k, i5k))

        W = self.sef_from_all_invariants(i1s, i2s, js, fiber_inv_pairs or None)

        dW = [diff(W, s) for s in inv_syms]
        param_syms = list(self._symbols.values())
        fn = sym_lambdify((*inv_syms, *param_syms), dW, modules="numpy")

        i1 = Kinematics.isochoric_invariant1(c_batch)
        i2 = Kinematics.isochoric_invariant2(c_batch)
        j = np.sqrt(Kinematics.det_invariant(c_batch))

        param_vals = list(self._params.values())
        n = len(c_batch)

        inv_vals: list[np.ndarray] = [i1, i2, j]
        for a0 in self._fiber_directions:
            inv_vals.append(Kinematics.fiber_invariant4(c_batch, a0))
            inv_vals.append(Kinematics.fiber_invariant5(c_batch, a0))

        results = fn(*inv_vals, *param_vals)

        return np.column_stack([np.broadcast_to(np.asarray(r, dtype=float), (n,)) for r in results])

    def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
        """Evaluate strain energy for (N,3,3) C tensors. Returns (N,)."""
        from sympy import lambdify as sym_lambdify

        c_syms = self._handler.c_symbols()
        sef_func = sym_lambdify((c_syms, *self._symbols.values()), self.sef, modules="numpy")
        param_values = list(self._params.values())
        results = []
        for c in c_batch:
            result = sef_func(c.flatten(), *param_values)
            results.append(float(result))
        return np.array(results)

    # --- Symbolic accessors for UMAT generation ---

    def cauchy_voigt(self, f: Matrix) -> Matrix:
        """Voigt-reduced Cauchy stress (6x1) in symbolic form."""
        sigma = self._handler.cauchy(self.sef, f)
        return SymbolicHandler.to_voigt_2(sigma)

    def tangent_voigt(self, f: Matrix, use_jaumann_rate: bool = False) -> Matrix:
        """Voigt-reduced tangent (6x6) in symbolic form."""
        smat = self._handler.spatial_tangent(self.pk2_expr, f)
        if use_jaumann_rate:
            sigma = self._handler.cauchy(self.sef, f)
            smat = smat + self._handler.jaumann_correction(sigma)
        return SymbolicHandler.to_voigt_4(smat)

fiber_direction property

First fiber direction (backward compat). None if isotropic.

fiber_directions property

All fiber directions. Empty list if isotropic.

input_dim property

Invariant input dimension: 3 (isotropic) + 2 per fiber family.

is_anisotropic property

Whether this material depends on fiber invariants (I4, I5).

num_fiber_families property

Number of fiber families.

cauchy_voigt(f)

Voigt-reduced Cauchy stress (6x1) in symbolic form.

Source code in hyper_surrogate/mechanics/materials.py

def cauchy_voigt(self, f: Matrix) -> Matrix:
    """Voigt-reduced Cauchy stress (6x1) in symbolic form."""
    sigma = self._handler.cauchy(self.sef, f)
    return SymbolicHandler.to_voigt_2(sigma)

evaluate_cmat(c_batch)

Vectorized CMAT evaluation over (N,3,3) C tensors.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_cmat(self, c_batch: np.ndarray) -> np.ndarray:
    """Vectorized CMAT evaluation over (N,3,3) C tensors."""
    param_values = list(self._params.values())
    results = []
    for c in c_batch:
        result = self.cmat_func(c.flatten(), *param_values)
        results.append(np.array(result, dtype=float))
    return np.array(results)

evaluate_energy(c_batch)

Evaluate strain energy for (N,3,3) C tensors. Returns (N,).

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
    """Evaluate strain energy for (N,3,3) C tensors. Returns (N,)."""
    from sympy import lambdify as sym_lambdify

    c_syms = self._handler.c_symbols()
    sef_func = sym_lambdify((c_syms, *self._symbols.values()), self.sef, modules="numpy")
    param_values = list(self._params.values())
    results = []
    for c in c_batch:
        result = sef_func(c.flatten(), *param_values)
        results.append(float(result))
    return np.array(results)

evaluate_energy_grad_invariants(c_batch)

Compute dW/d(invariants) for (N,3,3) C tensors.

Returns (N, input_dim) where input_dim = 3 + 2*num_fiber_families.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
    """Compute dW/d(invariants) for (N,3,3) C tensors.

    Returns (N, input_dim) where input_dim = 3 + 2*num_fiber_families.
    """
    from sympy import diff
    from sympy import lambdify as sym_lambdify

    from hyper_surrogate.mechanics.kinematics import Kinematics

    i1s, i2s, js = Symbol("I1b"), Symbol("I2b"), Symbol("J")

    inv_syms: list[Symbol] = [i1s, i2s, js]
    fiber_inv_pairs: list[tuple[Symbol, Symbol]] = []

    for k in range(self.num_fiber_families):
        i4k = Symbol(f"I4_{k}")
        i5k = Symbol(f"I5_{k}")
        inv_syms.extend([i4k, i5k])
        fiber_inv_pairs.append((i4k, i5k))

    W = self.sef_from_all_invariants(i1s, i2s, js, fiber_inv_pairs or None)

    dW = [diff(W, s) for s in inv_syms]
    param_syms = list(self._symbols.values())
    fn = sym_lambdify((*inv_syms, *param_syms), dW, modules="numpy")

    i1 = Kinematics.isochoric_invariant1(c_batch)
    i2 = Kinematics.isochoric_invariant2(c_batch)
    j = np.sqrt(Kinematics.det_invariant(c_batch))

    param_vals = list(self._params.values())
    n = len(c_batch)

    inv_vals: list[np.ndarray] = [i1, i2, j]
    for a0 in self._fiber_directions:
        inv_vals.append(Kinematics.fiber_invariant4(c_batch, a0))
        inv_vals.append(Kinematics.fiber_invariant5(c_batch, a0))

    results = fn(*inv_vals, *param_vals)

    return np.column_stack([np.broadcast_to(np.asarray(r, dtype=float), (n,)) for r in results])

evaluate_pk2(c_batch)

Vectorized PK2 evaluation over (N,3,3) C tensors.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_pk2(self, c_batch: np.ndarray) -> np.ndarray:
    """Vectorized PK2 evaluation over (N,3,3) C tensors."""
    param_values = list(self._params.values())
    results = []
    for c in c_batch:
        result = self.pk2_func(c.flatten(), *param_values)
        results.append(np.array(result, dtype=float))
    return np.array(results)

sef_from_all_invariants(i1_bar, i2_bar, j, fiber_invariants=None)

Return SEF as a function of all invariant symbols including multi-fiber.

fiber_invariants is a list of (I4_k, I5_k) tuples, one per fiber family. Default implementation delegates to :meth:sef_from_invariants for 0 or 1 fibers. Override in multi-fiber subclasses.

Source code in hyper_surrogate/mechanics/materials.py

def sef_from_all_invariants(
    self,
    i1_bar: Symbol,
    i2_bar: Symbol,
    j: Symbol,
    fiber_invariants: list[tuple[Symbol, Symbol]] | None = None,
) -> Expr:
    """Return SEF as a function of all invariant symbols including multi-fiber.

    ``fiber_invariants`` is a list of (I4_k, I5_k) tuples, one per fiber family.
    Default implementation delegates to :meth:`sef_from_invariants` for 0 or 1 fibers.
    Override in multi-fiber subclasses.
    """
    if not fiber_invariants:
        return self.sef_from_invariants(i1_bar, i2_bar, j)
    if len(fiber_invariants) == 1:
        i4, i5 = fiber_invariants[0]
        return self.sef_from_invariants(i1_bar, i2_bar, j, i4, i5)
    msg = f"{type(self).__name__} does not support {len(fiber_invariants)} fiber families"
    raise NotImplementedError(msg)

sef_from_invariants(i1_bar, i2_bar, j, i4=None, i5=None)

Return SEF as a function of invariant symbols (I1_bar, I2_bar, J[, I4, I5]).

Override in subclasses to enable energy gradient computation w.r.t. invariants.

Source code in hyper_surrogate/mechanics/materials.py

def sef_from_invariants(
    self,
    i1_bar: Symbol,
    i2_bar: Symbol,
    j: Symbol,
    i4: Symbol | None = None,
    i5: Symbol | None = None,
) -> Expr:
    """Return SEF as a function of invariant symbols (I1_bar, I2_bar, J[, I4, I5]).

    Override in subclasses to enable energy gradient computation w.r.t. invariants.
    """
    raise NotImplementedError

tangent_voigt(f, use_jaumann_rate=False)

Voigt-reduced tangent (6x6) in symbolic form.

Source code in hyper_surrogate/mechanics/materials.py

def tangent_voigt(self, f: Matrix, use_jaumann_rate: bool = False) -> Matrix:
    """Voigt-reduced tangent (6x6) in symbolic form."""
    smat = self._handler.spatial_tangent(self.pk2_expr, f)
    if use_jaumann_rate:
        sigma = self._handler.cauchy(self.sef, f)
        smat = smat + self._handler.jaumann_correction(sigma)
    return SymbolicHandler.to_voigt_4(smat)

Ogden

Bases: Material

Ogden hyperelastic model (N-term, principal stretch formulation).

W = sum_p (mu_p / alpha_p) * (lambda1^alpha_p + lambda2^alpha_p + lambda3^alpha_p - 3) + vol(J)

This model is expressed in principal stretches, not invariants. The sef_from_invariants method approximates via I1_bar and I2_bar using the relationship between invariants and symmetric functions of stretches.

Source code in hyper_surrogate/mechanics/materials.py

class Ogden(Material):
    """Ogden hyperelastic model (N-term, principal stretch formulation).

    W = sum_p (mu_p / alpha_p) * (lambda1^alpha_p + lambda2^alpha_p + lambda3^alpha_p - 3) + vol(J)

    This model is expressed in principal stretches, not invariants.
    The sef_from_invariants method approximates via I1_bar and I2_bar using
    the relationship between invariants and symmetric functions of stretches.
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "mu1": 1.0,
        "alpha1": 2.0,
        "KBULK": 1000.0,
    }

    def __init__(self, parameters: dict[str, float] | None = None) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        # Determine number of Ogden terms from parameter keys
        n = 1
        while f"mu{n + 1}" in params and f"alpha{n + 1}" in params:
            n += 1
        self._n_terms = n
        super().__init__(params, fiber_direction=None)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        msg = "Ogden.sef requires principal stretches; use evaluate_energy instead"
        raise NotImplementedError(msg)

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        msg = "Ogden model is not naturally expressed in invariants; use evaluate_energy instead"
        raise NotImplementedError(msg)

    def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
        """Evaluate Ogden energy via eigenvalues of C."""
        from hyper_surrogate.mechanics.kinematics import Kinematics

        n_samples = len(c_batch)
        j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
        # Principal stretches squared = eigenvalues of C
        eigvals = np.linalg.eigvalsh(c_batch)  # (N, 3), sorted ascending
        # Isochoric principal stretches: lambda_bar_i = J^(-1/3) * lambda_i
        lambda_bar_sq = eigvals * (j_vals ** (-2.0 / 3.0))[:, None]
        lambda_bar = np.sqrt(np.maximum(lambda_bar_sq, 1e-30))

        W = np.zeros(n_samples)
        for p in range(1, self._n_terms + 1):
            mu_p = self._params[f"mu{p}"]
            alpha_p = self._params[f"alpha{p}"]
            W += (mu_p / alpha_p) * (
                lambda_bar[:, 0] ** alpha_p + lambda_bar[:, 1] ** alpha_p + lambda_bar[:, 2] ** alpha_p - 3
            )

        # Volumetric
        KBULK = self._params["KBULK"]
        j2 = j_vals**2
        W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

        return W

    def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
        """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

        Returns (N, 6) array. Not to be confused with invariant gradients.
        """
        n = len(c_batch)
        eps = 1e-7
        W0 = self.evaluate_energy(c_batch)
        grad = np.zeros((n, 6))
        idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
        for k, (i, j) in enumerate(idx_map):
            c_p = c_batch.copy()
            c_p[:, i, j] += eps
            c_p[:, j, i] += eps  # symmetry
            W_p = self.evaluate_energy(c_p)
            grad[:, k] = (W_p - W0) / eps
        return grad

    def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
        """Not available for Ogden — uses principal stretch formulation, not invariants."""
        msg = "Ogden uses principal stretches; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
        raise NotImplementedError(msg)

evaluate_energy(c_batch)

Evaluate Ogden energy via eigenvalues of C.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy(self, c_batch: np.ndarray) -> np.ndarray:
    """Evaluate Ogden energy via eigenvalues of C."""
    from hyper_surrogate.mechanics.kinematics import Kinematics

    n_samples = len(c_batch)
    j_vals = np.sqrt(Kinematics.det_invariant(c_batch))
    # Principal stretches squared = eigenvalues of C
    eigvals = np.linalg.eigvalsh(c_batch)  # (N, 3), sorted ascending
    # Isochoric principal stretches: lambda_bar_i = J^(-1/3) * lambda_i
    lambda_bar_sq = eigvals * (j_vals ** (-2.0 / 3.0))[:, None]
    lambda_bar = np.sqrt(np.maximum(lambda_bar_sq, 1e-30))

    W = np.zeros(n_samples)
    for p in range(1, self._n_terms + 1):
        mu_p = self._params[f"mu{p}"]
        alpha_p = self._params[f"alpha{p}"]
        W += (mu_p / alpha_p) * (
            lambda_bar[:, 0] ** alpha_p + lambda_bar[:, 1] ** alpha_p + lambda_bar[:, 2] ** alpha_p - 3
        )

    # Volumetric
    KBULK = self._params["KBULK"]
    j2 = j_vals**2
    W += 0.25 * KBULK * (j2 - 1 - 2 * np.log(j_vals))

    return W

evaluate_energy_grad_invariants(c_batch)

Not available for Ogden — uses principal stretch formulation, not invariants.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_invariants(self, c_batch: np.ndarray) -> np.ndarray:
    """Not available for Ogden — uses principal stretch formulation, not invariants."""
    msg = "Ogden uses principal stretches; use evaluate_energy_grad_voigt for dW/dC Voigt gradients"
    raise NotImplementedError(msg)

evaluate_energy_grad_voigt(c_batch)

Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

Returns (N, 6) array. Not to be confused with invariant gradients.

Source code in hyper_surrogate/mechanics/materials.py

def evaluate_energy_grad_voigt(self, c_batch: np.ndarray) -> np.ndarray:
    """Numerical gradient dW/dC via finite differences on the 6 Voigt components of C.

    Returns (N, 6) array. Not to be confused with invariant gradients.
    """
    n = len(c_batch)
    eps = 1e-7
    W0 = self.evaluate_energy(c_batch)
    grad = np.zeros((n, 6))
    idx_map = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]
    for k, (i, j) in enumerate(idx_map):
        c_p = c_batch.copy()
        c_p[:, i, j] += eps
        c_p[:, j, i] += eps  # symmetry
        W_p = self.evaluate_energy(c_p)
        grad[:, k] = (W_p - W0) / eps
    return grad

Yeoh

Bases: Material

Yeoh hyperelastic model (third-order reduced polynomial).

W = C10(I1_bar - 3) + C20(I1_bar - 3)^2 + C30*(I1_bar - 3)^3 + vol(J)

Source code in hyper_surrogate/mechanics/materials.py

class Yeoh(Material):
    """Yeoh hyperelastic model (third-order reduced polynomial).

    W = C10*(I1_bar - 3) + C20*(I1_bar - 3)^2 + C30*(I1_bar - 3)^3 + vol(J)
    """

    DEFAULT_PARAMS: ClassVar[dict[str, float]] = {
        "C10": 0.5,
        "C20": -0.01,
        "C30": 0.001,
        "KBULK": 1000.0,
    }

    def __init__(self, parameters: dict[str, float] | None = None) -> None:
        params = {**self.DEFAULT_PARAMS, **(parameters or {})}
        super().__init__(params, fiber_direction=None)

    def _volumetric(self, j: Symbol) -> Expr:
        KBULK = self._symbols["KBULK"]
        i3 = j**2
        return Rational(1, 4) * KBULK * (i3 - 1 - 2 * log(i3 ** Rational(1, 2)))

    @property
    def sef(self) -> Expr:
        h = self._handler
        C10 = self._symbols["C10"]
        C20 = self._symbols["C20"]
        C30 = self._symbols["C30"]
        KBULK = self._symbols["KBULK"]
        i1_dev = h.isochoric_invariant1 - 3
        return (
            C10 * i1_dev
            + C20 * i1_dev**2
            + C30 * i1_dev**3
            + Rational(1, 4) * KBULK * (h.invariant3 - 1 - 2 * log(h.invariant3 ** Rational(1, 2)))
        )

    def sef_from_invariants(
        self,
        i1_bar: Symbol,
        i2_bar: Symbol,
        j: Symbol,
        i4: Symbol | None = None,
        i5: Symbol | None = None,
    ) -> Expr:
        C10 = self._symbols["C10"]
        C20 = self._symbols["C20"]
        C30 = self._symbols["C30"]
        i1_dev = i1_bar - 3
        return C10 * i1_dev + C20 * i1_dev**2 + C30 * i1_dev**3 + self._volumetric(j)

Data

DeformationGenerator

Generates physically valid deformation gradients for training data.

Source code in hyper_surrogate/data/deformation.py

class DeformationGenerator:
    """Generates physically valid deformation gradients for training data."""

    def __init__(self, seed: int = 42) -> None:
        self._rng = np.random.default_rng(seed)

    def uniaxial(self, n: int, stretch_range: tuple[float, float] = (0.4, 3.0)) -> np.ndarray:
        stretch = self._rng.uniform(*stretch_range, size=n)
        stretch_t = stretch**-0.5
        result = np.zeros((n, 3, 3))
        result[:, 0, 0] = stretch
        result[:, 1, 1] = stretch_t
        result[:, 2, 2] = stretch_t
        return result

    def biaxial(self, n: int, stretch_range: tuple[float, float] = (0.4, 3.0)) -> np.ndarray:
        s1 = self._rng.uniform(*stretch_range, size=n)
        s2 = self._rng.uniform(*stretch_range, size=n)
        s3 = (s1 * s2) ** -1.0
        result = np.zeros((n, 3, 3))
        result[:, 0, 0] = s1
        result[:, 1, 1] = s2
        result[:, 2, 2] = s3
        return result

    def shear(self, n: int, shear_range: tuple[float, float] = (-1.0, 1.0)) -> np.ndarray:
        gamma = self._rng.uniform(*shear_range, size=n)
        result = np.repeat(np.eye(3)[np.newaxis, :, :], n, axis=0)
        result[:, 0, 1] = gamma
        return result

    def random_rotation(self, n: int) -> np.ndarray:
        axes = self._rng.integers(0, 3, size=n)
        angles = self._rng.uniform(0, np.pi, size=n)
        rotations = []
        for ax, ang in zip(axes, angles, strict=False):
            c, s = np.cos(ang), np.sin(ang)
            if ax == 0:
                R = np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
            elif ax == 1:
                R = np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]])
            else:
                R = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
            rotations.append(R)
        return np.array(rotations)

    @staticmethod
    def _rotate(F: np.ndarray, R: np.ndarray) -> np.ndarray:
        return np.einsum("nij,njk,nlk->nil", R, F, R)  # type: ignore[no-any-return]

    def fiber_directions(
        self,
        n: int,
        preferred: np.ndarray | None = None,
        dispersion: float = 0.0,
    ) -> np.ndarray:
        """Generate fiber direction vectors.

        Args:
            n: Number of samples.
            preferred: Preferred fiber direction (3,). Default: [1, 0, 0].
            dispersion: Cone half-angle in radians. 0 = all aligned.

        Returns:
            Fiber directions (N, 3), unit vectors.
        """
        if preferred is None:
            preferred = np.array([1.0, 0.0, 0.0])
        preferred = preferred / np.linalg.norm(preferred)

        if dispersion <= 0.0:
            return np.tile(preferred, (n, 1))

        # Sample directions in a cone around preferred
        # Build local frame: preferred = e3, find orthogonal e1, e2
        if abs(preferred[2]) < 0.9:
            t = np.cross(preferred, np.array([0.0, 0.0, 1.0]))
        else:
            t = np.cross(preferred, np.array([1.0, 0.0, 0.0]))
        e1 = t / np.linalg.norm(t)
        e2 = np.cross(preferred, e1)

        phi = self._rng.uniform(0, 2 * np.pi, size=n)
        theta = self._rng.uniform(0, dispersion, size=n)

        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)

        dirs = x[:, None] * e1 + y[:, None] * e2 + z[:, None] * preferred
        norms = np.linalg.norm(dirs, axis=1, keepdims=True)
        result: np.ndarray = dirs / norms
        return result

    def combined(
        self,
        n: int,
        stretch_range: tuple[float, float] = (0.4, 3.0),
        shear_range: tuple[float, float] = (-1.0, 1.0),
    ) -> np.ndarray:
        fu = self.uniaxial(n, stretch_range)
        fs = self.shear(n, shear_range)
        fb = self.biaxial(n, stretch_range)
        r1, r2, r3 = self.random_rotation(n), self.random_rotation(n), self.random_rotation(n)
        fu = self._rotate(fu, r1)
        fs = self._rotate(fs, r2)
        fb = self._rotate(fb, r3)
        return np.matmul(np.matmul(fb, fu), fs)  # type: ignore[no-any-return]

fiber_directions(n, preferred=None, dispersion=0.0)

Generate fiber direction vectors.

Parameters:

Name	Type	Description	Default
n	int	Number of samples.	required
preferred	ndarray | None	Preferred fiber direction (3,). Default: [1, 0, 0].	None
dispersion	float	Cone half-angle in radians. 0 = all aligned.	0.0

Returns:

Type	Description
ndarray	Fiber directions (N, 3), unit vectors.

Source code in hyper_surrogate/data/deformation.py

def fiber_directions(
    self,
    n: int,
    preferred: np.ndarray | None = None,
    dispersion: float = 0.0,
) -> np.ndarray:
    """Generate fiber direction vectors.

    Args:
        n: Number of samples.
        preferred: Preferred fiber direction (3,). Default: [1, 0, 0].
        dispersion: Cone half-angle in radians. 0 = all aligned.

    Returns:
        Fiber directions (N, 3), unit vectors.
    """
    if preferred is None:
        preferred = np.array([1.0, 0.0, 0.0])
    preferred = preferred / np.linalg.norm(preferred)

    if dispersion <= 0.0:
        return np.tile(preferred, (n, 1))

    # Sample directions in a cone around preferred
    # Build local frame: preferred = e3, find orthogonal e1, e2
    if abs(preferred[2]) < 0.9:
        t = np.cross(preferred, np.array([0.0, 0.0, 1.0]))
    else:
        t = np.cross(preferred, np.array([1.0, 0.0, 0.0]))
    e1 = t / np.linalg.norm(t)
    e2 = np.cross(preferred, e1)

    phi = self._rng.uniform(0, 2 * np.pi, size=n)
    theta = self._rng.uniform(0, dispersion, size=n)

    x = np.sin(theta) * np.cos(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(theta)

    dirs = x[:, None] * e1 + y[:, None] * e2 + z[:, None] * preferred
    norms = np.linalg.norm(dirs, axis=1, keepdims=True)
    result: np.ndarray = dirs / norms
    return result

MaterialDataset

Bases: Dataset

Wraps (input, target) pairs for training.

Source code in hyper_surrogate/data/dataset.py

class MaterialDataset(Dataset):
    """Wraps (input, target) pairs for training."""

    def __init__(self, inputs: np.ndarray, targets: Any) -> None:
        self.inputs = inputs
        self.targets = targets

    def __len__(self) -> int:
        return len(self.inputs)

    def __getitem__(self, idx: int) -> tuple:
        x = self.inputs[idx]
        y = tuple(t[idx] for t in self.targets) if isinstance(self.targets, tuple) else self.targets[idx]
        return x, y

Normalizer

Standard (zero-mean, unit-variance) normalization with export support.

Source code in hyper_surrogate/data/dataset.py

class Normalizer:
    """Standard (zero-mean, unit-variance) normalization with export support."""

    def __init__(self) -> None:
        self._mean: np.ndarray | None = None
        self._std: np.ndarray | None = None

    def fit(self, data: np.ndarray) -> Normalizer:
        self._mean = data.mean(axis=0)
        self._std = data.std(axis=0)
        self._std[self._std < 1e-12] = 1.0  # type: ignore[index,operator]  # avoid division by zero
        return self

    def transform(self, data: np.ndarray) -> np.ndarray:
        return (data - self._mean) / self._std  # type: ignore[no-any-return]

    def inverse_transform(self, data: np.ndarray) -> np.ndarray:
        return data * self._std + self._mean  # type: ignore[no-any-return]

    @property
    def params(self) -> dict[str, np.ndarray]:
        return {"mean": self._mean, "std": self._std}  # type: ignore[dict-item]

create_datasets(material, n_samples, input_type='invariants', target_type='pk2_voigt', val_fraction=0.15, seed=42)

Generate data, normalize, split, and wrap in datasets.

Source code in hyper_surrogate/data/dataset.py

def create_datasets(
    material: Any,
    n_samples: int,
    input_type: Literal["invariants", "cauchy_green"] = "invariants",
    target_type: Literal["energy", "pk2_voigt", "pk2_voigt+cmat_voigt"] = "pk2_voigt",
    val_fraction: float = 0.15,
    seed: int = 42,
) -> tuple[MaterialDataset, MaterialDataset, Normalizer, Normalizer]:
    """Generate data, normalize, split, and wrap in datasets."""
    from hyper_surrogate.data.deformation import DeformationGenerator
    from hyper_surrogate.mechanics.kinematics import Kinematics

    # Generate deformations
    gen = DeformationGenerator(seed=seed)
    F = gen.combined(n_samples)
    C = Kinematics.right_cauchy_green(F)

    # Compute inputs
    if input_type == "invariants":
        i1 = Kinematics.isochoric_invariant1(C)
        i2 = Kinematics.isochoric_invariant2(C)
        j = np.sqrt(Kinematics.det_invariant(C))  # J = sqrt(det(C))
        if hasattr(material, "is_anisotropic") and material.is_anisotropic:
            num_fibers = getattr(material, "num_fiber_families", 1)
            if num_fibers > 1:
                fiber_invs = Kinematics.fiber_invariants_multi(C, material.fiber_directions)
                inputs = np.column_stack([i1, i2, j, fiber_invs])
            else:
                i4 = Kinematics.fiber_invariant4(C, material.fiber_direction)
                i5 = Kinematics.fiber_invariant5(C, material.fiber_direction)
                inputs = np.column_stack([i1, i2, j, i4, i5])
        else:
            inputs = np.column_stack([i1, i2, j])
    else:  # cauchy_green
        # 6 unique Voigt components: C11, C22, C33, C12, C13, C23
        inputs = np.column_stack([
            C[:, 0, 0],
            C[:, 1, 1],
            C[:, 2, 2],
            C[:, 0, 1],
            C[:, 0, 2],
            C[:, 1, 2],
        ])

    if target_type == "energy":
        energy = material.evaluate_energy(C)  # (N,)
        targets_raw = energy.reshape(-1, 1)

        # Stress target must match input dimensionality for EnergyStressLoss.
        # Invariant inputs (3D+) → dW/d(invariants); cauchy_green (6D) → PK2 Voigt.
        if input_type == "invariants":
            stress_target = material.evaluate_energy_grad_invariants(C)
        else:
            stress_target = _pk2_voigt(material, C)

        in_norm = Normalizer().fit(inputs)
        inputs_normed = in_norm.transform(inputs)

        # EnergyStressLoss computes dW/d(x_norm) via autograd. Since
        # x_norm = (I - mean) / std, the chain rule gives
        # dW/d(x_norm) = dW/dI * std. Scale the stress targets to match.
        stress_target_scaled = stress_target * in_norm.params["std"]

        # Split
        n_val = int(n_samples * val_fraction)
        rng = np.random.default_rng(seed)
        idx = rng.permutation(n_samples)
        train_idx, val_idx = idx[n_val:], idx[:n_val]

        # Energy normalizer is still returned for Fortran export (denormalization),
        # but targets stored raw. Trainer/loss handles raw values.
        energy_norm = Normalizer().fit(targets_raw)

        train_ds = MaterialDataset(
            inputs_normed[train_idx].astype(np.float32),
            (targets_raw[train_idx].astype(np.float32), stress_target_scaled[train_idx].astype(np.float32)),
        )
        val_ds = MaterialDataset(
            inputs_normed[val_idx].astype(np.float32),
            (targets_raw[val_idx].astype(np.float32), stress_target_scaled[val_idx].astype(np.float32)),
        )
        return train_ds, val_ds, in_norm, energy_norm

    elif target_type == "pk2_voigt":
        targets_raw = _pk2_voigt(material, C)
    elif target_type == "pk2_voigt+cmat_voigt":
        cmat_batch = material.evaluate_cmat(C)  # (N, 3, 3, 3, 3)
        # Extract 21 unique Voigt components (upper triangle of 6x6)
        ii1 = [0, 1, 2, 0, 0, 1]
        ii2 = [0, 1, 2, 1, 2, 2]
        cmat_voigt = np.zeros((n_samples, 21))
        k = 0
        for i in range(6):
            for j in range(i, 6):
                cmat_voigt[:, k] = 0.5 * (
                    cmat_batch[:, ii1[i], ii2[i], ii1[j], ii2[j]] + cmat_batch[:, ii1[i], ii2[i], ii2[j], ii1[j]]
                )
                k += 1
        targets_raw = np.column_stack([_pk2_voigt(material, C), cmat_voigt])
    else:
        msg = f"Unknown target_type: {target_type}"
        raise ValueError(msg)

    # Normalize
    in_norm = Normalizer().fit(inputs)
    out_norm = Normalizer().fit(targets_raw)
    inputs_normed = in_norm.transform(inputs)
    targets_normed = out_norm.transform(targets_raw)

    # Split
    n_val = int(n_samples * val_fraction)
    rng = np.random.default_rng(seed)
    idx = rng.permutation(n_samples)
    train_idx, val_idx = idx[n_val:], idx[:n_val]

    train_ds = MaterialDataset(
        inputs_normed[train_idx].astype(np.float32), targets_normed[train_idx].astype(np.float32)
    )
    val_ds = MaterialDataset(inputs_normed[val_idx].astype(np.float32), targets_normed[val_idx].astype(np.float32))

    return train_ds, val_ds, in_norm, out_norm

Models

BranchInfo dataclass

Describes one branch of a multi-branch model (e.g. PolyconvexICNN).

Source code in hyper_surrogate/models/base.py

@dataclass
class BranchInfo:
    """Describes one branch of a multi-branch model (e.g. PolyconvexICNN)."""

    name: str
    input_indices: list[int]
    layers: list[LayerInfo] = field(default_factory=list)

SurrogateModel

Bases: Module, ABC

Source code in hyper_surrogate/models/base.py

class SurrogateModel(nn.Module, ABC):
    @property
    @abstractmethod
    def input_dim(self) -> int: ...

    @property
    @abstractmethod
    def output_dim(self) -> int: ...

    @abstractmethod
    def layer_sequence(self) -> list[LayerInfo]: ...

    def branch_sequence(self) -> list[BranchInfo] | None:
        """Return branch info for multi-branch models. None for single-branch."""
        return None

    def export_weights(self) -> dict[str, np.ndarray]:
        return {k: v.detach().cpu().numpy() for k, v in self.state_dict().items()}

branch_sequence()

Return branch info for multi-branch models. None for single-branch.

Source code in hyper_surrogate/models/base.py

def branch_sequence(self) -> list[BranchInfo] | None:
    """Return branch info for multi-branch models. None for single-branch."""
    return None

ICNN

Bases: SurrogateModel

Input-Convex Neural Network (Amos+ 2017).

Guarantees convexity of output w.r.t. input via: - Non-negative weights on z-path (enforced via softplus) - Skip connections from input to every layer

Source code in hyper_surrogate/models/icnn.py

class ICNN(SurrogateModel):
    """Input-Convex Neural Network (Amos+ 2017).

    Guarantees convexity of output w.r.t. input via:
    - Non-negative weights on z-path (enforced via softplus)
    - Skip connections from input to every layer
    """

    def __init__(
        self,
        input_dim: int,
        hidden_dims: list[int] | None = None,
        activation: str = "softplus",
    ) -> None:
        super().__init__()
        if hidden_dims is None:
            hidden_dims = [64, 64]
        self._input_dim = input_dim
        self._output_dim = 1
        self._activation_name = activation
        self._hidden_dims = hidden_dims

        # First layer: only x-path
        self.wx_layers = nn.ModuleList()
        self.wz_layers = nn.ModuleList()

        # wx_0: input -> hidden_0
        self.wx_layers.append(nn.Linear(input_dim, hidden_dims[0]))

        # Hidden layers: wz (non-negative) + wx (skip)
        for i in range(1, len(hidden_dims)):
            self.wz_layers.append(nn.Linear(hidden_dims[i - 1], hidden_dims[i], bias=False))
            self.wx_layers.append(nn.Linear(input_dim, hidden_dims[i]))

        # Output layer
        self.wz_final = nn.Linear(hidden_dims[-1], 1, bias=False)
        self.wx_final = nn.Linear(input_dim, 1)

        # Activation
        act_map = {"softplus": nn.Softplus(), "relu": nn.ReLU(), "tanh": nn.Tanh()}
        self._activation = act_map[activation]

    @property
    def input_dim(self) -> int:
        return self._input_dim

    @property
    def output_dim(self) -> int:
        return self._output_dim

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # First hidden layer (x-path only)
        z = self._activation(self.wx_layers[0](x))

        # Subsequent hidden layers (z-path with non-neg weights + x skip)
        for wz, wx in zip(self.wz_layers, self.wx_layers[1:], strict=False):
            z = self._activation(F.linear(z, F.softplus(wz.weight)) + wx(x))  # type: ignore[arg-type]

        # Output layer
        return F.linear(z, F.softplus(self.wz_final.weight)) + self.wx_final(x)  # type: ignore[no-any-return]

    def layer_sequence(self) -> list[LayerInfo]:
        result = []
        # First wx layer
        result.append(
            LayerInfo(
                weights="wx_layers.0.weight",
                bias="wx_layers.0.bias",
                activation=self._activation_name,
            )
        )
        # Hidden wz + wx pairs
        for i in range(len(self.wz_layers)):
            result.append(
                LayerInfo(
                    weights=f"wz_layers.{i}.weight",
                    bias=f"wx_layers.{i + 1}.bias",
                    activation=self._activation_name,
                )
            )
        # Output
        result.append(
            LayerInfo(
                weights="wz_final.weight",
                bias="wx_final.bias",
                activation="identity",
            )
        )
        return result

Training

SparseLoss

Bases: Module

Energy-stress loss with L1 regularization on model weights for CANN model discovery.

Encourages sparsity in CANN basis function weights, so surviving terms reveal the minimal constitutive law.

Source code in hyper_surrogate/training/losses.py

class SparseLoss(nn.Module):
    """Energy-stress loss with L1 regularization on model weights for CANN model discovery.

    Encourages sparsity in CANN basis function weights, so surviving terms
    reveal the minimal constitutive law.
    """

    def __init__(self, base_loss: nn.Module, l1_lambda: float = 0.01, weight_param: str = "raw_weights") -> None:
        super().__init__()
        self.base_loss = base_loss
        self.l1_lambda = l1_lambda
        self.weight_param = weight_param

    def forward(self, *args: torch.Tensor, model: nn.Module | None = None, **kwargs: torch.Tensor) -> torch.Tensor:
        loss: torch.Tensor = self.base_loss(*args, **kwargs)
        if model is not None:
            for name, param in model.named_parameters():
                if self.weight_param in name:
                    loss = loss + self.l1_lambda * torch.abs(param).sum()
        return loss

StressTangentLoss

Bases: Module

Source code in hyper_surrogate/training/losses.py

class StressTangentLoss(nn.Module):
    def __init__(self, alpha: float = 1.0, beta: float = 0.1) -> None:
        super().__init__()
        self.alpha = alpha
        self.beta = beta

    def forward(self, pred: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
        """pred and target are both (N, 27): first 6 = stress, rest = tangent."""
        s_pred, c_pred = pred[:, :6], pred[:, 6:]
        s_true, c_true = target[:, :6], target[:, 6:]
        return self.alpha * F.mse_loss(s_pred, s_true) + self.beta * F.mse_loss(c_pred, c_true)

forward(pred, target)

pred and target are both (N, 27): first 6 = stress, rest = tangent.

Source code in hyper_surrogate/training/losses.py

def forward(self, pred: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
    """pred and target are both (N, 27): first 6 = stress, rest = tangent."""
    s_pred, c_pred = pred[:, :6], pred[:, 6:]
    s_true, c_true = target[:, :6], target[:, 6:]
    return self.alpha * F.mse_loss(s_pred, s_true) + self.beta * F.mse_loss(c_pred, c_true)

Export

FortranEmitter

Emits Fortran 90 code for neural network inference.

Source code in hyper_surrogate/export/fortran/emitter.py

class FortranEmitter:
    """Emits Fortran 90 code for neural network inference."""

    ACTIVATIONS: ClassVar[dict[str, str]] = {
        "relu": "max(0.0d0, {x})",
        "tanh": "tanh({x})",
        "sigmoid": "1.0d0 / (1.0d0 + exp(-({x})))",
        "softplus": "log(1.0d0 + exp({x}))",
        "identity": "{x}",
    }

    def __init__(self, exported: ExportedModel) -> None:
        self.exported = exported

    def _format_array_1d(self, arr: np.ndarray, name: str) -> str:
        n = arr.shape[0]
        values = ", ".join(f"{v:.15e}d0" for v in arr.flat)
        return f"  DOUBLE PRECISION, PARAMETER :: {name}({n}) = (/ {values} /)"

    def _format_array_2d(self, arr: np.ndarray, name: str) -> str:
        rows, cols = arr.shape
        values = ", ".join(f"{v:.15e}d0" for v in arr.T.flat)  # Fortran is column-major
        return f"  DOUBLE PRECISION, PARAMETER :: {name}({rows},{cols}) = RESHAPE((/ {values} /), (/ {rows}, {cols} /))"

    def _emit_activation(self, var: str, activation: str, size: int) -> list[str]:
        if activation == "identity":
            return []
        lines = []
        template = self.ACTIVATIONS[activation]
        lines.append(f"  DO i = 1, {size}")
        lines.append(f"    {var}(i) = {template.format(x=f"{var}(i)")}")
        lines.append("  END DO")
        return lines

    def emit_mlp(self) -> str:
        layers = self.exported.layers
        weights = self.exported.weights
        meta = self.exported.metadata
        in_dim = meta["input_dim"]
        out_dim = meta["output_dim"]

        lines = ["MODULE nn_surrogate", "  IMPLICIT NONE", ""]

        # Declare weight/bias arrays as parameters
        for i, layer in enumerate(layers):
            w = weights[layer.weights]
            b = weights[layer.bias]
            lines.append(self._format_array_2d(w, f"w{i}"))
            lines.append(self._format_array_1d(b, f"b{i}"))

        # Normalization parameters
        if self.exported.input_normalizer:
            lines.append(self._format_array_1d(self.exported.input_normalizer["mean"], "in_mean"))
            lines.append(self._format_array_1d(self.exported.input_normalizer["std"], "in_std"))
        if self.exported.output_normalizer:
            lines.append(self._format_array_1d(self.exported.output_normalizer["mean"], "out_mean"))
            lines.append(self._format_array_1d(self.exported.output_normalizer["std"], "out_std"))

        lines.extend(["", "CONTAINS", ""])

        # Subroutine
        lines.append("  SUBROUTINE nn_forward(input, output)")
        lines.append(f"    DOUBLE PRECISION, INTENT(IN) :: input({in_dim})")
        lines.append(f"    DOUBLE PRECISION, INTENT(OUT) :: output({out_dim})")

        # Local variables
        hidden_dims = [weights[layer.weights].shape[0] for layer in layers]
        for i in range(len(layers)):
            lines.append(f"    DOUBLE PRECISION :: z{i}({hidden_dims[i]})")
        lines.append(f"    DOUBLE PRECISION :: x_norm({in_dim})")
        lines.append("    INTEGER :: i")
        lines.append("")

        # Normalize input
        if self.exported.input_normalizer:
            lines.append("    ! Normalize input")
            lines.append("    x_norm = (input - in_mean) / in_std")
        else:
            lines.append("    x_norm = input")
        lines.append("")

        # Forward pass
        for i, layer in enumerate(layers):
            input_var = "x_norm" if i == 0 else f"z{i - 1}"
            lines.append(f"    ! Layer {i}")
            lines.append(f"    z{i} = MATMUL(w{i}, {input_var}) + b{i}")
            lines.extend(["  " + line for line in self._emit_activation(f"z{i}", layer.activation, hidden_dims[i])])
            lines.append("")

        # Copy output and denormalize
        last = f"z{len(layers) - 1}"
        if self.exported.output_normalizer:
            lines.append("    ! Denormalize output")
            lines.append(f"    output = {last} * out_std + out_mean")
        else:
            lines.append(f"    output = {last}")

        lines.extend(["", "  END SUBROUTINE nn_forward", "", "END MODULE nn_surrogate"])

        return "\n".join(lines)

    def emit_icnn(self) -> str:
        layers = self.exported.layers
        weights = self.exported.weights
        meta = self.exported.metadata
        in_dim = meta["input_dim"]

        lines = ["MODULE nn_surrogate", "  IMPLICIT NONE", ""]

        # Declare all weight arrays - ICNN wz weights need softplus pre-applied
        for i, layer in enumerate(layers):
            w_key = layer.weights
            w = weights[w_key]
            if "wz" in w_key:
                w = np.log(1.0 + np.exp(w))  # Pre-apply softplus
            lines.append(self._format_array_2d(w, f"w{i}"))
            b_key = layer.bias
            b = weights[b_key]
            lines.append(self._format_array_1d(b, f"b{i}"))

        # wx skip-connection weights
        wx_keys = sorted([k for k in weights if k.startswith("wx_layers") and "weight" in k])
        for i, wk in enumerate(wx_keys):
            if f"w{i}" not in "\n".join(lines):
                lines.append(self._format_array_2d(weights[wk], f"wx{i}"))

        if self.exported.input_normalizer:
            lines.append(self._format_array_1d(self.exported.input_normalizer["mean"], "in_mean"))
            lines.append(self._format_array_1d(self.exported.input_normalizer["std"], "in_std"))

        lines.extend(["", "CONTAINS", ""])

        # Forward + backward pass for energy and stress
        lines.append("  SUBROUTINE nn_forward(input, energy, stress)")
        lines.append(f"    DOUBLE PRECISION, INTENT(IN) :: input({in_dim})")
        lines.append("    DOUBLE PRECISION, INTENT(OUT) :: energy")
        lines.append(f"    DOUBLE PRECISION, INTENT(OUT) :: stress({in_dim})")
        lines.append(f"    DOUBLE PRECISION :: x_norm({in_dim})")

        hidden_dims = []
        for layer in layers[:-1]:
            w = weights[layer.weights]
            hidden_dims.append(w.shape[0])

        for i, hd in enumerate(hidden_dims):
            lines.append(f"    DOUBLE PRECISION :: z{i}({hd})")
            lines.append(f"    DOUBLE PRECISION :: dz{i}({hd})")

        lines.append(f"    DOUBLE PRECISION :: denergy({in_dim})")
        lines.append("    INTEGER :: i, j")
        lines.append("")

        if self.exported.input_normalizer:
            lines.append("    x_norm = (input - in_mean) / in_std")
        else:
            lines.append("    x_norm = input")
        lines.append("")

        # Forward pass
        lines.append("    ! Forward pass")
        lines.append("    ! Layer 0 (x-path only)")
        lines.append("    z0 = MATMUL(w0, x_norm) + b0")
        lines.extend(["  " + line for line in self._emit_activation("z0", layers[0].activation, hidden_dims[0])])
        lines.append("")

        for i in range(1, len(hidden_dims)):
            lines.append(f"    ! Layer {i} (wz + wx skip)")
            lines.append(f"    z{i} = MATMUL(w{i}, z{i - 1}) + MATMUL(wx{i}, x_norm) + b{i}")
            lines.extend(["  " + line for line in self._emit_activation(f"z{i}", layers[i].activation, hidden_dims[i])])
            lines.append("")

        # Output
        last_hidden = len(hidden_dims) - 1
        last_layer_idx = len(layers) - 1
        lines.append("    ! Output layer")
        lines.append(f"    energy = DOT_PRODUCT(w{last_layer_idx}(1,:), z{last_hidden}) + b{last_layer_idx}(1)")
        lines.append("")

        # Backward pass placeholder
        lines.append("    ! Backward pass (chain rule for stress = d_energy/d_input)")
        lines.append("    energy = 0.0d0")
        lines.append("    stress = 0.0d0")

        lines.extend(["", "  END SUBROUTINE nn_forward", "", "END MODULE nn_surrogate"])

        return "\n".join(lines)

    def emit_polyconvex(self) -> str:
        """Emit Fortran for PolyconvexICNN: per-branch forward + backward."""
        meta = self.exported.metadata
        in_dim = meta["input_dim"]
        branches = meta.get("branches", [])

        lines = ["MODULE nn_surrogate", "  IMPLICIT NONE", ""]

        # Declare weight/bias arrays per branch
        for bi, branch in enumerate(branches):
            self._emit_poly_branch_weights(lines, bi, branch)

        if self.exported.input_normalizer:
            lines.append(self._format_array_1d(self.exported.input_normalizer["mean"], "in_mean"))
            lines.append(self._format_array_1d(self.exported.input_normalizer["std"], "in_std"))

        lines.extend(["", "CONTAINS", ""])

        # Subroutine
        lines.append("  SUBROUTINE nn_forward(input, energy, stress)")
        lines.append(f"    DOUBLE PRECISION, INTENT(IN) :: input({in_dim})")
        lines.append("    DOUBLE PRECISION, INTENT(OUT) :: energy")
        lines.append(f"    DOUBLE PRECISION, INTENT(OUT) :: stress({in_dim})")
        lines.append(f"    DOUBLE PRECISION :: x_norm({in_dim})")
        lines.append("    DOUBLE PRECISION :: branch_energy")
        lines.append("    INTEGER :: i, j")

        # Declare branch-local arrays
        for bi, branch in enumerate(branches):
            self._emit_poly_branch_locals(lines, bi, branch)
        lines.append("")

        # Normalize
        if self.exported.input_normalizer:
            lines.append("    x_norm = (input - in_mean) / in_std")
        else:
            lines.append("    x_norm = input")
        lines.append("")

        lines.append("    energy = 0.0d0")
        lines.append("    stress = 0.0d0")
        lines.append("")

        # Per-branch forward + backward
        for bi, branch in enumerate(branches):
            self._emit_poly_branch_body(lines, bi, branch)

        lines.extend(["  END SUBROUTINE nn_forward", "", "END MODULE nn_surrogate"])
        return "\n".join(lines)

    def _emit_poly_branch_weights(self, lines: list[str], bi: int, branch: dict[str, Any]) -> None:
        """Emit weight/bias parameter arrays for one polyconvex branch."""
        weights = self.exported.weights
        branch_layers = branch["layers"]
        for li, layer in enumerate(branch_layers):
            w = weights[layer["weights"]]
            b = weights[layer["bias"]]
            if "wz" in layer["weights"]:
                w = np.log(1.0 + np.exp(w))
            lines.append(self._format_array_2d(w, f"w_b{bi}_{li}"))
            lines.append(self._format_array_1d(b, f"b_b{bi}_{li}"))
        prefix = f"branches.{bi}."
        wx_keys = sorted([
            k
            for k in weights
            if k.startswith(prefix + "wx_layers") and "weight" in k and k != branch_layers[0]["weights"]
        ])
        for wi, wk in enumerate(wx_keys):
            lines.append(self._format_array_2d(weights[wk], f"wx_b{bi}_{wi + 1}"))

    def _emit_poly_branch_locals(self, lines: list[str], bi: int, branch: dict[str, Any]) -> None:
        """Emit local variable declarations for one polyconvex branch."""
        weights = self.exported.weights
        branch_layers = branch["layers"]
        b_in_dim = len(branch["input_indices"])
        lines.append(f"    DOUBLE PRECISION :: x_b{bi}({b_in_dim})")
        lines.append(f"    DOUBLE PRECISION :: grad_b{bi}({b_in_dim})")
        for li, layer in enumerate(branch_layers[:-1]):
            hd = weights[layer["weights"]].shape[0]
            lines.append(f"    DOUBLE PRECISION :: z_b{bi}_{li}({hd})")
            lines.append(f"    DOUBLE PRECISION :: dz_b{bi}_{li}({hd})")

    def _emit_poly_branch_body(self, lines: list[str], bi: int, branch: dict[str, Any]) -> None:
        """Emit forward + backward pass for one polyconvex branch."""
        weights = self.exported.weights
        branch_layers = branch["layers"]
        indices = branch["input_indices"]
        b_in_dim = len(indices)
        n_hidden = len(branch_layers) - 1

        lines.append(f"    ! === Branch {bi}: inputs [{", ".join(str(i + 1) for i in indices)}] ===")
        for si, idx in enumerate(indices):
            lines.append(f"    x_b{bi}({si + 1}) = x_norm({idx + 1})")
        lines.append("")

        # Forward pass
        lines.append(f"    ! Forward pass branch {bi}")
        lines.append(f"    z_b{bi}_0 = MATMUL(w_b{bi}_0, x_b{bi}) + b_b{bi}_0")
        lines.extend(
            "  " + line
            for line in self._emit_activation(
                f"z_b{bi}_0", branch_layers[0]["activation"], weights[branch_layers[0]["weights"]].shape[0]
            )
        )
        lines.append("")

        for li in range(1, n_hidden):
            layer = branch_layers[li]
            hd = weights[layer["weights"]].shape[0]
            lines.append(
                f"    z_b{bi}_{li} = MATMUL(w_b{bi}_{li}, z_b{bi}_{li - 1}) + MATMUL(wx_b{bi}_{li}, x_b{bi}) + b_b{bi}_{li}"
            )
            lines.extend("  " + line for line in self._emit_activation(f"z_b{bi}_{li}", layer["activation"], hd))
            lines.append("")

        last_li = len(branch_layers) - 1
        lines.append(
            f"    branch_energy = DOT_PRODUCT(w_b{bi}_{last_li}(1,:), z_b{bi}_{n_hidden - 1}) + b_b{bi}_{last_li}(1)"
        )
        lines.append("    energy = energy + branch_energy")
        lines.append("")

        # Backward pass
        lines.append(f"    ! Backward pass branch {bi}")
        lines.append(f"    dz_b{bi}_{n_hidden - 1} = w_b{bi}_{last_li}(1,:)")
        lines.append("")

        for li in range(n_hidden - 1, 0, -1):
            layer = branch_layers[li]
            hd = weights[layer["weights"]].shape[0]
            self._emit_dact_multiply(lines, f"dz_b{bi}_{li}", f"z_b{bi}_{li}", layer["activation"], hd, bi, li)
            prev_hd = weights[branch_layers[li - 1]["weights"]].shape[0]
            lines.append(f"    DO i = 1, {prev_hd}")
            lines.append(f"      dz_b{bi}_{li - 1}(i) = 0.0d0")
            lines.append(f"      DO j = 1, {hd}")
            lines.append(f"        dz_b{bi}_{li - 1}(i) = dz_b{bi}_{li - 1}(i) + w_b{bi}_{li}(j, i) * dz_b{bi}_{li}(j)")
            lines.append("      END DO")
            lines.append("    END DO")
            lines.append("")

        layer0 = branch_layers[0]
        hd0 = weights[layer0["weights"]].shape[0]
        self._emit_dact_multiply(lines, f"dz_b{bi}_0", f"z_b{bi}_0", layer0["activation"], hd0, bi, 0)

        lines.append(f"    DO i = 1, {b_in_dim}")
        lines.append(f"      grad_b{bi}(i) = 0.0d0")
        lines.append(f"      DO j = 1, {hd0}")
        lines.append(f"        grad_b{bi}(i) = grad_b{bi}(i) + w_b{bi}_0(j, i) * dz_b{bi}_0(j)")
        lines.append("      END DO")
        lines.append("    END DO")
        lines.append("")

        for si, idx in enumerate(indices):
            lines.append(f"    stress({idx + 1}) = stress({idx + 1}) + grad_b{bi}({si + 1})")
        lines.append("")

    @staticmethod
    def _emit_dact_multiply(
        lines: list[str], dz_var: str, z_var: str, activation: str, size: int, bi: int, li: int
    ) -> None:
        """Multiply dz by activation derivative in-place."""
        if activation == "identity":
            return
        if activation == "softplus":
            lines.append(f"    DO i = 1, {size}")
            lines.append(f"      {dz_var}(i) = {dz_var}(i) / (1.0d0 + exp(-{z_var}(i)))")
            lines.append("    END DO")
        elif activation == "tanh":
            lines.append(f"    DO i = 1, {size}")
            lines.append(f"      {dz_var}(i) = {dz_var}(i) * (1.0d0 - {z_var}(i)**2)")
            lines.append("    END DO")
        elif activation == "relu":
            lines.append(f"    DO i = 1, {size}")
            lines.append(f"      IF ({z_var}(i) <= 0.0d0) {dz_var}(i) = 0.0d0")
            lines.append("    END DO")

    def emit(self) -> str:
        arch = self.exported.metadata.get("architecture", "mlp")
        if arch == "mlp":
            return self.emit_mlp()
        elif arch == "icnn":
            return self.emit_icnn()
        elif arch == "polyconvexicnn":
            return self.emit_polyconvex()
        else:
            msg = f"Unknown architecture: {arch}"
            raise ValueError(msg)

    def write(self, path: str) -> None:
        with open(path, "w") as f:
            f.write(self.emit())

emit_polyconvex()

Emit Fortran for PolyconvexICNN: per-branch forward + backward.

Source code in hyper_surrogate/export/fortran/emitter.py

def emit_polyconvex(self) -> str:
    """Emit Fortran for PolyconvexICNN: per-branch forward + backward."""
    meta = self.exported.metadata
    in_dim = meta["input_dim"]
    branches = meta.get("branches", [])

    lines = ["MODULE nn_surrogate", "  IMPLICIT NONE", ""]

    # Declare weight/bias arrays per branch
    for bi, branch in enumerate(branches):
        self._emit_poly_branch_weights(lines, bi, branch)

    if self.exported.input_normalizer:
        lines.append(self._format_array_1d(self.exported.input_normalizer["mean"], "in_mean"))
        lines.append(self._format_array_1d(self.exported.input_normalizer["std"], "in_std"))

    lines.extend(["", "CONTAINS", ""])

    # Subroutine
    lines.append("  SUBROUTINE nn_forward(input, energy, stress)")
    lines.append(f"    DOUBLE PRECISION, INTENT(IN) :: input({in_dim})")
    lines.append("    DOUBLE PRECISION, INTENT(OUT) :: energy")
    lines.append(f"    DOUBLE PRECISION, INTENT(OUT) :: stress({in_dim})")
    lines.append(f"    DOUBLE PRECISION :: x_norm({in_dim})")
    lines.append("    DOUBLE PRECISION :: branch_energy")
    lines.append("    INTEGER :: i, j")

    # Declare branch-local arrays
    for bi, branch in enumerate(branches):
        self._emit_poly_branch_locals(lines, bi, branch)
    lines.append("")

    # Normalize
    if self.exported.input_normalizer:
        lines.append("    x_norm = (input - in_mean) / in_std")
    else:
        lines.append("    x_norm = input")
    lines.append("")

    lines.append("    energy = 0.0d0")
    lines.append("    stress = 0.0d0")
    lines.append("")

    # Per-branch forward + backward
    for bi, branch in enumerate(branches):
        self._emit_poly_branch_body(lines, bi, branch)

    lines.extend(["  END SUBROUTINE nn_forward", "", "END MODULE nn_surrogate"])
    return "\n".join(lines)

Hybrid UMAT generator: NN-based SEF with analytical continuum mechanics.

The NN learns W(invariants). Everything else — kinematics, stress, and tangent — is computed analytically in Fortran:

DFGRD1 → C → invariants → NN(W) → backprop(dW/dI, d²W/dI²)
       → PK2 → Cauchy stress → spatial tangent + Jaumann correction

Supports

• Isotropic: W(I1_bar, I2_bar, J) — input_dim=3
• Anisotropic: W(I1_bar, I2_bar, J, I4, I5) — input_dim=5

HybridUMATEmitter

Emit a complete Abaqus UMAT subroutine with NN-based strain energy.

Source code in hyper_surrogate/export/fortran/hybrid.py

class HybridUMATEmitter:
    """Emit a complete Abaqus UMAT subroutine with NN-based strain energy."""

    SUPPORTED_ARCHITECTURES = ("mlp", "polyconvexicnn")

    def __init__(self, exported: ExportedModel) -> None:
        arch = exported.metadata.get("architecture")
        if arch not in self.SUPPORTED_ARCHITECTURES:
            msg = f"HybridUMATEmitter supports {self.SUPPORTED_ARCHITECTURES}, got '{arch}'"
            raise ValueError(msg)
        if exported.metadata.get("output_dim") != 1:
            msg = "HybridUMATEmitter requires a scalar (output_dim=1) energy model"
            raise ValueError(msg)
        self.exported = exported

    # ------------------------------------------------------------------
    # Fortran formatting helpers
    # ------------------------------------------------------------------

    @staticmethod
    def _fmt_1d(arr: np.ndarray, name: str) -> str:
        n = arr.shape[0]
        vals = ", &\n    ".join(", ".join(f"{v:.15e}d0" for v in arr[i : i + 4]) for i in range(0, n, 4))
        return f"DOUBLE PRECISION, PARAMETER :: {name}({n}) = (/ &\n    {vals} /)"

    @staticmethod
    def _fmt_2d(arr: np.ndarray, name: str) -> str:
        rows, cols = arr.shape
        # Fortran is column-major
        vals = ", &\n    ".join(
            ", ".join(f"{v:.15e}d0" for v in arr.T.flat[i : i + 4]) for i in range(0, rows * cols, 4)
        )
        return (
            f"DOUBLE PRECISION, PARAMETER :: {name}({rows},{cols}) = RESHAPE((/ &\n    {vals} /), (/ {rows}, {cols} /))"
        )

    # ------------------------------------------------------------------
    # Code generation
    # ------------------------------------------------------------------

    def _emit_nn_parameters(self) -> str:
        """Emit weight/bias arrays and normalizer constants."""
        arch = self.exported.metadata.get("architecture")
        if arch == "polyconvexicnn":
            return self._emit_poly_nn_parameters()
        return self._emit_mlp_nn_parameters()

    def _emit_mlp_nn_parameters(self) -> str:
        """Emit MLP weight/bias arrays."""
        lines: list[str] = []
        layers = self.exported.layers
        weights = self.exported.weights

        for i, layer in enumerate(layers):
            w = weights[layer.weights]
            b = weights[layer.bias]
            lines.append(self._fmt_2d(w, f"w{i}"))
            lines.append(self._fmt_1d(b, f"b{i}"))

        if self.exported.input_normalizer:
            lines.append(self._fmt_1d(self.exported.input_normalizer["mean"], "in_mean"))
            lines.append(self._fmt_1d(self.exported.input_normalizer["std"], "in_std"))

        return "\n".join(lines)

    def _emit_poly_nn_parameters(self) -> str:
        """Emit per-branch ICNN weight/bias arrays."""
        lines: list[str] = []
        weights = self.exported.weights
        branches = self.exported.metadata["branches"]

        for bi, branch in enumerate(branches):
            branch_layers = branch["layers"]
            for li, layer in enumerate(branch_layers):
                w = weights[layer["weights"]]
                b = weights[layer["bias"]]
                # Pre-apply softplus to wz weights (non-negative constraint)
                if "wz" in layer["weights"]:
                    w = np.log(1.0 + np.exp(w))
                lines.append(self._fmt_2d(w, f"w_b{bi}_{li}"))
                lines.append(self._fmt_1d(b, f"b_b{bi}_{li}"))

            # wx skip-connection weights (layers 1..L-1 have wx_layers)
            prefix = f"branches.{bi}."
            wx_keys = sorted([
                k
                for k in weights
                if k.startswith(prefix + "wx_layers") and "weight" in k and k != branch_layers[0]["weights"]
            ])
            for wi, wk in enumerate(wx_keys):
                lines.append(self._fmt_2d(weights[wk], f"wx_b{bi}_{wi + 1}"))
            # wx_final skip
            wx_final_key = prefix + "wx_final.weight"
            if wx_final_key in weights:
                lines.append(self._fmt_2d(weights[wx_final_key], f"wxf_b{bi}"))

        if self.exported.input_normalizer:
            lines.append(self._fmt_1d(self.exported.input_normalizer["mean"], "in_mean"))
            lines.append(self._fmt_1d(self.exported.input_normalizer["std"], "in_std"))

        return "\n".join(lines)

    def _emit_nn_forward_and_backward(self) -> str:  # noqa: C901
        """Emit NN forward pass + analytical backward pass for dW/dI and d²W/dI²."""
        arch = self.exported.metadata.get("architecture")
        if arch == "polyconvexicnn":
            return self._emit_poly_nn_forward_and_backward()
        layers = self.exported.layers
        weights = self.exported.weights
        n_layers = len(layers)
        hidden_dims = [weights[layer.weights].shape[0] for layer in layers]
        in_dim = self.exported.metadata["input_dim"]

        lines: list[str] = []

        # Local variables for forward pass
        for i in range(n_layers):
            lines.append(f"DOUBLE PRECISION :: z{i}({hidden_dims[i]})")
            lines.append(f"DOUBLE PRECISION :: a{i}({hidden_dims[i]})")  # pre-activation
            lines.append(f"DOUBLE PRECISION :: dact{i}({hidden_dims[i]})")  # activation derivative
            lines.append(f"DOUBLE PRECISION :: d2act{i}({hidden_dims[i]})")  # second derivative
        lines.append(f"DOUBLE PRECISION :: x_norm({in_dim})")
        lines.append("")

        # Local variables for backward pass
        for i in range(n_layers - 1, -1, -1):
            lines.append(f"DOUBLE PRECISION :: delta{i}({hidden_dims[i]})")
        lines.append(f"DOUBLE PRECISION :: grad_x({in_dim})")
        lines.append("")

        # For Jacobian propagation and Hessian
        for i in range(n_layers):
            lines.append(f"DOUBLE PRECISION :: P{i}({hidden_dims[i]},{in_dim})")
            lines.append(f"DOUBLE PRECISION :: J{i}({hidden_dims[i]},{in_dim})")
        lines.append(f"DOUBLE PRECISION :: d2W_dx2({in_dim},{in_dim})")
        lines.append("DOUBLE PRECISION :: coeff")
        lines.append("")

        # --- Normalize input ---
        lines.append("! Normalize invariants")
        lines.append("x_norm = (nn_input - in_mean) / in_std")
        lines.append("")

        # --- Forward pass with d2act ---
        lines.append("! Forward pass")
        for i, layer in enumerate(layers):
            input_var = "x_norm" if i == 0 else f"z{i - 1}"
            lines.append(f"a{i} = MATMUL(w{i}, {input_var}) + b{i}")
            act = layer.activation
            if act == "identity":
                lines.append(f"z{i} = a{i}")
                lines.append(f"dact{i} = 1.0d0")
                lines.append(f"d2act{i} = 0.0d0")
            elif act == "tanh":
                lines.append(f"z{i} = tanh(a{i})")
                lines.append(f"dact{i} = 1.0d0 - z{i}**2")
                lines.append(f"d2act{i} = -2.0d0 * z{i} * dact{i}")
            elif act == "softplus":
                lines.append(f"z{i} = log(1.0d0 + exp(a{i}))")
                lines.append(f"dact{i} = 1.0d0 / (1.0d0 + exp(-a{i}))")
                lines.append(f"d2act{i} = dact{i} * (1.0d0 - dact{i})")
            elif act == "relu":
                lines.append(f"DO ii = 1, {hidden_dims[i]}")
                lines.append(f"  z{i}(ii) = max(0.0d0, a{i}(ii))")
                lines.append(f"  IF (a{i}(ii) > 0.0d0) THEN")
                lines.append(f"    dact{i}(ii) = 1.0d0")
                lines.append("  ELSE")
                lines.append(f"    dact{i}(ii) = 0.0d0")
                lines.append("  END IF")
                lines.append(f"  d2act{i}(ii) = 0.0d0")
                lines.append("END DO")
            elif act == "sigmoid":
                lines.append(f"z{i} = 1.0d0 / (1.0d0 + exp(-a{i}))")
                lines.append(f"dact{i} = z{i} * (1.0d0 - z{i})")
                lines.append(f"d2act{i} = dact{i} * (1.0d0 - 2.0d0 * z{i})")
            lines.append("")

        # W = z_{last}(1) (scalar output)
        last = n_layers - 1
        lines.append(f"W_nn = z{last}(1)")
        lines.append("")

        # --- Backward pass: dW/dx_norm ---
        lines.append("! Backward pass: dW/d(x_norm)")
        # delta for output layer
        lines.append(f"delta{last}(1) = dact{last}(1)")
        # Backpropagate
        for i in range(last - 1, -1, -1):
            lines.append(f"DO ii = 1, {hidden_dims[i]}")
            lines.append(f"  delta{i}(ii) = 0.0d0")
            lines.append(f"  DO jj = 1, {hidden_dims[i + 1]}")
            lines.append(f"    delta{i}(ii) = delta{i}(ii) + w{i + 1}(jj, ii) * delta{i + 1}(jj)")
            lines.append("  END DO")
            lines.append(f"  delta{i}(ii) = delta{i}(ii) * dact{i}(ii)")
            lines.append("END DO")

        # grad_x = W_0^T * delta_0
        lines.append(f"DO ii = 1, {in_dim}")
        lines.append("  grad_x(ii) = 0.0d0")
        lines.append(f"  DO jj = 1, {hidden_dims[0]}")
        lines.append("    grad_x(ii) = grad_x(ii) + w0(jj, ii) * delta0(jj)")
        lines.append("  END DO")
        lines.append("END DO")
        lines.append("")

        # dW/dI = dW/dx_norm / std (chain rule for normalization)
        lines.append("! Convert gradient to raw invariant space: dW/dI = dW/dx_norm / std")
        for k in range(in_dim):
            lines.append(f"dW_dI({k + 1}) = grad_x({k + 1}) / in_std({k + 1})")
        lines.append("")

        # --- Jacobian propagation: P_i = W_i * J_{i-1}, J_i = diag(dact_i) * P_i ---
        lines.append("! Jacobian propagation (forward mode)")
        # P0 = W0, J0 = diag(dact0) * W0
        lines.append(f"DO ii = 1, {hidden_dims[0]}")
        lines.append(f"  DO jj = 1, {in_dim}")
        lines.append("    P0(ii, jj) = w0(ii, jj)")
        lines.append("    J0(ii, jj) = dact0(ii) * w0(ii, jj)")
        lines.append("  END DO")
        lines.append("END DO")

        for i in range(1, n_layers):
            # P_i = W_i * J_{i-1}  (pre-activation Jacobian)
            # J_i = diag(dact_i) * P_i
            lines.append(f"DO ii = 1, {hidden_dims[i]}")
            lines.append(f"  DO jj = 1, {in_dim}")
            lines.append(f"    P{i}(ii, jj) = 0.0d0")
            lines.append(f"    DO kk = 1, {hidden_dims[i - 1]}")
            lines.append(f"      P{i}(ii, jj) = P{i}(ii, jj) + w{i}(ii, kk) * J{i - 1}(kk, jj)")
            lines.append("    END DO")
            lines.append(f"    J{i}(ii, jj) = dact{i}(ii) * P{i}(ii, jj)")
            lines.append("  END DO")
            lines.append("END DO")
        lines.append("")

        # --- Exact analytical Hessian: d²W/dx² = Σ_i P_i^T diag(beta_i * d2act_i) P_i ---
        # beta_i = delta_i / dact_i = dW/dz_i (sensitivity to post-activation)
        lines.append("! Exact analytical Hessian: d²W/dx_norm²")
        lines.append("d2W_dx2 = 0.0d0")
        for i in range(n_layers):
            act = layers[i].activation
            # Skip layers with zero d2act (relu, identity) — they contribute nothing
            if act in ("relu", "identity"):
                lines.append(f"! Layer {i} ({act}): d2act=0, no Hessian contribution")
                continue
            lines.append(f"! Layer {i} ({act})")
            lines.append(f"DO ii = 1, {hidden_dims[i]}")
            # coeff = beta_i(ii) * d2act_i(ii) = (delta_i(ii) / dact_i(ii)) * d2act_i(ii)
            lines.append(f"  coeff = delta{i}(ii) / dact{i}(ii) * d2act{i}(ii)")
            lines.append(f"  DO jj = 1, {in_dim}")
            lines.append(f"    DO kk = 1, {in_dim}")
            lines.append(f"      d2W_dx2(jj, kk) = d2W_dx2(jj, kk) + coeff * P{i}(ii, jj) * P{i}(ii, kk)")
            lines.append("    END DO")
            lines.append("  END DO")
            lines.append("END DO")
        lines.append("")

        # Convert to raw invariant space: d2W/dI2(k,l) = d2W/dx2(k,l) / (std(k) * std(l))
        lines.append("! Convert Hessian to raw invariant space: d2W/dI2(k,l) = d2W/dx2(k,l) / (std(k)*std(l))")
        lines.append(f"DO ii = 1, {in_dim}")
        lines.append(f"  DO jj = 1, {in_dim}")
        lines.append("    d2W_dI2(ii, jj) = d2W_dx2(ii, jj) / (in_std(ii) * in_std(jj))")
        lines.append("  END DO")
        lines.append("END DO")

        return "\n".join(lines)

    # ------------------------------------------------------------------
    # Anisotropic Fortran snippets (generalized for N fiber families)
    # ------------------------------------------------------------------

    @staticmethod
    def _emit_aniso_declarations(num_fibers: int) -> str:
        """Extra local variable declarations for anisotropic UMAT."""
        lines = ["  ! Fiber invariants (anisotropic)"]
        for k in range(num_fibers):
            lines.append(f"  DOUBLE PRECISION :: a0_{k + 1}(3), Ca0_{k + 1}(3), I4_{k + 1}, I5_{k + 1}")
            lines.append(f"  DOUBLE PRECISION :: dI4_{k + 1}_dC(3,3), dI5_{k + 1}_dC(3,3)")
        return "\n".join(lines) + "\n"

    @staticmethod
    def _emit_aniso_invariants(num_fibers: int) -> str:
        """Compute fiber invariants I4, I5 from C and a0 for each fiber family."""
        lines: list[str] = []
        for k in range(num_fibers):
            prop_start = k * 3 + 1
            lines.append(f"""
  ! Fiber family {k + 1}: direction from props({prop_start}:{prop_start + 2})
  a0_{k + 1}(1) = props({prop_start})
  a0_{k + 1}(2) = props({prop_start + 1})
  a0_{k + 1}(3) = props({prop_start + 2})

  ! Ca0_{k + 1} = C * a0_{k + 1}
  DO ii = 1, 3
    Ca0_{k + 1}(ii) = 0.0d0
    DO jj = 1, 3
      Ca0_{k + 1}(ii) = Ca0_{k + 1}(ii) + C(ii,jj) * a0_{k + 1}(jj)
    END DO
  END DO

  ! I4_{k + 1} = a0_{k + 1} . C . a0_{k + 1}
  I4_{k + 1} = 0.0d0
  DO ii = 1, 3
    I4_{k + 1} = I4_{k + 1} + a0_{k + 1}(ii) * Ca0_{k + 1}(ii)
  END DO

  ! I5_{k + 1} = a0_{k + 1} . C^2 . a0_{k + 1} = Ca0_{k + 1} . Ca0_{k + 1}
  I5_{k + 1} = 0.0d0
  DO ii = 1, 3
    I5_{k + 1} = I5_{k + 1} + Ca0_{k + 1}(ii) * Ca0_{k + 1}(ii)
  END DO""")
        return "\n".join(lines)

    @staticmethod
    def _emit_aniso_nn_input(num_fibers: int) -> str:
        """Set fiber invariant entries in nn_input."""
        lines: list[str] = []
        for k in range(num_fibers):
            idx_i4 = 3 + 2 * k + 1  # Fortran 1-indexed
            idx_i5 = 3 + 2 * k + 2
            lines.append(f"  nn_input({idx_i4}) = I4_{k + 1}")
            lines.append(f"  nn_input({idx_i5}) = I5_{k + 1}")
        return "\n".join(lines) + "\n"

    @staticmethod
    def _emit_aniso_didC(num_fibers: int) -> str:
        """Compute dI4/dC and dI5/dC for each fiber family."""
        lines: list[str] = []
        for k in range(num_fibers):
            lines.append(f"""
  ! dI4_{k + 1}/dC = a0_{k + 1} (x) a0_{k + 1}
  DO ii = 1, 3
    DO jj = 1, 3
      dI4_{k + 1}_dC(ii,jj) = a0_{k + 1}(ii) * a0_{k + 1}(jj)
    END DO
  END DO

  ! dI5_{k + 1}/dC = a0_{k + 1} (x) Ca0_{k + 1} + Ca0_{k + 1} (x) a0_{k + 1}
  DO ii = 1, 3
    DO jj = 1, 3
      dI5_{k + 1}_dC(ii,jj) = a0_{k + 1}(ii) * Ca0_{k + 1}(jj) + Ca0_{k + 1}(ii) * a0_{k + 1}(jj)
    END DO
  END DO""")
        return "\n".join(lines)

    @staticmethod
    def _emit_aniso_pk2_terms(num_fibers: int) -> str:
        """Additional PK2 terms for fiber invariants."""
        lines: list[str] = []
        for k in range(num_fibers):
            idx_i4 = 3 + 2 * k + 1  # Fortran 1-indexed
            idx_i5 = 3 + 2 * k + 2
            lines.append(f" &\n        + dW_dI({idx_i4}) * dI4_{k + 1}_dC(ii,jj)")
            lines.append(f" &\n        + dW_dI({idx_i5}) * dI5_{k + 1}_dC(ii,jj)")
        return "".join(lines)

    @staticmethod
    def _emit_aniso_dIdC_pack(num_fibers: int) -> str:
        """Pack fiber dI/dC into dIdC array."""
        lines = ["    DO ii = 1, 3", "      DO jj = 1, 3"]
        for k in range(num_fibers):
            idx_i4 = 3 + 2 * k + 1
            idx_i5 = 3 + 2 * k + 2
            lines.append(f"        dIdC(ii, jj, {idx_i4}) = dI4_{k + 1}_dC(ii, jj)")
            lines.append(f"        dIdC(ii, jj, {idx_i5}) = dI5_{k + 1}_dC(ii, jj)")
        lines.extend(["      END DO", "    END DO"])
        return "\n".join(lines) + "\n"

    @staticmethod
    def _emit_aniso_d2IdC2(num_fibers: int) -> str:
        """Second derivatives d²I4/dCdC=0, d²I5/dCdC for tangent."""
        lines: list[str] = []
        for k in range(num_fibers):
            idx_i5 = 3 + 2 * k + 2  # Fortran 1-indexed
            lines.append(f"""
    ! d²I4_{k + 1}/dCdC = 0 (dI4/dC = a0 (x) a0 is constant w.r.t. C)

    ! d²I5_{k + 1}/(dC_AB dC_CD)
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = 0.5d0 * ( &
                a0_{k + 1}(ii)*a0_{k + 1}(ll)*eye3(jj,kk) + a0_{k + 1}(ii)*a0_{k + 1}(kk)*eye3(jj,ll) &
              + a0_{k + 1}(jj)*a0_{k + 1}(ll)*eye3(ii,kk) + a0_{k + 1}(jj)*a0_{k + 1}(kk)*eye3(ii,ll))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI({idx_i5}) * val
          END DO
        END DO
      END DO
    END DO""")
        return "\n".join(lines)

    def _emit_poly_nn_forward_and_backward(self) -> str:  # noqa: C901
        """Emit polyconvex ICNN forward/backward with per-branch Hessian."""
        weights = self.exported.weights
        branches = self.exported.metadata["branches"]
        in_dim = self.exported.metadata["input_dim"]

        lines: list[str] = []

        # Local variables
        lines.append(f"DOUBLE PRECISION :: x_norm({in_dim})")
        for bi, branch in enumerate(branches):
            b_layers = branch["layers"]
            b_in = len(branch["input_indices"])
            n_hidden = len(b_layers) - 1
            lines.append(f"DOUBLE PRECISION :: xb{bi}({b_in})")
            for li in range(n_hidden):
                hd = weights[b_layers[li]["weights"]].shape[0]
                lines.append(f"DOUBLE PRECISION :: z_b{bi}_{li}({hd})")
                lines.append(f"DOUBLE PRECISION :: a_b{bi}_{li}({hd})")
                lines.append(f"DOUBLE PRECISION :: dact_b{bi}_{li}({hd})")
                lines.append(f"DOUBLE PRECISION :: d2act_b{bi}_{li}({hd})")
                lines.append(f"DOUBLE PRECISION :: delta_b{bi}_{li}({hd})")
                lines.append(f"DOUBLE PRECISION :: P_b{bi}_{li}({hd},{b_in})")
                lines.append(f"DOUBLE PRECISION :: J_b{bi}_{li}({hd},{b_in})")
            lines.append(f"DOUBLE PRECISION :: grad_b{bi}({b_in})")
            lines.append(f"DOUBLE PRECISION :: d2W_b{bi}({b_in},{b_in})")
        lines.append("DOUBLE PRECISION :: coeff, branch_W")
        lines.append("")

        # Normalize input
        lines.append("! Normalize invariants")
        lines.append("x_norm = (nn_input - in_mean) / in_std")
        lines.append("")

        # Initialize outputs
        lines.append("W_nn = 0.0d0")
        lines.append("dW_dI = 0.0d0")
        lines.append("d2W_dI2 = 0.0d0")
        lines.append("")

        # Per-branch forward + backward + Hessian
        for bi, branch in enumerate(branches):
            b_layers = branch["layers"]
            indices = branch["input_indices"]
            b_in = len(indices)
            n_hidden = len(b_layers) - 1

            lines.append(f"! ---- Branch {bi}: inputs [{", ".join(str(i + 1) for i in indices)}] ----")

            # Slice input
            for si, idx in enumerate(indices):
                lines.append(f"xb{bi}({si + 1}) = x_norm({idx + 1})")
            lines.append("")

            # Forward pass
            lines.append(f"! Forward pass branch {bi}")
            # Layer 0: x-path only
            lines.append(f"a_b{bi}_0 = MATMUL(w_b{bi}_0, xb{bi}) + b_b{bi}_0")
            act0 = b_layers[0]["activation"]
            hd0 = weights[b_layers[0]["weights"]].shape[0]
            self._emit_icnn_act(lines, bi, 0, act0, hd0)
            lines.append("")

            # Hidden layers with wz + wx skip
            for li in range(1, n_hidden):
                layer = b_layers[li]
                hd = weights[layer["weights"]].shape[0]
                act = layer["activation"]
                lines.append(
                    f"a_b{bi}_{li} = MATMUL(w_b{bi}_{li}, z_b{bi}_{li - 1}) + MATMUL(wx_b{bi}_{li}, xb{bi}) + b_b{bi}_{li}"
                )
                self._emit_icnn_act(lines, bi, li, act, hd)
                lines.append("")

            # Output: identity activation
            last_li = len(b_layers) - 1
            last_hidden = n_hidden - 1
            lines.append(
                f"branch_W = DOT_PRODUCT(w_b{bi}_{last_li}(1,:), z_b{bi}_{last_hidden}) + DOT_PRODUCT(wxf_b{bi}(1,:), xb{bi}) + b_b{bi}_{last_li}(1)"
            )
            lines.append("W_nn = W_nn + branch_W")
            lines.append("")

            # Backward pass
            lines.append(f"! Backward pass branch {bi}")
            # delta for last hidden layer
            hd_last = weights[b_layers[last_hidden]["weights"]].shape[0]
            lines.append(f"DO ii = 1, {hd_last}")
            lines.append(f"  delta_b{bi}_{last_hidden}(ii) = w_b{bi}_{last_li}(1, ii) * dact_b{bi}_{last_hidden}(ii)")
            lines.append("END DO")

            # Backprop through hidden layers
            for li in range(last_hidden - 1, -1, -1):
                hd = weights[b_layers[li]["weights"]].shape[0]
                hd_next = weights[b_layers[li + 1]["weights"]].shape[0]
                lines.append(f"DO ii = 1, {hd}")
                lines.append(f"  delta_b{bi}_{li}(ii) = 0.0d0")
                lines.append(f"  DO jj = 1, {hd_next}")
                lines.append(
                    f"    delta_b{bi}_{li}(ii) = delta_b{bi}_{li}(ii) + w_b{bi}_{li + 1}(jj, ii) * delta_b{bi}_{li + 1}(jj)"
                )
                lines.append("  END DO")
                lines.append(f"  delta_b{bi}_{li}(ii) = delta_b{bi}_{li}(ii) * dact_b{bi}_{li}(ii)")
                lines.append("END DO")

            # grad_b = W0^T * delta_0 + wx_final^T
            lines.append(f"DO ii = 1, {b_in}")
            lines.append(f"  grad_b{bi}(ii) = wxf_b{bi}(1, ii)")
            lines.append(f"  DO jj = 1, {hd0}")
            lines.append(f"    grad_b{bi}(ii) = grad_b{bi}(ii) + w_b{bi}_0(jj, ii) * delta_b{bi}_0(jj)")
            lines.append("  END DO")
            lines.append("END DO")
            lines.append("")

            # Scatter gradient
            for si, idx in enumerate(indices):
                lines.append(f"dW_dI({idx + 1}) = dW_dI({idx + 1}) + grad_b{bi}({si + 1}) / in_std({idx + 1})")
            lines.append("")

            # Jacobian propagation for Hessian
            lines.append(f"! Jacobian propagation branch {bi}")
            # P0 = W0, J0 = diag(dact0) * W0
            lines.append(f"DO ii = 1, {hd0}")
            lines.append(f"  DO jj = 1, {b_in}")
            lines.append(f"    P_b{bi}_0(ii, jj) = w_b{bi}_0(ii, jj)")
            lines.append(f"    J_b{bi}_0(ii, jj) = dact_b{bi}_0(ii) * w_b{bi}_0(ii, jj)")
            lines.append("  END DO")
            lines.append("END DO")

            for li in range(1, n_hidden):
                hd = weights[b_layers[li]["weights"]].shape[0]
                hd_prev = weights[b_layers[li - 1]["weights"]].shape[0]
                # P_i = Wz_i * J_{i-1} + Wx_i  (ICNN skip connection)
                lines.append(f"DO ii = 1, {hd}")
                lines.append(f"  DO jj = 1, {b_in}")
                lines.append(f"    P_b{bi}_{li}(ii, jj) = wx_b{bi}_{li}(ii, jj)")
                lines.append(f"    DO kk = 1, {hd_prev}")
                lines.append(
                    f"      P_b{bi}_{li}(ii, jj) = P_b{bi}_{li}(ii, jj) + w_b{bi}_{li}(ii, kk) * J_b{bi}_{li - 1}(kk, jj)"
                )
                lines.append("    END DO")
                lines.append(f"    J_b{bi}_{li}(ii, jj) = dact_b{bi}_{li}(ii) * P_b{bi}_{li}(ii, jj)")
                lines.append("  END DO")
                lines.append("END DO")
            lines.append("")

            # Hessian: d2W_b = sum_i P_i^T diag(beta_i * d2act_i) P_i
            lines.append(f"! Hessian branch {bi}")
            lines.append(f"d2W_b{bi} = 0.0d0")
            for li in range(n_hidden):
                act = b_layers[li]["activation"]
                if act in ("relu", "identity"):
                    lines.append(f"! Layer {li} ({act}): d2act=0, skip")
                    continue
                hd = weights[b_layers[li]["weights"]].shape[0]
                lines.append(f"DO ii = 1, {hd}")
                lines.append(f"  coeff = delta_b{bi}_{li}(ii) / dact_b{bi}_{li}(ii) * d2act_b{bi}_{li}(ii)")
                lines.append(f"  DO jj = 1, {b_in}")
                lines.append(f"    DO kk = 1, {b_in}")
                lines.append(
                    f"      d2W_b{bi}(jj, kk) = d2W_b{bi}(jj, kk) + coeff * P_b{bi}_{li}(ii, jj) * P_b{bi}_{li}(ii, kk)"
                )
                lines.append("    END DO")
                lines.append("  END DO")
                lines.append("END DO")
            lines.append("")

            # Scatter Hessian (block-diagonal)
            for si, idx_i in enumerate(indices):
                for sj, idx_j in enumerate(indices):
                    lines.append(
                        f"d2W_dI2({idx_i + 1},{idx_j + 1}) = d2W_dI2({idx_i + 1},{idx_j + 1})"
                        f" + d2W_b{bi}({si + 1},{sj + 1}) / (in_std({idx_i + 1}) * in_std({idx_j + 1}))"
                    )
            lines.append("")

        return "\n".join(lines)

    @staticmethod
    def _emit_icnn_act(lines: list[str], bi: int, li: int, act: str, hd: int) -> None:
        """Emit activation + dact + d2act for ICNN branch layer."""
        if act == "identity":
            lines.append(f"z_b{bi}_{li} = a_b{bi}_{li}")
            lines.append(f"dact_b{bi}_{li} = 1.0d0")
            lines.append(f"d2act_b{bi}_{li} = 0.0d0")
        elif act == "softplus":
            lines.append(f"z_b{bi}_{li} = log(1.0d0 + exp(a_b{bi}_{li}))")
            lines.append(f"dact_b{bi}_{li} = 1.0d0 / (1.0d0 + exp(-a_b{bi}_{li}))")
            lines.append(f"d2act_b{bi}_{li} = dact_b{bi}_{li} * (1.0d0 - dact_b{bi}_{li})")
        elif act == "tanh":
            lines.append(f"z_b{bi}_{li} = tanh(a_b{bi}_{li})")
            lines.append(f"dact_b{bi}_{li} = 1.0d0 - z_b{bi}_{li}**2")
            lines.append(f"d2act_b{bi}_{li} = -2.0d0 * z_b{bi}_{li} * dact_b{bi}_{li}")
        elif act == "relu":
            lines.append(f"DO ii = 1, {hd}")
            lines.append(f"  z_b{bi}_{li}(ii) = max(0.0d0, a_b{bi}_{li}(ii))")
            lines.append(f"  IF (a_b{bi}_{li}(ii) > 0.0d0) THEN")
            lines.append(f"    dact_b{bi}_{li}(ii) = 1.0d0")
            lines.append("  ELSE")
            lines.append(f"    dact_b{bi}_{li}(ii) = 0.0d0")
            lines.append("  END IF")
            lines.append(f"  d2act_b{bi}_{li}(ii) = 0.0d0")
            lines.append("END DO")

    def emit(self) -> str:
        """Generate the complete hybrid UMAT Fortran code."""
        today = datetime.datetime.now().strftime("%Y-%m-%d")
        nn_params = self._emit_nn_parameters()
        layers = self.exported.layers
        weights = self.exported.weights
        hidden_dims = [weights[layer.weights].shape[0] for layer in layers]
        in_dim = self.exported.metadata["input_dim"]
        num_fibers = (in_dim - 3) // 2
        aniso = num_fibers > 0
        n_inv = in_dim

        if aniso:
            inv_parts = ["I1_bar", "I2_bar", "J"]
            for k in range(num_fibers):
                inv_parts.extend([f"I4_{k + 1}", f"I5_{k + 1}"])
            input_desc = ", ".join(inv_parts)
        else:
            input_desc = "I1_bar, I2_bar, J"

        code = f"""\
!>********************************************************************
!> Hybrid UMAT: NN-based strain energy with analytical mechanics
!>
!> Generated: {today}
!> Architecture: MLP {" x ".join(str(d) for d in hidden_dims)}
!> Input: {input_desc}  ->  Output: W (strain energy)
!>
!> The neural network provides W({input_desc}).
!> Stress and tangent are derived analytically via chain rule.
!>   Cauchy = (1/J) * F * PK2 * F^T
!>   Tangent = push-forward of material tangent + Jaumann correction
{"!> Fiber directions: " + ", ".join(f"a0_{k + 1} from props({k * 3 + 1}:{k * 3 + 3})" for k in range(num_fibers)) if aniso else ""}
!>********************************************************************

MODULE nn_sef
  IMPLICIT NONE

  ! NN weights and biases
  {nn_params}

CONTAINS

  SUBROUTINE nn_eval(nn_input, W_nn, dW_dI, d2W_dI2)
    !> Evaluate NN: W({input_desc}), dW/dI, and d²W/dI² via backprop
    DOUBLE PRECISION, INTENT(IN)  :: nn_input({n_inv})
    DOUBLE PRECISION, INTENT(OUT) :: W_nn
    DOUBLE PRECISION, INTENT(OUT) :: dW_dI({n_inv})
    DOUBLE PRECISION, INTENT(OUT) :: d2W_dI2({n_inv},{n_inv})

    ! Local variables
    INTEGER :: ii, jj, kk
    {self._emit_nn_forward_and_backward()}

  END SUBROUTINE nn_eval

END MODULE nn_sef


SUBROUTINE umat(stress, statev, ddsdde, sse, spd, scd, rpl, &
    ddsddt, drplde, drpldt, stran, dstran, time, dtime, temp, &
    dtemp, predef, dpred, cmname, ndi, nshr, ntens, nstatev, &
    props, nprops, coords, drot, pnewdt, celent, dfgrd0, &
    dfgrd1, noel, npt, layer, kspt, kstep, kinc)

  USE nn_sef
  IMPLICIT NONE

  ! Standard UMAT interface
  INTEGER, INTENT(IN OUT) :: noel, npt, layer, kspt, kstep, kinc
  INTEGER, INTENT(IN OUT) :: ndi, nshr, ntens, nstatev, nprops
  DOUBLE PRECISION, INTENT(IN OUT) :: sse, spd, scd, rpl, dtime
  DOUBLE PRECISION, INTENT(IN OUT) :: drpldt, temp, dtemp, pnewdt, celent
  CHARACTER(LEN=8), INTENT(IN OUT) :: cmname
  DOUBLE PRECISION, INTENT(IN OUT) :: stress(ntens), statev(nstatev)
  DOUBLE PRECISION, INTENT(IN OUT) :: ddsdde(ntens, ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: ddsddt(ntens), drplde(ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: stran(ntens), dstran(ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: time(2), predef(1), dpred(1)
  DOUBLE PRECISION, INTENT(IN)     :: props(nprops)
  DOUBLE PRECISION, INTENT(IN OUT) :: coords(3), drot(3, 3)
  DOUBLE PRECISION, INTENT(IN OUT) :: dfgrd0(3, 3), dfgrd1(3, 3)

  ! Local variables
  DOUBLE PRECISION :: C(3,3), invC(3,3), detC, detF
  DOUBLE PRECISION :: I1, I2, I1_bar, I2_bar, Jac
  DOUBLE PRECISION :: trC, trC2
  DOUBLE PRECISION :: nn_input({n_inv}), W_nn, dW_dI({n_inv})
  DOUBLE PRECISION :: PK2(3,3), sigma_full(3,3)
  DOUBLE PRECISION :: dI1_dC(3,3), dI2_dC(3,3), dJ_dC(3,3)
  DOUBLE PRECISION :: eye3(3,3)
  DOUBLE PRECISION :: Jm23, Jm43
  INTEGER :: ii, jj, kk, ll
{self._emit_aniso_declarations(num_fibers) if aniso else ""}
  ! d²W/dI² from analytical NN Hessian
  DOUBLE PRECISION :: d2W_dI2({n_inv},{n_inv})

  ! For tangent: material tangent in reference config
  DOUBLE PRECISION :: dPK2_dC(3,3,3,3)
  ! Spatial tangent + Jaumann
  DOUBLE PRECISION :: cmat_spatial(3,3,3,3)

  ! Voigt indices
  INTEGER :: v1(6), v2(6)
  DATA v1 /1, 2, 3, 1, 1, 2/
  DATA v2 /1, 2, 3, 2, 3, 3/

  ! Identity tensor
  eye3 = 0.0d0
  eye3(1,1) = 1.0d0
  eye3(2,2) = 1.0d0
  eye3(3,3) = 1.0d0

  ! ================================================================
  ! 1. Kinematics: F -> C -> invariants
  ! ================================================================
  ! Right Cauchy-Green: C = F^T F
  DO ii = 1, 3
    DO jj = 1, 3
      C(ii,jj) = 0.0d0
      DO kk = 1, 3
        C(ii,jj) = C(ii,jj) + dfgrd1(kk,ii) * dfgrd1(kk,jj)
      END DO
    END DO
  END DO

  ! Determinants
  detC = C(1,1)*(C(2,2)*C(3,3) - C(2,3)*C(3,2)) &
       - C(1,2)*(C(2,1)*C(3,3) - C(2,3)*C(3,1)) &
       + C(1,3)*(C(2,1)*C(3,2) - C(2,2)*C(3,1))
  Jac = sqrt(detC)
  detF = Jac

  ! Invariants
  trC = C(1,1) + C(2,2) + C(3,3)
  trC2 = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      trC2 = trC2 + C(ii,jj) * C(jj,ii)
    END DO
  END DO

  Jm23 = detC**(-1.0d0/3.0d0)
  Jm43 = detC**(-2.0d0/3.0d0)

  I1_bar = trC * Jm23
  I2_bar = 0.5d0 * (trC**2 - trC2) * Jm43
{self._emit_aniso_invariants(num_fibers) if aniso else ""}
  ! ================================================================
  ! 2. NN evaluation: W({input_desc}) and dW/dI
  ! ================================================================
  nn_input(1) = I1_bar
  nn_input(2) = I2_bar
  nn_input(3) = Jac
{self._emit_aniso_nn_input(num_fibers) if aniso else ""}  CALL nn_eval(nn_input, W_nn, dW_dI, d2W_dI2)

  ! Store strain energy
  sse = W_nn

  ! ================================================================
  ! 3. Analytical derivatives: dI/dC
  ! ================================================================
  ! Inverse of C (for dJ/dC)
  invC(1,1) = (C(2,2)*C(3,3) - C(2,3)*C(3,2)) / detC
  invC(2,2) = (C(1,1)*C(3,3) - C(1,3)*C(3,1)) / detC
  invC(3,3) = (C(1,1)*C(2,2) - C(1,2)*C(2,1)) / detC
  invC(1,2) = (C(1,3)*C(3,2) - C(1,2)*C(3,3)) / detC
  invC(2,1) = invC(1,2)
  invC(1,3) = (C(1,2)*C(2,3) - C(1,3)*C(2,2)) / detC
  invC(3,1) = invC(1,3)
  invC(2,3) = (C(1,3)*C(2,1) - C(1,1)*C(2,3)) / detC
  invC(3,2) = invC(2,3)

  ! dI1_bar/dC = J^(-2/3) * (I - (1/3)*I1_bar * C^-1)
  DO ii = 1, 3
    DO jj = 1, 3
      dI1_dC(ii,jj) = Jm23 * (eye3(ii,jj) - (1.0d0/3.0d0) * trC * invC(ii,jj))
    END DO
  END DO

  ! dI2_bar/dC = J^(-4/3) * (trC*I - C - (2/3)*I2_bar*J^(4/3) * C^-1)
  DO ii = 1, 3
    DO jj = 1, 3
      dI2_dC(ii,jj) = Jm43 * (trC * eye3(ii,jj) - C(ii,jj)) &
                     - (2.0d0/3.0d0) * I2_bar * invC(ii,jj)
    END DO
  END DO

  ! dJ/dC = (1/2) * J * C^-1
  DO ii = 1, 3
    DO jj = 1, 3
      dJ_dC(ii,jj) = 0.5d0 * Jac * invC(ii,jj)
    END DO
  END DO
{self._emit_aniso_didC(num_fibers) if aniso else ""}
  ! ================================================================
  ! 4. PK2 stress: S = 2 * dW/dC = 2 * sum_k (dW/dIk * dIk/dC)
  ! ================================================================
  DO ii = 1, 3
    DO jj = 1, 3
      PK2(ii,jj) = 2.0d0 * ( &
          dW_dI(1) * dI1_dC(ii,jj) &
        + dW_dI(2) * dI2_dC(ii,jj) &
        + dW_dI(3) * dJ_dC(ii,jj){self._emit_aniso_pk2_terms(num_fibers) if aniso else ""} )
    END DO
  END DO

  ! ================================================================
  ! 5. Cauchy stress: sigma = (1/J) * F * S * F^T
  ! ================================================================
  sigma_full = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          sigma_full(ii,jj) = sigma_full(ii,jj) &
            + dfgrd1(ii,kk) * PK2(kk,ll) * dfgrd1(jj,ll)
        END DO
      END DO
    END DO
  END DO
  sigma_full = sigma_full / detF

  ! To Voigt: stress(1..6) = [s11, s22, s33, s12, s13, s23]
  DO ii = 1, ntens
    stress(ii) = sigma_full(v1(ii), v2(ii))
  END DO

  ! ================================================================
  ! 6. Material tangent: C_ABCD = 4 * d²W/dC²
  !    d²W/dCdC = sum_k sum_l (d²W/dIk*dIl) * (dIk/dC) x (dIl/dC)
  !             + sum_k (dW/dIk) * (d²Ik/dCdC)
  ! ================================================================
  dPK2_dC = 0.0d0

  BLOCK
    DOUBLE PRECISION :: dIdC(3, 3, {n_inv})  ! dIdC(:,:,k) = dIk/dC
    DOUBLE PRECISION :: val

    DO ii = 1, 3
      DO jj = 1, 3
        dIdC(ii, jj, 1) = dI1_dC(ii, jj)
        dIdC(ii, jj, 2) = dI2_dC(ii, jj)
        dIdC(ii, jj, 3) = dJ_dC(ii, jj)
      END DO
    END DO
{self._emit_aniso_dIdC_pack(num_fibers) if aniso else ""}
    ! C_ABCD += 4 * sum_k sum_l d²W/dIk*dIl * (dIk/dC)_AB * (dIl/dC)_CD
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = 0.0d0
            DO mm = 1, {n_inv}
              DO nn = 1, {n_inv}
                val = val + d2W_dI2(mm, nn) * dIdC(ii, jj, mm) * dIdC(kk, ll, nn)
              END DO
            END DO
            dPK2_dC(ii, jj, kk, ll) = 4.0d0 * val
          END DO
        END DO
      END DO
    END DO

    ! Term 2: 4 * sum_k dW/dIk * d²Ik/dCdC

    ! d²J/(dC_AB dC_CD) = J/4 * (invC_AB*invC_CD - invC_AC*invC_BD - invC_AD*invC_BC)
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = 0.25d0 * Jac * ( &
                invC(ii,jj) * invC(kk,ll) &
              - 0.5d0 * (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(3) * val
          END DO
        END DO
      END DO
    END DO

    ! d²I1_bar/dCdC
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = Jm23 * ( &
              (-1.0d0/3.0d0) * (eye3(ii,jj)*invC(kk,ll) + invC(ii,jj)*eye3(kk,ll)) &
              + (1.0d0/9.0d0) * trC * invC(ii,jj)*invC(kk,ll) &
              + (1.0d0/3.0d0) * trC * 0.5d0 * &
                  (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(1) * val
          END DO
        END DO
      END DO
    END DO

    ! d²I2_bar/dCdC
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = Jm43 * ( &
              eye3(ii,jj)*eye3(kk,ll) &
              - 0.5d0*(eye3(ii,kk)*eye3(jj,ll) + eye3(ii,ll)*eye3(jj,kk)) &
              + (-2.0d0/3.0d0) * (trC*eye3(ii,jj) - C(ii,jj)) * invC(kk,ll) &
              + (-2.0d0/3.0d0) * invC(ii,jj) * (trC*eye3(kk,ll) - C(kk,ll)) &
              + (4.0d0/9.0d0) * I2_bar / Jm43 * invC(ii,jj) * invC(kk,ll) &
              + (2.0d0/3.0d0) * I2_bar / Jm43 * 0.5d0 * &
                  (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(2) * val
          END DO
        END DO
      END DO
    END DO
{self._emit_aniso_d2IdC2(num_fibers) if aniso else ""}
  END BLOCK

  ! ================================================================
  ! 7. Push forward to spatial tangent: c_ijkl = (1/J) F_iA F_jB F_kC F_lD C_ABCD
  ! ================================================================
  cmat_spatial = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          val = 0.0d0
          DO mm = 1, 3
            DO nn = 1, 3
              DO pp = 1, 3
                DO qq = 1, 3
                  val = val + dfgrd1(ii,mm)*dfgrd1(jj,nn)*dfgrd1(kk,pp)*dfgrd1(ll,qq) &
                            * dPK2_dC(mm,nn,pp,qq)
                END DO
              END DO
            END DO
          END DO
          cmat_spatial(ii,jj,kk,ll) = val / detF
        END DO
      END DO
    END DO
  END DO

  ! Jaumann correction: c_ijkl += 0.5*(delta_ik*sigma_jl + sigma_ik*delta_jl
  !                                   + delta_il*sigma_jk + sigma_il*delta_jk)
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          cmat_spatial(ii,jj,kk,ll) = cmat_spatial(ii,jj,kk,ll) + 0.5d0 * ( &
            eye3(ii,kk)*sigma_full(jj,ll) + sigma_full(ii,kk)*eye3(jj,ll) &
          + eye3(ii,ll)*sigma_full(jj,kk) + sigma_full(ii,ll)*eye3(jj,kk))
        END DO
      END DO
    END DO
  END DO

  ! To Voigt: ddsdde(6,6)
  DO ii = 1, ntens
    DO jj = 1, ntens
      ddsdde(ii, jj) = cmat_spatial(v1(ii), v2(ii), v1(jj), v2(jj))
    END DO
  END DO

RETURN
END SUBROUTINE umat
"""
        return code

    def write(self, path: str | Path) -> None:
        """Write the hybrid UMAT to a file."""
        Path(path).write_text(self.emit())

emit()

Generate the complete hybrid UMAT Fortran code.

Source code in hyper_surrogate/export/fortran/hybrid.py

    def emit(self) -> str:
        """Generate the complete hybrid UMAT Fortran code."""
        today = datetime.datetime.now().strftime("%Y-%m-%d")
        nn_params = self._emit_nn_parameters()
        layers = self.exported.layers
        weights = self.exported.weights
        hidden_dims = [weights[layer.weights].shape[0] for layer in layers]
        in_dim = self.exported.metadata["input_dim"]
        num_fibers = (in_dim - 3) // 2
        aniso = num_fibers > 0
        n_inv = in_dim

        if aniso:
            inv_parts = ["I1_bar", "I2_bar", "J"]
            for k in range(num_fibers):
                inv_parts.extend([f"I4_{k + 1}", f"I5_{k + 1}"])
            input_desc = ", ".join(inv_parts)
        else:
            input_desc = "I1_bar, I2_bar, J"

        code = f"""\
!>********************************************************************
!> Hybrid UMAT: NN-based strain energy with analytical mechanics
!>
!> Generated: {today}
!> Architecture: MLP {" x ".join(str(d) for d in hidden_dims)}
!> Input: {input_desc}  ->  Output: W (strain energy)
!>
!> The neural network provides W({input_desc}).
!> Stress and tangent are derived analytically via chain rule.
!>   Cauchy = (1/J) * F * PK2 * F^T
!>   Tangent = push-forward of material tangent + Jaumann correction
{"!> Fiber directions: " + ", ".join(f"a0_{k + 1} from props({k * 3 + 1}:{k * 3 + 3})" for k in range(num_fibers)) if aniso else ""}
!>********************************************************************

MODULE nn_sef
  IMPLICIT NONE

  ! NN weights and biases
  {nn_params}

CONTAINS

  SUBROUTINE nn_eval(nn_input, W_nn, dW_dI, d2W_dI2)
    !> Evaluate NN: W({input_desc}), dW/dI, and d²W/dI² via backprop
    DOUBLE PRECISION, INTENT(IN)  :: nn_input({n_inv})
    DOUBLE PRECISION, INTENT(OUT) :: W_nn
    DOUBLE PRECISION, INTENT(OUT) :: dW_dI({n_inv})
    DOUBLE PRECISION, INTENT(OUT) :: d2W_dI2({n_inv},{n_inv})

    ! Local variables
    INTEGER :: ii, jj, kk
    {self._emit_nn_forward_and_backward()}

  END SUBROUTINE nn_eval

END MODULE nn_sef


SUBROUTINE umat(stress, statev, ddsdde, sse, spd, scd, rpl, &
    ddsddt, drplde, drpldt, stran, dstran, time, dtime, temp, &
    dtemp, predef, dpred, cmname, ndi, nshr, ntens, nstatev, &
    props, nprops, coords, drot, pnewdt, celent, dfgrd0, &
    dfgrd1, noel, npt, layer, kspt, kstep, kinc)

  USE nn_sef
  IMPLICIT NONE

  ! Standard UMAT interface
  INTEGER, INTENT(IN OUT) :: noel, npt, layer, kspt, kstep, kinc
  INTEGER, INTENT(IN OUT) :: ndi, nshr, ntens, nstatev, nprops
  DOUBLE PRECISION, INTENT(IN OUT) :: sse, spd, scd, rpl, dtime
  DOUBLE PRECISION, INTENT(IN OUT) :: drpldt, temp, dtemp, pnewdt, celent
  CHARACTER(LEN=8), INTENT(IN OUT) :: cmname
  DOUBLE PRECISION, INTENT(IN OUT) :: stress(ntens), statev(nstatev)
  DOUBLE PRECISION, INTENT(IN OUT) :: ddsdde(ntens, ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: ddsddt(ntens), drplde(ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: stran(ntens), dstran(ntens)
  DOUBLE PRECISION, INTENT(IN OUT) :: time(2), predef(1), dpred(1)
  DOUBLE PRECISION, INTENT(IN)     :: props(nprops)
  DOUBLE PRECISION, INTENT(IN OUT) :: coords(3), drot(3, 3)
  DOUBLE PRECISION, INTENT(IN OUT) :: dfgrd0(3, 3), dfgrd1(3, 3)

  ! Local variables
  DOUBLE PRECISION :: C(3,3), invC(3,3), detC, detF
  DOUBLE PRECISION :: I1, I2, I1_bar, I2_bar, Jac
  DOUBLE PRECISION :: trC, trC2
  DOUBLE PRECISION :: nn_input({n_inv}), W_nn, dW_dI({n_inv})
  DOUBLE PRECISION :: PK2(3,3), sigma_full(3,3)
  DOUBLE PRECISION :: dI1_dC(3,3), dI2_dC(3,3), dJ_dC(3,3)
  DOUBLE PRECISION :: eye3(3,3)
  DOUBLE PRECISION :: Jm23, Jm43
  INTEGER :: ii, jj, kk, ll
{self._emit_aniso_declarations(num_fibers) if aniso else ""}
  ! d²W/dI² from analytical NN Hessian
  DOUBLE PRECISION :: d2W_dI2({n_inv},{n_inv})

  ! For tangent: material tangent in reference config
  DOUBLE PRECISION :: dPK2_dC(3,3,3,3)
  ! Spatial tangent + Jaumann
  DOUBLE PRECISION :: cmat_spatial(3,3,3,3)

  ! Voigt indices
  INTEGER :: v1(6), v2(6)
  DATA v1 /1, 2, 3, 1, 1, 2/
  DATA v2 /1, 2, 3, 2, 3, 3/

  ! Identity tensor
  eye3 = 0.0d0
  eye3(1,1) = 1.0d0
  eye3(2,2) = 1.0d0
  eye3(3,3) = 1.0d0

  ! ================================================================
  ! 1. Kinematics: F -> C -> invariants
  ! ================================================================
  ! Right Cauchy-Green: C = F^T F
  DO ii = 1, 3
    DO jj = 1, 3
      C(ii,jj) = 0.0d0
      DO kk = 1, 3
        C(ii,jj) = C(ii,jj) + dfgrd1(kk,ii) * dfgrd1(kk,jj)
      END DO
    END DO
  END DO

  ! Determinants
  detC = C(1,1)*(C(2,2)*C(3,3) - C(2,3)*C(3,2)) &
       - C(1,2)*(C(2,1)*C(3,3) - C(2,3)*C(3,1)) &
       + C(1,3)*(C(2,1)*C(3,2) - C(2,2)*C(3,1))
  Jac = sqrt(detC)
  detF = Jac

  ! Invariants
  trC = C(1,1) + C(2,2) + C(3,3)
  trC2 = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      trC2 = trC2 + C(ii,jj) * C(jj,ii)
    END DO
  END DO

  Jm23 = detC**(-1.0d0/3.0d0)
  Jm43 = detC**(-2.0d0/3.0d0)

  I1_bar = trC * Jm23
  I2_bar = 0.5d0 * (trC**2 - trC2) * Jm43
{self._emit_aniso_invariants(num_fibers) if aniso else ""}
  ! ================================================================
  ! 2. NN evaluation: W({input_desc}) and dW/dI
  ! ================================================================
  nn_input(1) = I1_bar
  nn_input(2) = I2_bar
  nn_input(3) = Jac
{self._emit_aniso_nn_input(num_fibers) if aniso else ""}  CALL nn_eval(nn_input, W_nn, dW_dI, d2W_dI2)

  ! Store strain energy
  sse = W_nn

  ! ================================================================
  ! 3. Analytical derivatives: dI/dC
  ! ================================================================
  ! Inverse of C (for dJ/dC)
  invC(1,1) = (C(2,2)*C(3,3) - C(2,3)*C(3,2)) / detC
  invC(2,2) = (C(1,1)*C(3,3) - C(1,3)*C(3,1)) / detC
  invC(3,3) = (C(1,1)*C(2,2) - C(1,2)*C(2,1)) / detC
  invC(1,2) = (C(1,3)*C(3,2) - C(1,2)*C(3,3)) / detC
  invC(2,1) = invC(1,2)
  invC(1,3) = (C(1,2)*C(2,3) - C(1,3)*C(2,2)) / detC
  invC(3,1) = invC(1,3)
  invC(2,3) = (C(1,3)*C(2,1) - C(1,1)*C(2,3)) / detC
  invC(3,2) = invC(2,3)

  ! dI1_bar/dC = J^(-2/3) * (I - (1/3)*I1_bar * C^-1)
  DO ii = 1, 3
    DO jj = 1, 3
      dI1_dC(ii,jj) = Jm23 * (eye3(ii,jj) - (1.0d0/3.0d0) * trC * invC(ii,jj))
    END DO
  END DO

  ! dI2_bar/dC = J^(-4/3) * (trC*I - C - (2/3)*I2_bar*J^(4/3) * C^-1)
  DO ii = 1, 3
    DO jj = 1, 3
      dI2_dC(ii,jj) = Jm43 * (trC * eye3(ii,jj) - C(ii,jj)) &
                     - (2.0d0/3.0d0) * I2_bar * invC(ii,jj)
    END DO
  END DO

  ! dJ/dC = (1/2) * J * C^-1
  DO ii = 1, 3
    DO jj = 1, 3
      dJ_dC(ii,jj) = 0.5d0 * Jac * invC(ii,jj)
    END DO
  END DO
{self._emit_aniso_didC(num_fibers) if aniso else ""}
  ! ================================================================
  ! 4. PK2 stress: S = 2 * dW/dC = 2 * sum_k (dW/dIk * dIk/dC)
  ! ================================================================
  DO ii = 1, 3
    DO jj = 1, 3
      PK2(ii,jj) = 2.0d0 * ( &
          dW_dI(1) * dI1_dC(ii,jj) &
        + dW_dI(2) * dI2_dC(ii,jj) &
        + dW_dI(3) * dJ_dC(ii,jj){self._emit_aniso_pk2_terms(num_fibers) if aniso else ""} )
    END DO
  END DO

  ! ================================================================
  ! 5. Cauchy stress: sigma = (1/J) * F * S * F^T
  ! ================================================================
  sigma_full = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          sigma_full(ii,jj) = sigma_full(ii,jj) &
            + dfgrd1(ii,kk) * PK2(kk,ll) * dfgrd1(jj,ll)
        END DO
      END DO
    END DO
  END DO
  sigma_full = sigma_full / detF

  ! To Voigt: stress(1..6) = [s11, s22, s33, s12, s13, s23]
  DO ii = 1, ntens
    stress(ii) = sigma_full(v1(ii), v2(ii))
  END DO

  ! ================================================================
  ! 6. Material tangent: C_ABCD = 4 * d²W/dC²
  !    d²W/dCdC = sum_k sum_l (d²W/dIk*dIl) * (dIk/dC) x (dIl/dC)
  !             + sum_k (dW/dIk) * (d²Ik/dCdC)
  ! ================================================================
  dPK2_dC = 0.0d0

  BLOCK
    DOUBLE PRECISION :: dIdC(3, 3, {n_inv})  ! dIdC(:,:,k) = dIk/dC
    DOUBLE PRECISION :: val

    DO ii = 1, 3
      DO jj = 1, 3
        dIdC(ii, jj, 1) = dI1_dC(ii, jj)
        dIdC(ii, jj, 2) = dI2_dC(ii, jj)
        dIdC(ii, jj, 3) = dJ_dC(ii, jj)
      END DO
    END DO
{self._emit_aniso_dIdC_pack(num_fibers) if aniso else ""}
    ! C_ABCD += 4 * sum_k sum_l d²W/dIk*dIl * (dIk/dC)_AB * (dIl/dC)_CD
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = 0.0d0
            DO mm = 1, {n_inv}
              DO nn = 1, {n_inv}
                val = val + d2W_dI2(mm, nn) * dIdC(ii, jj, mm) * dIdC(kk, ll, nn)
              END DO
            END DO
            dPK2_dC(ii, jj, kk, ll) = 4.0d0 * val
          END DO
        END DO
      END DO
    END DO

    ! Term 2: 4 * sum_k dW/dIk * d²Ik/dCdC

    ! d²J/(dC_AB dC_CD) = J/4 * (invC_AB*invC_CD - invC_AC*invC_BD - invC_AD*invC_BC)
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = 0.25d0 * Jac * ( &
                invC(ii,jj) * invC(kk,ll) &
              - 0.5d0 * (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(3) * val
          END DO
        END DO
      END DO
    END DO

    ! d²I1_bar/dCdC
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = Jm23 * ( &
              (-1.0d0/3.0d0) * (eye3(ii,jj)*invC(kk,ll) + invC(ii,jj)*eye3(kk,ll)) &
              + (1.0d0/9.0d0) * trC * invC(ii,jj)*invC(kk,ll) &
              + (1.0d0/3.0d0) * trC * 0.5d0 * &
                  (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(1) * val
          END DO
        END DO
      END DO
    END DO

    ! d²I2_bar/dCdC
    DO ii = 1, 3
      DO jj = 1, 3
        DO kk = 1, 3
          DO ll = 1, 3
            val = Jm43 * ( &
              eye3(ii,jj)*eye3(kk,ll) &
              - 0.5d0*(eye3(ii,kk)*eye3(jj,ll) + eye3(ii,ll)*eye3(jj,kk)) &
              + (-2.0d0/3.0d0) * (trC*eye3(ii,jj) - C(ii,jj)) * invC(kk,ll) &
              + (-2.0d0/3.0d0) * invC(ii,jj) * (trC*eye3(kk,ll) - C(kk,ll)) &
              + (4.0d0/9.0d0) * I2_bar / Jm43 * invC(ii,jj) * invC(kk,ll) &
              + (2.0d0/3.0d0) * I2_bar / Jm43 * 0.5d0 * &
                  (invC(ii,kk)*invC(jj,ll) + invC(ii,ll)*invC(jj,kk)))
            dPK2_dC(ii,jj,kk,ll) = dPK2_dC(ii,jj,kk,ll) + 4.0d0 * dW_dI(2) * val
          END DO
        END DO
      END DO
    END DO
{self._emit_aniso_d2IdC2(num_fibers) if aniso else ""}
  END BLOCK

  ! ================================================================
  ! 7. Push forward to spatial tangent: c_ijkl = (1/J) F_iA F_jB F_kC F_lD C_ABCD
  ! ================================================================
  cmat_spatial = 0.0d0
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          val = 0.0d0
          DO mm = 1, 3
            DO nn = 1, 3
              DO pp = 1, 3
                DO qq = 1, 3
                  val = val + dfgrd1(ii,mm)*dfgrd1(jj,nn)*dfgrd1(kk,pp)*dfgrd1(ll,qq) &
                            * dPK2_dC(mm,nn,pp,qq)
                END DO
              END DO
            END DO
          END DO
          cmat_spatial(ii,jj,kk,ll) = val / detF
        END DO
      END DO
    END DO
  END DO

  ! Jaumann correction: c_ijkl += 0.5*(delta_ik*sigma_jl + sigma_ik*delta_jl
  !                                   + delta_il*sigma_jk + sigma_il*delta_jk)
  DO ii = 1, 3
    DO jj = 1, 3
      DO kk = 1, 3
        DO ll = 1, 3
          cmat_spatial(ii,jj,kk,ll) = cmat_spatial(ii,jj,kk,ll) + 0.5d0 * ( &
            eye3(ii,kk)*sigma_full(jj,ll) + sigma_full(ii,kk)*eye3(jj,ll) &
          + eye3(ii,ll)*sigma_full(jj,kk) + sigma_full(ii,ll)*eye3(jj,kk))
        END DO
      END DO
    END DO
  END DO

  ! To Voigt: ddsdde(6,6)
  DO ii = 1, ntens
    DO jj = 1, ntens
      ddsdde(ii, jj) = cmat_spatial(v1(ii), v2(ii), v1(jj), v2(jj))
    END DO
  END DO

RETURN
END SUBROUTINE umat
"""
        return code

write(path)

Write the hybrid UMAT to a file.

Source code in hyper_surrogate/export/fortran/hybrid.py

def write(self, path: str | Path) -> None:
    """Write the hybrid UMAT to a file."""
    Path(path).write_text(self.emit())

UMATHandler

Source code in hyper_surrogate/export/fortran/analytical.py

class UMATHandler:
    def __init__(self, material_model: Material) -> None:
        """
        Initialize the UMAT handler with a specific material model.

        Args:
            material_model: The material model (e.g., NeoHooke) to generate UMAT code for.
        """
        self.material = material_model
        self.sigma_code = None
        self.smat_code = None

    @staticmethod
    def common_subexpressions(tensor: sym.Matrix, var_name: str) -> list[str]:
        """
        Perform common subexpression elimination on a vector or matrix and generate Fortran code.

        Args:
            var_name (str): The base name for the variables in the Fortran code.

        Returns:
            tuple: A tuple containing Fortran code for auxiliary variables and reduced expressions.
        """
        # Extract individual components
        # tensor_components = [tensor[i] for i in range(tensor.shape[0])]
        tensor_components = [tensor[i, j] for i in range(tensor.shape[0]) for j in range(tensor.shape[1])]
        # Convert to a matrix to check shape
        tensor_matrix = sym.Matrix(tensor)
        # Perform common subexpression elimination
        replacements, reduced_exprs = sym.cse(tensor_components)

        # Generate Fortran code for auxiliary variables (replacements)
        aux_code = [
            sym.fcode(
                expr,
                standard=90,
                source_format="free",
                assign_to=sym.fcode(var, standard=90, source_format="free"),
            )
            for var, expr in replacements
        ]

        # Generate Fortran code for reduced expressions
        if tensor_matrix.shape[1] == 1:  # If vector
            reduced_code = [
                sym.fcode(
                    expr,
                    standard=90,
                    source_format="free",
                    assign_to=f"{var_name}({i + 1})",
                )
                for i, expr in enumerate(reduced_exprs)
            ]
        else:  # If matrix
            _, cols = tensor.shape
            reduced_code = [
                sym.fcode(
                    expr,
                    standard=90,
                    source_format="free",
                    assign_to=f"{var_name}({i // cols + 1},{i % cols + 1})",
                )
                for i, expr in enumerate(reduced_exprs)
            ]

        return aux_code + reduced_code

    @property
    def f(self) -> sym.Matrix:
        """The deformation gradient tensor."""
        return sym.Matrix(3, 3, lambda i, j: sym.Symbol(f"DFGRD1({i + 1},{j + 1})"))

    @property
    def sub_exp(self) -> dict:
        """Substitution expressions for the right Cauchy-Green tensor."""
        c = self.f.T * self.f
        return {self.material.handler.c_tensor[i, j]: c[i, j] for i in range(3) for j in range(3)}

    def generate(self, filename: Path) -> None:
        """
        Generate the UMAT code for the material model and write it to a file.

        Args:
            filename (str): The file path where the UMAT code will be written.
        """
        props_code = self.generate_props_code()
        sigma_code = self.generate_expression(self.cauchy, "stress")
        smat_code = self.generate_expression(self.tangent, "ddsdde")
        props_code_str = self.code_as_string(props_code)
        sigma_code_str = self.code_as_string(sigma_code)
        smat_code_str = self.code_as_string(smat_code)
        self.write_umat_code(props_code_str, sigma_code_str, smat_code_str, filename)

    @property
    def cauchy(self) -> sym.Matrix:
        """
        Generate the symbolic expression for the Cauchy stress tensor.
        """
        logging.info("Generating Cauchy stress tensor...")
        return self.material.cauchy_voigt(self.f).subs(self.sub_exp)

    @property
    def tangent(self) -> sym.Matrix:
        """
        Generate the symbolic expression for the tangent matrix.
        """
        logging.info("Generating tangent matrix...")
        return self.material.tangent_voigt(self.f, use_jaumann_rate=True).subs(self.sub_exp)

    def generate_props_code(self) -> list[str]:
        """
        Generate the Fortran code for material properties.

        Returns:
            list: The Fortran code for material properties.
        """
        logging.info("Generating material properties code...")
        parameters_names = list(self.material._params.keys())
        props_init = [f"DOUBLE PRECISION :: {param}" for param in parameters_names]
        props_code = [f"{param} = PROPS({i + 1})" for i, param in enumerate(parameters_names)]
        return props_init + props_code

    def generate_expression(self, tensor: sym.Matrix, var_name: str) -> list[str]:
        logging.info(f"Generating CSE for {var_name}...")
        return self.common_subexpressions(tensor, var_name)

    @staticmethod
    def code_as_string(code: list) -> str:
        """
        Convert a list of code lines into a single string.

        Args:
            code (list): The list of code lines.

        Returns:
            str: The code as a single string.
        """
        return "\n".join(code)

    def write_umat_code(self, props_code_str: str, sigma_code_str: str, smat_code_str: str, filename: Path) -> None:
        """
        Write the generated Fortran code into a UMAT subroutine file.

        Args:
            filename (Path): The file path where the UMAT code will be written.
        """
        today = datetime.datetime.now().strftime("%Y-%m-%d")
        description = f"Automatically generated code for the UMAT subroutine using {self.material.__class__.__name__}."
        umat_code = f"""
!>********************************************************************
!> Record of revisions:                                              |
!>        Date        Programmer        Description of change        |
!>        ====        ==========        =====================        |
!>     {today}    Automatic Code      {description}                  |
!>--------------------------------------------------------------------
!C>
!C>   Material model: {self.material.__class__.__name__}
!C>
!---------------------------------------------------------------------

SUBROUTINE umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, drplde, drpldt,  &
    stran, dstran, time, dtime, temp, dtemp, predef, dpred, cmname,  &
    ndi, nshr, ntens, nstatev, props, nprops, coords, drot, pnewdt,  &
    celent, dfgrd0, dfgrd1, noel, npt, layer, kspt, kstep, kinc)
!
!use global
IMPLICIT DOUBLE PRECISION (x)
!----------------------------------------------------------------------
!--------------------------- DECLARATIONS -----------------------------
!----------------------------------------------------------------------
INTEGER, INTENT(IN OUT)                  :: noel
INTEGER, INTENT(IN OUT)                  :: npt
INTEGER, INTENT(IN OUT)                  :: layer
INTEGER, INTENT(IN OUT)                  :: kspt
INTEGER, INTENT(IN OUT)                  :: kstep
INTEGER, INTENT(IN OUT)                  :: kinc
INTEGER, INTENT(IN OUT)                  :: ndi
INTEGER, INTENT(IN OUT)                  :: nshr
INTEGER, INTENT(IN OUT)                  :: ntens
INTEGER, INTENT(IN OUT)                  :: nstatev
INTEGER, INTENT(IN OUT)                  :: nprops
DOUBLE PRECISION, INTENT(IN OUT)         :: sse
DOUBLE PRECISION, INTENT(IN OUT)         :: spd
DOUBLE PRECISION, INTENT(IN OUT)         :: scd
DOUBLE PRECISION, INTENT(IN OUT)         :: rpl
DOUBLE PRECISION, INTENT(IN OUT)         :: dtime
DOUBLE PRECISION, INTENT(IN OUT)         :: drpldt
DOUBLE PRECISION, INTENT(IN OUT)         :: temp
DOUBLE PRECISION, INTENT(IN OUT)         :: dtemp
CHARACTER (LEN=8), INTENT(IN OUT)        :: cmname
DOUBLE PRECISION, INTENT(IN OUT)         :: pnewdt
DOUBLE PRECISION, INTENT(IN OUT)         :: celent

DOUBLE PRECISION, INTENT(IN OUT)         :: stress(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: statev(nstatev)
DOUBLE PRECISION, INTENT(IN OUT)         :: ddsdde(ntens, ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: ddsddt(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: drplde(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: stran(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: dstran(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: time(2)
DOUBLE PRECISION, INTENT(IN OUT)         :: predef(1)
DOUBLE PRECISION, INTENT(IN OUT)         :: dpred(1)
DOUBLE PRECISION, INTENT(IN)             :: props(nprops)
DOUBLE PRECISION, INTENT(IN OUT)         :: coords(3)
DOUBLE PRECISION, INTENT(IN OUT)         :: drot(3, 3)
DOUBLE PRECISION, INTENT(IN OUT)         :: dfgrd0(3, 3)
DOUBLE PRECISION, INTENT(IN OUT)         :: dfgrd1(3, 3)

! Initialize material properties
{props_code_str}

! Define the stress calculation from the pk2 symbolic expression
{sigma_code_str}

! Define the tangent matrix calculation from the smat symbolic expression
{smat_code_str}

RETURN
END SUBROUTINE umat
"""

        with open(filename, "w") as file:
            file.write(umat_code)
        logging.info(f"UMAT subroutine written to {filename}")

cauchy property

Generate the symbolic expression for the Cauchy stress tensor.

f property

The deformation gradient tensor.

sub_exp property

Substitution expressions for the right Cauchy-Green tensor.

tangent property

Generate the symbolic expression for the tangent matrix.

__init__(material_model)

Initialize the UMAT handler with a specific material model.

Parameters:

Name	Type	Description	Default
material_model	Material	The material model (e.g., NeoHooke) to generate UMAT code for.	required

Source code in hyper_surrogate/export/fortran/analytical.py

def __init__(self, material_model: Material) -> None:
    """
    Initialize the UMAT handler with a specific material model.

    Args:
        material_model: The material model (e.g., NeoHooke) to generate UMAT code for.
    """
    self.material = material_model
    self.sigma_code = None
    self.smat_code = None

code_as_string(code) staticmethod

Convert a list of code lines into a single string.

Parameters:

Name	Type	Description	Default
code	list	The list of code lines.	required

Returns:

Name	Type	Description
str	str	The code as a single string.

Source code in hyper_surrogate/export/fortran/analytical.py

@staticmethod
def code_as_string(code: list) -> str:
    """
    Convert a list of code lines into a single string.

    Args:
        code (list): The list of code lines.

    Returns:
        str: The code as a single string.
    """
    return "\n".join(code)

common_subexpressions(tensor, var_name) staticmethod

Perform common subexpression elimination on a vector or matrix and generate Fortran code.

Parameters:

Name	Type	Description	Default
var_name	str	The base name for the variables in the Fortran code.	required

Returns:

Name	Type	Description
tuple	list[str]	A tuple containing Fortran code for auxiliary variables and reduced expressions.

Source code in hyper_surrogate/export/fortran/analytical.py

@staticmethod
def common_subexpressions(tensor: sym.Matrix, var_name: str) -> list[str]:
    """
    Perform common subexpression elimination on a vector or matrix and generate Fortran code.

    Args:
        var_name (str): The base name for the variables in the Fortran code.

    Returns:
        tuple: A tuple containing Fortran code for auxiliary variables and reduced expressions.
    """
    # Extract individual components
    # tensor_components = [tensor[i] for i in range(tensor.shape[0])]
    tensor_components = [tensor[i, j] for i in range(tensor.shape[0]) for j in range(tensor.shape[1])]
    # Convert to a matrix to check shape
    tensor_matrix = sym.Matrix(tensor)
    # Perform common subexpression elimination
    replacements, reduced_exprs = sym.cse(tensor_components)

    # Generate Fortran code for auxiliary variables (replacements)
    aux_code = [
        sym.fcode(
            expr,
            standard=90,
            source_format="free",
            assign_to=sym.fcode(var, standard=90, source_format="free"),
        )
        for var, expr in replacements
    ]

    # Generate Fortran code for reduced expressions
    if tensor_matrix.shape[1] == 1:  # If vector
        reduced_code = [
            sym.fcode(
                expr,
                standard=90,
                source_format="free",
                assign_to=f"{var_name}({i + 1})",
            )
            for i, expr in enumerate(reduced_exprs)
        ]
    else:  # If matrix
        _, cols = tensor.shape
        reduced_code = [
            sym.fcode(
                expr,
                standard=90,
                source_format="free",
                assign_to=f"{var_name}({i // cols + 1},{i % cols + 1})",
            )
            for i, expr in enumerate(reduced_exprs)
        ]

    return aux_code + reduced_code

generate(filename)

Generate the UMAT code for the material model and write it to a file.

Parameters:

Name	Type	Description	Default
filename	str	The file path where the UMAT code will be written.	required

Source code in hyper_surrogate/export/fortran/analytical.py

def generate(self, filename: Path) -> None:
    """
    Generate the UMAT code for the material model and write it to a file.

    Args:
        filename (str): The file path where the UMAT code will be written.
    """
    props_code = self.generate_props_code()
    sigma_code = self.generate_expression(self.cauchy, "stress")
    smat_code = self.generate_expression(self.tangent, "ddsdde")
    props_code_str = self.code_as_string(props_code)
    sigma_code_str = self.code_as_string(sigma_code)
    smat_code_str = self.code_as_string(smat_code)
    self.write_umat_code(props_code_str, sigma_code_str, smat_code_str, filename)

generate_props_code()

Generate the Fortran code for material properties.

Returns:

Name	Type	Description
list	list[str]	The Fortran code for material properties.

Source code in hyper_surrogate/export/fortran/analytical.py

def generate_props_code(self) -> list[str]:
    """
    Generate the Fortran code for material properties.

    Returns:
        list: The Fortran code for material properties.
    """
    logging.info("Generating material properties code...")
    parameters_names = list(self.material._params.keys())
    props_init = [f"DOUBLE PRECISION :: {param}" for param in parameters_names]
    props_code = [f"{param} = PROPS({i + 1})" for i, param in enumerate(parameters_names)]
    return props_init + props_code

write_umat_code(props_code_str, sigma_code_str, smat_code_str, filename)

Write the generated Fortran code into a UMAT subroutine file.

Parameters:

Name	Type	Description	Default
filename	Path	The file path where the UMAT code will be written.	required

Source code in hyper_surrogate/export/fortran/analytical.py

    def write_umat_code(self, props_code_str: str, sigma_code_str: str, smat_code_str: str, filename: Path) -> None:
        """
        Write the generated Fortran code into a UMAT subroutine file.

        Args:
            filename (Path): The file path where the UMAT code will be written.
        """
        today = datetime.datetime.now().strftime("%Y-%m-%d")
        description = f"Automatically generated code for the UMAT subroutine using {self.material.__class__.__name__}."
        umat_code = f"""
!>********************************************************************
!> Record of revisions:                                              |
!>        Date        Programmer        Description of change        |
!>        ====        ==========        =====================        |
!>     {today}    Automatic Code      {description}                  |
!>--------------------------------------------------------------------
!C>
!C>   Material model: {self.material.__class__.__name__}
!C>
!---------------------------------------------------------------------

SUBROUTINE umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, drplde, drpldt,  &
    stran, dstran, time, dtime, temp, dtemp, predef, dpred, cmname,  &
    ndi, nshr, ntens, nstatev, props, nprops, coords, drot, pnewdt,  &
    celent, dfgrd0, dfgrd1, noel, npt, layer, kspt, kstep, kinc)
!
!use global
IMPLICIT DOUBLE PRECISION (x)
!----------------------------------------------------------------------
!--------------------------- DECLARATIONS -----------------------------
!----------------------------------------------------------------------
INTEGER, INTENT(IN OUT)                  :: noel
INTEGER, INTENT(IN OUT)                  :: npt
INTEGER, INTENT(IN OUT)                  :: layer
INTEGER, INTENT(IN OUT)                  :: kspt
INTEGER, INTENT(IN OUT)                  :: kstep
INTEGER, INTENT(IN OUT)                  :: kinc
INTEGER, INTENT(IN OUT)                  :: ndi
INTEGER, INTENT(IN OUT)                  :: nshr
INTEGER, INTENT(IN OUT)                  :: ntens
INTEGER, INTENT(IN OUT)                  :: nstatev
INTEGER, INTENT(IN OUT)                  :: nprops
DOUBLE PRECISION, INTENT(IN OUT)         :: sse
DOUBLE PRECISION, INTENT(IN OUT)         :: spd
DOUBLE PRECISION, INTENT(IN OUT)         :: scd
DOUBLE PRECISION, INTENT(IN OUT)         :: rpl
DOUBLE PRECISION, INTENT(IN OUT)         :: dtime
DOUBLE PRECISION, INTENT(IN OUT)         :: drpldt
DOUBLE PRECISION, INTENT(IN OUT)         :: temp
DOUBLE PRECISION, INTENT(IN OUT)         :: dtemp
CHARACTER (LEN=8), INTENT(IN OUT)        :: cmname
DOUBLE PRECISION, INTENT(IN OUT)         :: pnewdt
DOUBLE PRECISION, INTENT(IN OUT)         :: celent

DOUBLE PRECISION, INTENT(IN OUT)         :: stress(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: statev(nstatev)
DOUBLE PRECISION, INTENT(IN OUT)         :: ddsdde(ntens, ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: ddsddt(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: drplde(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: stran(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: dstran(ntens)
DOUBLE PRECISION, INTENT(IN OUT)         :: time(2)
DOUBLE PRECISION, INTENT(IN OUT)         :: predef(1)
DOUBLE PRECISION, INTENT(IN OUT)         :: dpred(1)
DOUBLE PRECISION, INTENT(IN)             :: props(nprops)
DOUBLE PRECISION, INTENT(IN OUT)         :: coords(3)
DOUBLE PRECISION, INTENT(IN OUT)         :: drot(3, 3)
DOUBLE PRECISION, INTENT(IN OUT)         :: dfgrd0(3, 3)
DOUBLE PRECISION, INTENT(IN OUT)         :: dfgrd1(3, 3)

! Initialize material properties
{props_code_str}

! Define the stress calculation from the pk2 symbolic expression
{sigma_code_str}

! Define the tangent matrix calculation from the smat symbolic expression
{smat_code_str}

RETURN
END SUBROUTINE umat
"""

        with open(filename, "w") as file:
            file.write(umat_code)
        logging.info(f"UMAT subroutine written to {filename}")

Reporting

Reporter

Generate a PDF report with deformation diagnostics.

Accepts a batch of (N, 3, 3) tensors — either deformation gradients F or right Cauchy-Green tensors C — and produces histograms of key continuum-mechanics quantities (eigenvalues, determinants, invariants, principal stretches).

Parameters:

Name	Type	Description	Default
tensor	ndarray	Array of shape (N, 3, 3).	required
tensor_type	str	"C" (right Cauchy-Green, default) or "F" (deformation gradient).	'C'

Example::

from hyper_surrogate.reporting.reporter import Reporter

reporter = Reporter(C)  # (N, 3, 3)
reporter.fig_eigenvalues()
reporter.fig_determinants()
reporter.fig_invariants()
reporter.fig_principal_stretches()
reporter.generate_report("report/")

Source code in hyper_surrogate/reporting/reporter.py

class Reporter:
    """Generate a PDF report with deformation diagnostics.

    Accepts a batch of (N, 3, 3) tensors — either deformation gradients **F**
    or right Cauchy-Green tensors **C** — and produces histograms of key
    continuum-mechanics quantities (eigenvalues, determinants, invariants,
    principal stretches).

    Args:
        tensor: Array of shape ``(N, 3, 3)``.
        tensor_type: ``"C"`` (right Cauchy-Green, default) or ``"F"``
            (deformation gradient).

    Example::

        from hyper_surrogate.reporting.reporter import Reporter

        reporter = Reporter(C)  # (N, 3, 3)
        reporter.fig_eigenvalues()
        reporter.fig_determinants()
        reporter.fig_invariants()
        reporter.fig_principal_stretches()
        reporter.generate_report("report/")
    """

    LAYOUT: ClassVar[list[str]] = ["standalone", "combined"]

    FIG_SIZE: ClassVar[tuple[float, float]] = (8.27, 11.69)

    REPORT_FIGURES: ClassVar[list[str]] = [
        "fig_eigenvalues",
        "fig_determinants",
        "fig_invariants",
        "fig_principal_stretches",
        "fig_volume_change",
    ]

    def __init__(
        self,
        tensor: np.ndarray,
        tensor_type: str = "C",
    ) -> None:
        if tensor.ndim != 3 or tensor.shape[1:] != (3, 3):
            msg = f"Expected tensor of shape (N, 3, 3), got {tensor.shape}"
            raise ValueError(msg)

        self.tensor_type = tensor_type.upper()
        self.F: np.ndarray | None = None
        if self.tensor_type == "F":
            self.F = tensor
            self.C = Kinematics.right_cauchy_green(tensor)
        elif self.tensor_type == "C":
            self.C = tensor
        else:
            msg = f"tensor_type must be 'C' or 'F', got {tensor_type!r}"
            raise ValueError(msg)

    # ------------------------------------------------------------------
    # Convenience properties
    # ------------------------------------------------------------------

    @property
    def n_samples(self) -> int:
        """Number of samples in the batch."""
        return int(self.C.shape[0])

    # ------------------------------------------------------------------
    # Statistics
    # ------------------------------------------------------------------

    def basic_statistics(self) -> dict[str, dict[str, float]]:
        """Per-quantity summary statistics.

        Returns a dict keyed by quantity name, each containing
        ``mean``, ``std``, ``min``, ``max``.
        """
        quantities: dict[str, np.ndarray] = {
            "det(C)": Kinematics.det_invariant(self.C),
            "I1_bar": Kinematics.isochoric_invariant1(self.C),
            "I2_bar": Kinematics.isochoric_invariant2(self.C),
            "J": np.sqrt(Kinematics.det_invariant(self.C)),
        }
        stats: dict[str, dict[str, float]] = {}
        for name, values in quantities.items():
            stats[name] = {
                "mean": float(np.mean(values)),
                "std": float(np.std(values)),
                "min": float(np.min(values)),
                "max": float(np.max(values)),
            }
        return stats

    # ------------------------------------------------------------------
    # Individual figure methods
    # ------------------------------------------------------------------

    def fig_eigenvalues(self) -> list[matplotlib.figure.Figure]:
        """Histogram of eigenvalues of C."""
        eigenvalues = np.linalg.eigvalsh(self.C).real
        fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
        labels = [r"$\lambda_1^2$", r"$\lambda_2^2$", r"$\lambda_3^2$"]
        for i, ax in enumerate(axes):
            ax.hist(eigenvalues[:, i], bins="auto", alpha=0.75, edgecolor="black")
            ax.set_title(f"Eigenvalue {labels[i]}")
            ax.set_xlabel("Value")
            ax.set_ylabel("Frequency")
        fig.suptitle("Eigenvalues of C", fontsize=14)
        fig.tight_layout()
        return [fig]

    def fig_determinants(self) -> list[matplotlib.figure.Figure]:
        """Histogram of det(C)."""
        determinants = Kinematics.det_invariant(self.C)
        fig, ax = plt.subplots(1, 1, figsize=self.FIG_SIZE)
        ax.hist(determinants, bins="auto", alpha=0.75, edgecolor="black")
        ax.set_title("Determinant of C")
        ax.set_xlabel("det(C)")
        ax.set_ylabel("Frequency")
        fig.tight_layout()
        return [fig]

    def fig_invariants(self) -> list[matplotlib.figure.Figure]:
        """Histograms of isochoric invariants I1_bar, I2_bar, and J."""
        i1 = Kinematics.isochoric_invariant1(self.C)
        i2 = Kinematics.isochoric_invariant2(self.C)
        j = np.sqrt(Kinematics.det_invariant(self.C))

        fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
        for ax, data, label in zip(
            axes,
            [i1, i2, j],
            [r"$\bar{I}_1$", r"$\bar{I}_2$", r"$J$"],
            strict=False,
        ):
            ax.hist(data, bins="auto", alpha=0.75, edgecolor="black")
            ax.set_title(label)
            ax.set_xlabel("Value")
            ax.set_ylabel("Frequency")
        fig.suptitle("Isochoric Invariants", fontsize=14)
        fig.tight_layout()
        return [fig]

    def fig_principal_stretches(self) -> list[matplotlib.figure.Figure]:
        """Histogram of principal stretches (sorted)."""
        eigenvalues = np.sort(np.linalg.eigvalsh(self.C).real, axis=1)
        stretches = np.sqrt(np.maximum(eigenvalues, 0.0))

        fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
        for i, ax in enumerate(axes):
            ax.hist(stretches[:, i], bins="auto", alpha=0.75, edgecolor="black")
            ax.set_title(rf"$\lambda_{{{i + 1}}}$")
            ax.set_xlabel("Stretch")
            ax.set_ylabel("Frequency")
        fig.suptitle("Principal Stretches", fontsize=14)
        fig.tight_layout()
        return [fig]

    def fig_volume_change(self) -> list[matplotlib.figure.Figure]:
        """Histogram of volume ratio J = sqrt(det(C))."""
        j = np.sqrt(Kinematics.det_invariant(self.C))
        fig, ax = plt.subplots(1, 1, figsize=self.FIG_SIZE)
        ax.hist(j, bins="auto", alpha=0.75, edgecolor="black")
        ax.axvline(1.0, color="red", linestyle="--", linewidth=1, label="J = 1")
        ax.set_title("Volume Ratio J")
        ax.set_xlabel("J")
        ax.set_ylabel("Frequency")
        ax.legend()
        fig.tight_layout()
        return [fig]

    # ------------------------------------------------------------------
    # Report generation
    # ------------------------------------------------------------------

    def generate_figures(self) -> list[matplotlib.figure.Figure]:
        """Generate all report figures."""
        fig_list: list[matplotlib.figure.Figure] = []
        for name in self.REPORT_FIGURES:
            func = getattr(self, name)
            fig_list.extend(func())
        return fig_list

    def generate_report(self, save_dir: str | Path, layout: str = "combined") -> None:
        """Create a PDF report.

        Args:
            save_dir: Directory to write the report into (created if needed).
            layout: ``"combined"`` (single PDF) or ``"standalone"`` (one PDF
                per figure).
        """
        save_path = Path(save_dir)
        save_path.mkdir(parents=True, exist_ok=True)

        metadata = {
            "Title": "Deformation Report",
            "Subject": f"Batch of {self.n_samples} tensors",
            "CreationDate": datetime.today(),
        }

        fig_list = self.generate_figures()

        if layout == "combined":
            with PdfPages(save_path / "report.pdf", metadata=metadata) as pp:
                for fig in fig_list:
                    fig.savefig(pp, format="pdf", bbox_inches="tight")
        else:
            for fig in fig_list:
                title = fig.axes[0].get_title() or "figure"
                safe_name = title.replace(" ", "_").replace("/", "_")
                with PdfPages(save_path / f"{safe_name}.pdf", metadata=metadata) as pp:
                    fig.savefig(pp, format="pdf", bbox_inches="tight")
        plt.close("all")

    # Keep backward compat alias
    create_report = generate_report

n_samples property

Number of samples in the batch.

basic_statistics()

Per-quantity summary statistics.

Returns a dict keyed by quantity name, each containing mean, std, min, max.

Source code in hyper_surrogate/reporting/reporter.py

def basic_statistics(self) -> dict[str, dict[str, float]]:
    """Per-quantity summary statistics.

    Returns a dict keyed by quantity name, each containing
    ``mean``, ``std``, ``min``, ``max``.
    """
    quantities: dict[str, np.ndarray] = {
        "det(C)": Kinematics.det_invariant(self.C),
        "I1_bar": Kinematics.isochoric_invariant1(self.C),
        "I2_bar": Kinematics.isochoric_invariant2(self.C),
        "J": np.sqrt(Kinematics.det_invariant(self.C)),
    }
    stats: dict[str, dict[str, float]] = {}
    for name, values in quantities.items():
        stats[name] = {
            "mean": float(np.mean(values)),
            "std": float(np.std(values)),
            "min": float(np.min(values)),
            "max": float(np.max(values)),
        }
    return stats

fig_determinants()

Histogram of det(C).

Source code in hyper_surrogate/reporting/reporter.py

def fig_determinants(self) -> list[matplotlib.figure.Figure]:
    """Histogram of det(C)."""
    determinants = Kinematics.det_invariant(self.C)
    fig, ax = plt.subplots(1, 1, figsize=self.FIG_SIZE)
    ax.hist(determinants, bins="auto", alpha=0.75, edgecolor="black")
    ax.set_title("Determinant of C")
    ax.set_xlabel("det(C)")
    ax.set_ylabel("Frequency")
    fig.tight_layout()
    return [fig]

fig_eigenvalues()

Histogram of eigenvalues of C.

Source code in hyper_surrogate/reporting/reporter.py

def fig_eigenvalues(self) -> list[matplotlib.figure.Figure]:
    """Histogram of eigenvalues of C."""
    eigenvalues = np.linalg.eigvalsh(self.C).real
    fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
    labels = [r"$\lambda_1^2$", r"$\lambda_2^2$", r"$\lambda_3^2$"]
    for i, ax in enumerate(axes):
        ax.hist(eigenvalues[:, i], bins="auto", alpha=0.75, edgecolor="black")
        ax.set_title(f"Eigenvalue {labels[i]}")
        ax.set_xlabel("Value")
        ax.set_ylabel("Frequency")
    fig.suptitle("Eigenvalues of C", fontsize=14)
    fig.tight_layout()
    return [fig]

fig_invariants()

Histograms of isochoric invariants I1_bar, I2_bar, and J.

Source code in hyper_surrogate/reporting/reporter.py

def fig_invariants(self) -> list[matplotlib.figure.Figure]:
    """Histograms of isochoric invariants I1_bar, I2_bar, and J."""
    i1 = Kinematics.isochoric_invariant1(self.C)
    i2 = Kinematics.isochoric_invariant2(self.C)
    j = np.sqrt(Kinematics.det_invariant(self.C))

    fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
    for ax, data, label in zip(
        axes,
        [i1, i2, j],
        [r"$\bar{I}_1$", r"$\bar{I}_2$", r"$J$"],
        strict=False,
    ):
        ax.hist(data, bins="auto", alpha=0.75, edgecolor="black")
        ax.set_title(label)
        ax.set_xlabel("Value")
        ax.set_ylabel("Frequency")
    fig.suptitle("Isochoric Invariants", fontsize=14)
    fig.tight_layout()
    return [fig]

fig_principal_stretches()

Histogram of principal stretches (sorted).

Source code in hyper_surrogate/reporting/reporter.py

def fig_principal_stretches(self) -> list[matplotlib.figure.Figure]:
    """Histogram of principal stretches (sorted)."""
    eigenvalues = np.sort(np.linalg.eigvalsh(self.C).real, axis=1)
    stretches = np.sqrt(np.maximum(eigenvalues, 0.0))

    fig, axes = plt.subplots(1, 3, figsize=self.FIG_SIZE)
    for i, ax in enumerate(axes):
        ax.hist(stretches[:, i], bins="auto", alpha=0.75, edgecolor="black")
        ax.set_title(rf"$\lambda_{{{i + 1}}}$")
        ax.set_xlabel("Stretch")
        ax.set_ylabel("Frequency")
    fig.suptitle("Principal Stretches", fontsize=14)
    fig.tight_layout()
    return [fig]

fig_volume_change()

Histogram of volume ratio J = sqrt(det(C)).

Source code in hyper_surrogate/reporting/reporter.py

def fig_volume_change(self) -> list[matplotlib.figure.Figure]:
    """Histogram of volume ratio J = sqrt(det(C))."""
    j = np.sqrt(Kinematics.det_invariant(self.C))
    fig, ax = plt.subplots(1, 1, figsize=self.FIG_SIZE)
    ax.hist(j, bins="auto", alpha=0.75, edgecolor="black")
    ax.axvline(1.0, color="red", linestyle="--", linewidth=1, label="J = 1")
    ax.set_title("Volume Ratio J")
    ax.set_xlabel("J")
    ax.set_ylabel("Frequency")
    ax.legend()
    fig.tight_layout()
    return [fig]

generate_figures()

Generate all report figures.

Source code in hyper_surrogate/reporting/reporter.py

def generate_figures(self) -> list[matplotlib.figure.Figure]:
    """Generate all report figures."""
    fig_list: list[matplotlib.figure.Figure] = []
    for name in self.REPORT_FIGURES:
        func = getattr(self, name)
        fig_list.extend(func())
    return fig_list

generate_report(save_dir, layout='combined')

Create a PDF report.

Parameters:

Name	Type	Description	Default
save_dir	str | Path	Directory to write the report into (created if needed).	required
layout	str	"combined" (single PDF) or "standalone" (one PDF per figure).	'combined'

Source code in hyper_surrogate/reporting/reporter.py

def generate_report(self, save_dir: str | Path, layout: str = "combined") -> None:
    """Create a PDF report.

    Args:
        save_dir: Directory to write the report into (created if needed).
        layout: ``"combined"`` (single PDF) or ``"standalone"`` (one PDF
            per figure).
    """
    save_path = Path(save_dir)
    save_path.mkdir(parents=True, exist_ok=True)

    metadata = {
        "Title": "Deformation Report",
        "Subject": f"Batch of {self.n_samples} tensors",
        "CreationDate": datetime.today(),
    }

    fig_list = self.generate_figures()

    if layout == "combined":
        with PdfPages(save_path / "report.pdf", metadata=metadata) as pp:
            for fig in fig_list:
                fig.savefig(pp, format="pdf", bbox_inches="tight")
    else:
        for fig in fig_list:
            title = fig.axes[0].get_title() or "figure"
            safe_name = title.replace(" ", "_").replace("/", "_")
            with PdfPages(save_path / f"{safe_name}.pdf", metadata=metadata) as pp:
                fig.savefig(pp, format="pdf", bbox_inches="tight")
    plt.close("all")