# Copyright (C) 2021-2025 Université Gustave Eiffel.
# This file is part of the EasyFEA project.
# EasyFEA is distributed under the terms of the GNU General Public License v3, see LICENSE.txt and CREDITS.md for more information.
"""Hyperelastic laws."""
from abc import ABC, abstractmethod
import numpy as np
from typing import Union
# utilities
from ..fem import Mesh, MatrixType, FeArray, TensorProd
from ._hyperelastic import HyperElasticState
# others
from ._utils import _IModel, ModelType
from ..utilities import _params, _types
# ----------------------------------------------
# Hyper Elastic
# ----------------------------------------------
class _HyperElas(_IModel, ABC):
"""HyperElastic material.\n
NeoHookean, MooneyRivlin and SaintVenantKirchhoff inherit from _HyperElas class.
"""
dim: float = _params.ParameterInValues([1, 2, 3])
thickness: float = _params.PositiveScalarParameter()
def __init__(self, dim: int, thickness: float):
self.dim = dim
self.thickness = thickness
@property
def modelType(self) -> ModelType:
return ModelType.hyperelastic
@property
def coef(self) -> float:
"""kelvin mandel coef -> sqrt(2)"""
return np.sqrt(2)
# Model
@staticmethod
def Available_Laws():
laws = [NeoHookean, MooneyRivlin, SaintVenantKirchhoff]
return laws
@property
def isHeterogeneous(self) -> bool:
return False
@abstractmethod
def Compute_W(
self, mesh: Mesh, u: _types.FloatArray, matrixType=MatrixType.rigi
) -> FeArray:
"""Computes the quadratic energy W(u).
Parameters
----------
mesh : Mesh
mesh
u : _types.FloatArray
discretized displacement field [u1, v1, w1, . . ., uN, vN, wN] of size Nn * dim
matrixType : MatrixType, optional
matrix type, by default MatrixType.rigi
Returns
-------
FeArray
We_e_pg of shape (Ne, pg)
"""
return None # type: ignore [return-value]
@abstractmethod
def Compute_dWde(
self, mesh: Mesh, u: _types.FloatArray, matrixType=MatrixType.rigi
) -> FeArray:
"""Computes the second Piola-Kirchhoff tensor Σ(u).
Parameters
----------
mesh : Mesh
mesh
u : _types.FloatArray
discretized displacement field [u1, v1, w1, . . ., uN, vN, wN] of size Nn * dim
matrixType : MatrixType, optional
matrix type, by default MatrixType.rigi
Returns
-------
FeArray
Σ_e_pg of shape (Ne, pg, 6)
Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy
"""
return None # type: ignore [return-value]
@abstractmethod
def Compute_d2Wde(
self, mesh: Mesh, u: _types.FloatArray, matrixType=MatrixType.rigi
) -> FeArray:
"""Computes dΣde.
Parameters
----------
mesh : Mesh
mesh
u : _types.FloatArray
discretized displacement field [u1, v1, w1, . . ., uN, vN, wN] of size Nn * dim
matrixType : MatrixType, optional
matrix type, by default MatrixType.rigi
Returns
-------
FeArray
dΣde_e_pg of shape (Ne, pg, 6, 6)
"""
return None # type: ignore [return-value]
# ----------------------------------------------
# Neo-Hookean
# ----------------------------------------------
[docs]
class NeoHookean(_HyperElas):
K: float = _params.PositiveScalarParameter()
"""Bulk modulus"""
def __init__(self, dim: int, K: Union[float, _types.FloatArray], thickness=1.0):
"""Creates an Neo-Hookean material.
Parameters
----------
dim : int
dimension (e.g 2 or 3)
K : float|_types.FloatArray
Bulk modulus
thickness : float, optional
thickness, by default 1.0
"""
_HyperElas.__init__(self, dim, thickness)
self.K = K
[docs]
def Compute_W(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I3 = hyperElasticState.Compute_I3()
W = K * (I1 / I3 ** (1 / 3) - 3)
return W
[docs]
def Compute_dWde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I3 = hyperElasticState.Compute_I3()
dWdI1 = K / I3 ** (1 / 3)
dWdI3 = -I1 * K / (3 * I3 ** (4 / 3))
dI1dC = hyperElasticState.Compute_dI1dC()
dI3dC = hyperElasticState.Compute_dI3dC()
dWdI1 = K / I3 ** (1 / 3)
dWdI3 = -I1 * K / (3 * I3 ** (4 / 3))
dW = 2 * (dWdI1 * dI1dC + dWdI3 * dI3dC)
return dW
[docs]
def Compute_d2Wde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I3 = hyperElasticState.Compute_I3()
dI1dC = hyperElasticState.Compute_dI1dC()
dI3dC = hyperElasticState.Compute_dI3dC()
d2I1dC = hyperElasticState.Compute_d2I1dC()
d2I3dC = hyperElasticState.Compute_d2I3dC()
dWdI1 = K / I3 ** (1 / 3)
d2WdI1dI3 = -K / (3 * I3 ** (4 / 3))
dWdI3 = -I1 * K / (3 * I3 ** (4 / 3))
d2WdI3dI1 = -K / (3 * I3 ** (4 / 3))
d2WdI3dI3 = 4 * I1 * K / (9 * I3 ** (7 / 3))
d2W = 4 * (dWdI1 * d2I1dC + dWdI3 * d2I3dC) + 4 * (
d2WdI1dI3 * TensorProd(dI1dC, dI3dC)
+ d2WdI3dI1 * TensorProd(dI3dC, dI1dC)
+ d2WdI3dI3 * TensorProd(dI3dC, dI3dC)
)
return d2W
# ----------------------------------------------
# Mooney-Rivlin
# ----------------------------------------------
[docs]
class MooneyRivlin(_HyperElas):
K1: float = _params.PositiveScalarParameter()
"""Kappa1"""
K2: float = _params.PositiveScalarParameter()
"""Kappa2"""
K: float = _params.PositiveScalarParameter()
"""Bulk modulus"""
def __init__(
self,
dim: int,
K1: Union[float, _types.FloatArray],
K2: Union[float, _types.FloatArray],
K: Union[float, _types.FloatArray] = 0.0,
thickness=1.0,
):
"""Creates an Mooney-Rivlin material.
Parameters
----------
dim : int
dimension (e.g 2 or 3)
K1 : float|_types.FloatArray
Kappa1
K2 : float|_types.FloatArray
Kappa2 -> Neo-Hoolkean if K2=0
K : float|_types.FloatArray, optional
Bulk modulus, by default 0.0
thickness : float, optional
thickness, by default 1.0
"""
_HyperElas.__init__(self, dim, thickness)
self.K1 = K1
self.K2 = K2
self.K = K
[docs]
def Compute_W(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
K1 = self.K1
K2 = self.K2
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
W = (
K * (np.sqrt(I3) - 1) ** 2
+ K1 * (I1 / I3 ** (1 / 3) - 3)
+ K2 * (I2 / I3 ** (2 / 3) - 3)
)
return W
[docs]
def Compute_dWde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
K1 = self.K1
K2 = self.K2
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
dI1dC = hyperElasticState.Compute_dI1dC()
dI2dC = hyperElasticState.Compute_dI2dC()
dI3dC = hyperElasticState.Compute_dI3dC()
dWdI1 = K1 / I3 ** (1 / 3)
dWdI2 = K2 / I3 ** (2 / 3)
dWdI3 = (
-I1 * K1 / (3 * I3 ** (4 / 3))
- 2 * I2 * K2 / (3 * I3 ** (5 / 3))
+ K * (np.sqrt(I3) - 1) / np.sqrt(I3)
)
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC)
return dW
[docs]
def Compute_d2Wde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
K = self.K
K1 = self.K1
K2 = self.K2
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
dI1dC = hyperElasticState.Compute_dI1dC()
dI2dC = hyperElasticState.Compute_dI2dC()
dI3dC = hyperElasticState.Compute_dI3dC()
d2I1dC = hyperElasticState.Compute_d2I1dC()
d2I2dC = hyperElasticState.Compute_d2I2dC()
d2I3dC = hyperElasticState.Compute_d2I3dC()
dWdI1 = K1 / I3 ** (1 / 3)
d2WdI1dI3 = -K1 / (3 * I3 ** (4 / 3))
dWdI2 = K2 / I3 ** (2 / 3)
d2WdI2dI3 = -2 * K2 / (3 * I3 ** (5 / 3))
dWdI3 = (
-I1 * K1 / (3 * I3 ** (4 / 3))
- 2 * I2 * K2 / (3 * I3 ** (5 / 3))
+ K * (np.sqrt(I3) - 1) / np.sqrt(I3)
)
d2WdI3dI1 = -K1 / (3 * I3 ** (4 / 3))
d2WdI3dI2 = -2 * K2 / (3 * I3 ** (5 / 3))
d2WdI3dI3 = (
4 * I1 * K1 / (9 * I3 ** (7 / 3))
+ 10 * I2 * K2 / (9 * I3 ** (8 / 3))
+ K / (2 * I3)
- K * (np.sqrt(I3) - 1) / (2 * I3 ** (3 / 2))
)
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC) + 4 * (
d2WdI1dI3 * TensorProd(dI1dC, dI3dC)
+ d2WdI2dI3 * TensorProd(dI2dC, dI3dC)
+ d2WdI3dI1 * TensorProd(dI3dC, dI1dC)
+ d2WdI3dI2 * TensorProd(dI3dC, dI2dC)
+ d2WdI3dI3 * TensorProd(dI3dC, dI3dC)
)
return d2W
# ----------------------------------------------
# Saint-Venant-Kirchhoff
# ----------------------------------------------
[docs]
class SaintVenantKirchhoff(_HyperElas):
lmbda: float = _params.ScalarParameter()
"""Lame's first parameter"""
mu: float = _params.PositiveScalarParameter()
"""Shear modulus"""
K: float = _params.PositiveScalarParameter()
"""Bulk modulus"""
def __init__(
self,
dim: int,
lmbda: Union[float, _types.FloatArray],
mu: Union[float, _types.FloatArray],
K: Union[float, _types.FloatArray] = 0.0,
thickness=1.0,
):
"""Creates Saint-Venant-Kirchhoff material.
Parameters
----------
dim : int
dimension (e.g 2 or 3)
lmbda : float|_types.FloatArray
Lame's first parameter
mu : float|_types.FloatArray
Shear modulus
K : float|_types.FloatArray, optional
Bulk modulus, by default 0.0
thickness : float, optional
thickness, by default 1.0
"""
_HyperElas.__init__(self, dim, thickness)
self.lmbda = lmbda
self.mu = mu
self.K = K
[docs]
def Compute_W(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
lmbda = self.lmbda
mu = self.mu
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
W = (
I1**2 * (lmbda / 8 + mu / 4)
- I1 * (3 * lmbda / 4 + mu / 2)
- I2 * mu / 2
+ 0.5 * K * (I3 - 1) ** 2
+ 9 * lmbda / 8
+ 3 * mu / 4
)
return W
[docs]
def Compute_dWde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
lmbda = self.lmbda
mu = self.mu
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I3 = hyperElasticState.Compute_I3()
dI1dC = hyperElasticState.Compute_dI1dC()
dI2dC = hyperElasticState.Compute_dI2dC()
dI3dC = hyperElasticState.Compute_dI3dC()
dWdI1 = 2 * I1 * (lmbda / 8 + mu / 4) - 3 * lmbda / 4 - mu / 2
dWdI2 = -mu / 2
dWdI3 = 0.5 * K * (2 * I3 - 2)
dW = 2 * (dWdI1 * dI1dC + dWdI2 * dI2dC + dWdI3 * dI3dC)
return dW
[docs]
def Compute_d2Wde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
lmbda = self.lmbda
mu = self.mu
K = self.K
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I3 = hyperElasticState.Compute_I3()
dI1dC = hyperElasticState.Compute_dI1dC()
dI3dC = hyperElasticState.Compute_dI3dC()
d2I1dC = hyperElasticState.Compute_d2I1dC()
d2I2dC = hyperElasticState.Compute_d2I2dC()
d2I3dC = hyperElasticState.Compute_d2I3dC()
dWdI1 = 2 * I1 * (lmbda / 8 + mu / 4) - 3 * lmbda / 4 - mu / 2
dWdI2 = -mu / 2
dWdI3 = 0.5 * K * (2 * I3 - 2)
d2WdI1dI1 = lmbda / 4 + mu / 2
d2WdI3dI3 = 1.0 * K
d2W = 4 * (dWdI1 * d2I1dC + dWdI2 * d2I2dC + dWdI3 * d2I3dC) + 4 * (
d2WdI1dI1 * TensorProd(dI1dC, dI1dC) + d2WdI3dI3 * TensorProd(dI3dC, dI3dC)
)
return d2W
# ----------------------------------------------
# Holzapfel-Ogden
# ----------------------------------------------
[docs]
class HolzapfelOgden(_HyperElas):
C0: float = _params.PositiveScalarParameter()
C1: float = _params.PositiveScalarParameter()
C2: float = _params.PositiveScalarParameter()
C3: float = _params.PositiveScalarParameter()
C4: float = _params.PositiveScalarParameter()
C5: float = _params.PositiveScalarParameter()
C6: float = _params.PositiveScalarParameter()
C7: float = _params.PositiveScalarParameter()
K: float = _params.PositiveScalarParameter()
"""Bulk modulus"""
Mu1: float = _params.PositiveScalarParameter()
Mu2: float = _params.PositiveScalarParameter()
T1 = _params.VectorParameter()
"""direction(s) 1, used for the invariants I4 and I8"""
T2 = _params.VectorParameter()
"""direction(s) 2, used for the invariants I6 and I8"""
__ks: float = _params.PositiveScalarParameter()
"""A positive constant used in the incompressibility penalty term, as proposed in http://dx.doi.org/10.1016/0045-7825(94)90051-5."""
def __init__(
self,
dim: int,
C0: float,
C1: float,
C2: float,
C3: float,
C4: float,
C5: float,
C6: float,
C7: float,
K: float,
Mu1: float,
Mu2: float,
T1: _types.FloatArray,
T2: _types.FloatArray,
ks: float = 100,
thickness=1.0,
):
"""Creates Holzapfel-Ogden material.
Parameters
----------
dim : int
dimension (e.g 2 or 3)
C0 : float
C0
C1 : float
C1
C2 : float
C2
C3 : float
C3
C4 : float
C4
C5 : float
C5
C6 : float
C6
C7 : float
C7
K : float
bulk modulus
Mu1 : float
Mu1
Mu2 : float
Mu2
T1 : _type.FloatArray
direction(s) 1, used for the invariants I4 and I8
T2 : _type.FloatArray
direction(s) 2, used for the invariants I6 and I8
Mu2 : float
Mu2
ks : float, optional
A positive constant used in the incompressibility penalty term.
thickness : float, optional
thickness, by default 1.0
"""
_HyperElas.__init__(self, dim, thickness)
self.C0 = C0
self.C1 = C1
self.C2 = C2
self.C3 = C3
self.C4 = C4
self.C5 = C5
self.C6 = C6
self.C7 = C7
self.K = K
self.Mu1 = Mu1
self.Mu2 = Mu2
self.T1 = T1
self.T2 = T2
self.__ks = ks
[docs]
def Compute_W(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
C0 = self.C0
C1 = self.C1
C2 = self.C2
C3 = self.C3
C4 = self.C4
C5 = self.C5
C6 = self.C6
C7 = self.C7
K = self.K
Mu1 = self.Mu1
Mu2 = self.Mu2
T1 = self.T1
T2 = self.T2
ks = self.__ks
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
I4 = hyperElasticState.Compute_I4(T1)
I6 = hyperElasticState.Compute_I6(T2)
I8 = hyperElasticState.Compute_I8(T1, T2)
W = (
C0 * (np.exp(C1 * (I1 / I3 ** (1 / 3) - 3)) - 1)
+ C2 * (np.exp(C3 * (I4 - 1) ** 2) - 1) / (1 + np.exp(-ks * (I4 - 1)))
+ C4 * (np.exp(C5 * (I6 - 1) ** 2) - 1) / (1 + np.exp(-ks * (I6 - 1)))
+ C6 * (np.exp(C7 * I8**2) - 1)
+ K * (I3 - 2 * np.log(np.sqrt(I3)) - 1) / 4
+ Mu1 * (I1 / I3 ** (1 / 3) - 3)
+ Mu2 * (I2 / I3 ** (2 / 3) - 3)
)
return W
[docs]
def Compute_dWde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
C0 = self.C0
C1 = self.C1
C2 = self.C2
C3 = self.C3
C4 = self.C4
C5 = self.C5
C6 = self.C6
C7 = self.C7
K = self.K
Mu1 = self.Mu1
Mu2 = self.Mu2
T1 = self.T1
T2 = self.T2
ks = self.__ks
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
I4 = hyperElasticState.Compute_I4(T1)
I6 = hyperElasticState.Compute_I6(T2)
I8 = hyperElasticState.Compute_I8(T1, T2)
dI1dC = hyperElasticState.Compute_dI1dC()
dI2dC = hyperElasticState.Compute_dI2dC()
dI3dC = hyperElasticState.Compute_dI3dC()
dI4dC = hyperElasticState.Compute_dI4dC(T1)
dI6dC = hyperElasticState.Compute_dI6dC(T2)
dI8dC = hyperElasticState.Compute_dI8dC(T1, T2)
# see: examples/HyperElastic/HyperElasticLaws.py
# TODO Optimize
# fmt: off
dWdI1 = C0*C1*np.exp(C1*(I1/I3**(1/3) - 3))/I3**(1/3) + Mu1/I3**(1/3)
dWdI2 = Mu2/I3**(2/3)
dWdI3 = -C0*C1*I1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - I1*Mu1/(3*I3**(4/3)) - 2*I2*Mu2/(3*I3**(5/3)) + K*(1 - 1/I3)/4
dWdI4 = C2*C3*(2*I4 - 2)*np.exp(C3*(I4 - 1)**2)/(1 + np.exp(-ks*(I4 - 1))) + C2*ks*(np.exp(C3*(I4 - 1)**2) - 1)*np.exp(-ks*(I4 - 1))/(1 + np.exp(-ks*(I4 - 1)))**2
dWdI6 = C4*C5*(2*I6 - 2)*np.exp(C5*(I6 - 1)**2)/(1 + np.exp(-ks*(I6 - 1))) + C4*ks*(np.exp(C5*(I6 - 1)**2) - 1)*np.exp(-ks*(I6 - 1))/(1 + np.exp(-ks*(I6 - 1)))**2
dWdI8 = 2*C6*C7*I8*np.exp(C7*I8**2)
# fmt: on
dW = 2 * (
dWdI1 * dI1dC
+ dWdI2 * dI2dC
+ dWdI3 * dI3dC
+ dWdI4 * dI4dC
+ dWdI6 * dI6dC
+ dWdI8 * dI8dC
)
return dW
[docs]
def Compute_d2Wde(self, mesh, u, matrixType=MatrixType.rigi) -> FeArray:
C0 = self.C0
C1 = self.C1
C2 = self.C2
C3 = self.C3
C4 = self.C4
C5 = self.C5
C6 = self.C6
C7 = self.C7
K = self.K
Mu1 = self.Mu1
Mu2 = self.Mu2
T1 = self.T1
T2 = self.T2
ks = self.__ks
hyperElasticState = HyperElasticState(mesh, u, matrixType)
I1 = hyperElasticState.Compute_I1()
I2 = hyperElasticState.Compute_I2()
I3 = hyperElasticState.Compute_I3()
I4 = hyperElasticState.Compute_I4(T1)
I6 = hyperElasticState.Compute_I6(T2)
I8 = hyperElasticState.Compute_I8(T1, T2)
dI1dC = hyperElasticState.Compute_dI1dC()
dI2dC = hyperElasticState.Compute_dI2dC()
dI3dC = hyperElasticState.Compute_dI3dC()
dI4dC = hyperElasticState.Compute_dI4dC(T1)
dI6dC = hyperElasticState.Compute_dI6dC(T2)
dI8dC = hyperElasticState.Compute_dI8dC(T1, T2)
d2I1dC = hyperElasticState.Compute_d2I1dC()
d2I2dC = hyperElasticState.Compute_d2I2dC()
d2I3dC = hyperElasticState.Compute_d2I3dC()
d2I4dC = hyperElasticState.Compute_d2I4dC()
d2I6dC = hyperElasticState.Compute_d2I6dC()
d2I8dC = hyperElasticState.Compute_d2I8dC()
# see: examples/HyperElastic/HyperElasticLaws.py
# TODO Optimize
# fmt: off
dWdI1 = C0*C1*np.exp(C1*(I1/I3**(1/3) - 3))/I3**(1/3) + Mu1/I3**(1/3)
d2WdI1dI1 = C0*C1**2*np.exp(C1*(I1/I3**(1/3) - 3))/I3**(2/3)
d2WdI1dI3 = -C0*C1**2*I1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(5/3)) - C0*C1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - Mu1/(3*I3**(4/3))
dWdI2 = Mu2/I3**(2/3)
d2WdI2dI3 = -2*Mu2/(3*I3**(5/3))
dWdI3 = -C0*C1*I1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - I1*Mu1/(3*I3**(4/3)) - 2*I2*Mu2/(3*I3**(5/3)) + K*(1 - 1/I3)/4
d2WdI3dI1 = -C0*C1**2*I1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(5/3)) - C0*C1*np.exp(C1*(I1/I3**(1/3) - 3))/(3*I3**(4/3)) - Mu1/(3*I3**(4/3))
d2WdI3dI2 = -2*Mu2/(3*I3**(5/3))
d2WdI3dI3 = C0*C1**2*I1**2*np.exp(C1*(I1/I3**(1/3) - 3))/(9*I3**(8/3)) + 4*C0*C1*I1*np.exp(C1*(I1/I3**(1/3) - 3))/(9*I3**(7/3)) + 4*I1*Mu1/(9*I3**(7/3)) + 10*I2*Mu2/(9*I3**(8/3)) + K/(4*I3**2)
dWdI4 = C2*C3*(2*I4 - 2)*np.exp(C3*(I4 - 1)**2)/(1 + np.exp(-ks*(I4 - 1))) + C2*ks*(np.exp(C3*(I4 - 1)**2) - 1)*np.exp(-ks*(I4 - 1))/(1 + np.exp(-ks*(I4 - 1)))**2
d2WdI4dI4 = C2*C3**2*(2*I4 - 2)**2*np.exp(C3*(I4 - 1)**2)/(1 + np.exp(-ks*(I4 - 1))) + 2*C2*C3*ks*(2*I4 - 2)*np.exp(C3*(I4 - 1)**2)*np.exp(-ks*(I4 - 1))/(1 + np.exp(-ks*(I4 - 1)))**2 + 2*C2*C3*np.exp(C3*(I4 - 1)**2)/(1 + np.exp(-ks*(I4 - 1))) - C2*ks**2*(np.exp(C3*(I4 - 1)**2) - 1)*np.exp(-ks*(I4 - 1))/(1 + np.exp(-ks*(I4 - 1)))**2 + 2*C2*ks**2*(np.exp(C3*(I4 - 1)**2) - 1)*np.exp(-2*ks*(I4 - 1))/(1 + np.exp(-ks*(I4 - 1)))**3
dWdI6 = C4*C5*(2*I6 - 2)*np.exp(C5*(I6 - 1)**2)/(1 + np.exp(-ks*(I6 - 1))) + C4*ks*(np.exp(C5*(I6 - 1)**2) - 1)*np.exp(-ks*(I6 - 1))/(1 + np.exp(-ks*(I6 - 1)))**2
d2WdI6dI6 = C4*C5**2*(2*I6 - 2)**2*np.exp(C5*(I6 - 1)**2)/(1 + np.exp(-ks*(I6 - 1))) + 2*C4*C5*ks*(2*I6 - 2)*np.exp(C5*(I6 - 1)**2)*np.exp(-ks*(I6 - 1))/(1 + np.exp(-ks*(I6 - 1)))**2 + 2*C4*C5*np.exp(C5*(I6 - 1)**2)/(1 + np.exp(-ks*(I6 - 1))) - C4*ks**2*(np.exp(C5*(I6 - 1)**2) - 1)*np.exp(-ks*(I6 - 1))/(1 + np.exp(-ks*(I6 - 1)))**2 + 2*C4*ks**2*(np.exp(C5*(I6 - 1)**2) - 1)*np.exp(-2*ks*(I6 - 1))/(1 + np.exp(-ks*(I6 - 1)))**3
dWdI8 = 2*C6*C7*I8*np.exp(C7*I8**2)
d2WdI8dI8 = 4*C6*C7**2*I8**2*np.exp(C7*I8**2) + 2*C6*C7*np.exp(C7*I8**2)
# fmt: on
d2W = 4 * (
dWdI1 * d2I1dC
+ dWdI2 * d2I2dC
+ dWdI3 * d2I3dC
+ dWdI4 * d2I4dC
+ dWdI6 * d2I6dC
+ dWdI8 * d2I8dC
) + 4 * (
d2WdI1dI1 * TensorProd(dI1dC, dI1dC)
+ d2WdI1dI3 * TensorProd(dI1dC, dI3dC)
+ d2WdI2dI3 * TensorProd(dI2dC, dI3dC)
+ d2WdI3dI1 * TensorProd(dI3dC, dI1dC)
+ d2WdI3dI2 * TensorProd(dI3dC, dI2dC)
+ d2WdI3dI3 * TensorProd(dI3dC, dI3dC)
+ d2WdI4dI4 * TensorProd(dI4dC, dI4dC)
+ d2WdI6dI6 * TensorProd(dI6dC, dI6dC)
+ d2WdI8dI8 * TensorProd(dI8dC, dI8dC)
)
return d2W