Models

The EasyFEA.Models module provides essential tools for creating and managing models. These models are used to build _Simu instances and mainly contain material parameters.

In the simulation workflow, Models is the third step: once the mesh exists, a model encapsulates the physics and material constants (Young’s modulus, thermal conductivity, fracture toughness, …) before being passed to the simulation constructor.

With this module, you can construct:

• Linear elastic materials, such as Isotropic, TransverselyIsotropic, Orthotropic, and Anisotropic, in Models.Elastic.

• Nonlinear hyperelastic materials, such as NeoHookean, CiarletGeymonat, MooneyRivlin, SaintVenantKirchhoff, and HolzapfelOgden, in Models.HyperElastic.

• Elastic beams with Isotropic, BeamStructure, in Models.Beam.

• Phase-field materials with PhaseField.

• Thermal materials with Thermal.

• Weak forms with WeakForms.

See also

• Create a model

Models API

Module implementing constitutive laws used in simulations.

class EasyFEA.Models.ModelType(value)

Bases: str, Enum

Model types.

classmethod types()

beam = 'beam'

damage = 'damage'

elastic = 'elastic'

hyperelastic = 'hyperelastic'

thermal = 'thermal'

weakForm = 'weakForm'

class EasyFEA.Models.PhaseField(

material: _Elastic,

split: SplitType,

regularization: ReguType,

Gc: float | ndarray[tuple[Any, ...], dtype[floating]],

l0: float,

solver=SolverType.History,

A: ndarray[tuple[Any, ...], dtype[floating]] | None = None,

)

Bases: _IModel

PhaseField class.

class ReguType(value)

Bases: str, Enum

Regularization Models.

AT1 = 'AT1'

AT2 = 'AT2'

class SolverType(value)

Bases: str, Enum

Solver used to manage crack irreversibility.

BoundConstrain = 'BoundConstrain'

History = 'History'

HistoryDamage = 'HistoryDamage'

class SplitType(value)

Bases: str, Enum

Split Models.

Amor = 'Amor'

AnisotStrain = 'AnisotStrain'

AnisotStrain_MP = 'AnisotStrain_MP'

AnisotStrain_NoCross = 'AnisotStrain_NoCross'

AnisotStrain_PM = 'AnisotStrain_PM'

AnisotStress = 'AnisotStress'

AnisotStress_MP = 'AnisotStress_MP'

AnisotStress_NoCross = 'AnisotStress_NoCross'

AnisotStress_PM = 'AnisotStress_PM'

Bourdin = 'Bourdin'

He = 'He'

Miehe = 'Miehe'

Stress = 'Stress'

Zhang = 'Zhang'

static Get_regularizations() → list[ReguType]

Returns regularizations available

static Get_solvers() → list[SolverType]

Returns available solvers used to manage crack irreversibility

static Get_splits() → list[SplitType]

Returns splits available

Calc_C(

Epsilon_e_pg: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

verif=False,

) → tuple[FeArray | ndarray[tuple[Any, ...], dtype[Any]], FeArray | ndarray[tuple[Any, ...], dtype[Any]]]

Computes the splited stifness matrices for the given strain field.

Parameters:

Epsilon_e_pg (FeArray.FeArrayALike) – strains field (e, p, D)

Returns:

cP_e_pg, cM_e_pg: positive and negative stifness matrices (e, p, D, D)

Return type:

FeArray

Calc_Sigma_e_pg(

Epsilon_e_pg: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

) → tuple[FeArray | ndarray[tuple[Any, ...], dtype[Any]], FeArray | ndarray[tuple[Any, ...], dtype[Any]]]

Computes the Stress field using the strains and the split such that:

SigmaP_e_pg = cP_e_pg * Epsilon_e_pg

SigmaM_e_pg = cM_e_pg * Epsilon_e_pg

Parameters:

Epsilon_e_pg (FeArray.FeArrayALike) – strains field (e, p, D)

Returns:

SigmaP_e_pg, SigmaM_e_pg: positive and negative stress fields (e, p, D)

Return type:

FeArray

Calc_psi_e_pg(

Epsilon_e_pg: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

) → tuple[FeArray | ndarray[tuple[Any, ...], dtype[Any]], FeArray | ndarray[tuple[Any, ...], dtype[Any]]]

Computes the elastic energy densities.

psiP_e_pg = 1/2 SigmaP_e_pg * Epsilon_e_pg

psiM_e_pg = 1/2 SigmaM_e_pg * Epsilon_e_pg

Such as :

SigmaP_e_pg = cP_e_pg * Epsilon_e_pg

SigmaM_e_pg = cM_e_pg * Epsilon_e_pg

Get_f_e_pg(

PsiP_e_pg: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Returns source therm

Get_g_e_pg(

d_n: ndarray[tuple[Any, ...], dtype[floating]],

groupElem: _GroupElem,

matrixType: MatrixType,

k_res=1e-12,

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Returns degradation function

Get_r_e_pg(

PsiP_e_pg: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Returns reaction therm

_Eigen_values_vectors_projectors(

vector_e_pg: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

verif=False,

) → tuple[FeArray | ndarray[tuple[Any, ...], dtype[Any]], list[FeArray | ndarray[tuple[Any, ...], dtype[Any]]], list[FeArray | ndarray[tuple[Any, ...], dtype[Any]]]]

Computes the eigen values and eigen projectors of a second-order tensor (as a vector).

A

matrix characterizing the weak anisotropy in the crack surface density function

Gc: float

critical energy release rate

property c_w

scaling parameter for accurate dissipation of crack energy

property dim: int

model dimension

property isHeterogeneous: bool

indicates whether the model has heterogeneous parameters

property k: float | ndarray[tuple[Any, ...], dtype[floating]]

get diffusion therm

l0: float

half crack width

property material: _Elastic

elastic material

property modelType: ModelType

model type

regularization: ReguType

crack regularization model

solver: SolverType

solver used to manage crack irreversibility

split: SplitType

split used to decompose the elastic energy density

property thickness: float

thickness used in the model

class EasyFEA.Models.Thermal(k: float, c=0.0, thickness: float = 1.0)

Bases: _IModel

Thermal class.

c

mass heat capacity [J K^-1 kg^-1]

dim: int

property isHeterogeneous: bool

indicates whether the model has heterogeneous parameters

k

thermal conductivity [W m^-1]

property modelType: ModelType

model type

thickness: float

class EasyFEA.Models.WeakForms(

field: Field,

computeK: BiLinearForm,

computeC: BiLinearForm | None = None,

computeM: BiLinearForm | None = None,

computeF: LinearForm | None = None,

thickness: float = 1.0,

)

Bases: _IModel

Class responsible for computing the finite element matrices used in the system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\).

property computeC: BiLinearForm | None

Function used to build stiffness matrix \(\Crm\).

property computeF: LinearForm | None

Function used to build force vector \(\Frm\).

property computeK: BiLinearForm | None

Function used to build stiffness matrix \(\Krm\).

property computeM: BiLinearForm | None

Function used to build stiffness matrix \(\Mrm\).

property dim: int

model dimension

property field: Field

Finite element field.

property modelType: ModelType

model type

thickness: float

class EasyFEA.Models._IModel

Bases: Observable, Updatable, ABC

Model interface.

Need_Update(value=True)

Indicates whether the object needs to be updated.

abstract property dim: int

model dimension

property isHeterogeneous: bool

indicates whether the model has heterogeneous parameters

abstract property modelType: ModelType

model type

abstract property thickness: float

thickness used in the model

EasyFEA.Models.Apply_Pmat(

P: ndarray[tuple[Any, ...], dtype[floating]],

M: ndarray[tuple[Any, ...], dtype[floating]],

toGlobal=True,

) → ndarray[tuple[Any, ...], dtype[floating]]

Performs a basis transformation from the material’s coordinate system to the (x,y,z) coordinate system to orient the material in space.

Caution: P must be in Kelvin mandel notation

Parameters:

• P (_types.FloatArray) – P in mandel notation obtained with Get_Pmat

• M (_types.FloatArray) – 3x3 or 6x6 matrix

• toGlobal (bool, optional) –

  sets wheter you want to get matrix in global or material coordinates, by default True

  if toGlobal:

  Matrix_global = P * C_material * P’

  else:

  Matrix_material = P’ * Matrix_global * P

Returns:

new matrix

Return type:

_types.FloatArray

EasyFEA.Models.Get_Pmat(

axis_1: ndarray[tuple[Any, ...], dtype[floating]],

axis_2: ndarray[tuple[Any, ...], dtype[floating]],

useMandel=True,

)

Constructs Pmat to pass from the material coordinates (x,y,z) to the global coordinate (X,Y,Z) such that:

if useMandel:

return [Pm]

else:

return [Ps], [Pe]

In Kelvin Mandel notation:

Sig & Eps en [11, 22, 33, sqrt(2)*23, sqrt(2)*13, sqrt(2)*12]

[C_global] = [Pm] * [C_material] * [Pm]^T & [C_material] = [Pm]^T * [C_global] * [Pm]

[S_global] = [Pm] * [S_material] * [Pm]^T & [S_material] = [Pm]^T * [S_global] * [Pm]

Sig_global = [Pm] * Sig_material & Sig_material = [Pm]^T * Sig_global

Eps_global = [Pm] * Eps_material & Eps_material = [Pm]^T * Eps_global

In Voigt’s notation:

Sig [S11, S22, S33, S23, S13, S12]

Eps [E11, E22, E33, 2*E23, 2*E13, 2*E12]

[C_global] = [Ps] * [C_material] * [Ps]^T & [C_material] = [Pe]^T * [C_global] * [Pe]

S_global = [Pe] * [S_material] * [Pe]^T & [S_material] = [Ps]^T * S_global * [Ps]

Sig_global = [Ps] * Sig_material & Sig_material = [Pe]^T * Sig_global

Eps_global = [Pe] * Eps_material & Eps_material = [Ps]^T * Eps_global

P matrices are orhogonal such that: inv([P]) = [P]^T

Here we use “Chevalier 1988 : Comportements élastique et viscoélastique des composites”

EasyFEA.Models.Heterogeneous_Array(

array: ndarray[tuple[Any, ...], dtype[floating]],

)

Builds a heterogeneous array.

EasyFEA.Models.KelvinMandel_Matrix(

dim: int,

M: ndarray[tuple[Any, ...], dtype[floating]],

) → ndarray[tuple[Any, ...], dtype[floating]]

Apply Kelvin Mandel coefficient to constitutive laws.

In 2D:

[1, 1, r2]

[1, 1, r2]

[r2, r2, 2]]

In 3D:

[1,1,1,r2,r2,r2]

[1,1,1,r2,r2,r2]

[1,1,1,r2,r2,r2]

[r2,r2,r2,2,2,2]

[r2,r2,r2,2,2,2]

[r2,r2,r2,2,2,2]]

EasyFEA.Models.Project_Kelvin(

A: ndarray[tuple[Any, ...], dtype[floating]],

orderA: int | None = None,

) → ndarray[tuple[Any, ...], dtype[floating]]

Projects the tensor A in Kelvin Mandel notation.

Parameters:

• A (_types.FloatArray) – tensor A (2 or 4 order tensor)

• orderA (int, optional) – tensor order, by default None

Returns:

Projected tensor

Return type:

_types.FloatArray

EasyFEA.Models.Project_matrix_to_vector(

matrix: ndarray[tuple[Any, ...], dtype[floating]],

coef=np.float64(1.4142135623730951),

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

EasyFEA.Models.Project_vector_to_matrix(

vector: ndarray[tuple[Any, ...], dtype[floating]],

coef=np.float64(1.4142135623730951),

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

EasyFEA.Models.Reshape_variable(

variable: int | float | ndarray[tuple[Any, ...], dtype[Any]],

Ne: int,

nPg: int,

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Resizes variable to (Ne, nPg, …) shape.

EasyFEA.Models.Result_in_Strain_or_Stress_field(

field_e_pg: FeArray,

result: str,

coef=np.float64(1.4142135623730951),

) → ndarray[tuple[Any, ...], dtype[floating]]

Extracts a specific result from a 2D or 3D strain or stress field.

Parameters:

• field_e_pg (_types.FloatArray) – Strain or stress field in each element and gauss points.

• result (str) –

  Desired result/value to extract:

  2D: [xx, yy, xy, vm, Strain, Stress]

  3D: [xx, yy, zz, yz, xz, xy, vm, Strain, Stress]

• coef (float, optional) – Coefficient used to scale cross components in the field (e.g., xy/coef in dim=2, or yz/coef, xz/coef, xy/coef if dim=3).

Returns:

The extracted field corresponding to the specified result.

Return type:

_types.FloatArray

Elastic module.

class EasyFEA.Models.Elastic.Anisotropic(

dim: int,

C: ndarray[tuple[Any, ...], dtype[floating]],

useVoigtNotation: bool,

axis1: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (1, 0, 0),

axis2: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (0, 1, 0),

thickness=1.0,

)

Bases: _Elastic

Anisotropic Linearized Elastic material.

Set_C(

C: ndarray[tuple[Any, ...], dtype[floating]],

useVoigtNotation=True,

update_S=True,

)

Updates the constitutives laws by updating the C stiffness and S compliance matrices in Kelvin Mandel notation.

Parameters:

• C (_types.FloatArray) – Stifness matrix (Lamé’s law)

• useVoigtNotation (bool, optional) – uses Kevin Mandel’s notation, by default True

• update_S (bool, optional) – updates the compliance matrix (Hooke’s law), by default True

Walpole_Decomposition() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Walpole’s decomposition in Kelvin Mandel notation such that:

C = sum(ci * Ei).

returns ci, Ei

_Behavior(

C: ndarray[tuple[Any, ...], dtype[floating]],

useVoigtNotation: bool,

) → ndarray[tuple[Any, ...], dtype[floating]]

_Update() → None

Updates the constitutives laws by updating the C stiffness and S compliance matrices. in Kelvin Mandel notation

property axis1: ndarray[tuple[Any, ...], dtype[floating]]

axis1 vector

property axis2: ndarray[tuple[Any, ...], dtype[floating]]

axis2 vector

class EasyFEA.Models.Elastic.Isotropic(dim: int, E=210000.0, v=0.3, planeStress=True, thickness=1.0)

Bases: _Elastic

Isotropic Linearized Elastic material.

Walpole_Decomposition() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Walpole’s decomposition in Kelvin Mandel notation such that:

C = sum(ci * Ei).

returns ci, Ei

_Behavior(dim: int | None = None)

Updates the constitutives laws by updating the C stiffness and S compliance matrices in Kelvin Mandel notation.

In 2D:

C -> C : Epsilon = Sigma [Sxx Syy sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy sqrt(2)*Exy]

In 3D:

C -> C : Epsilon = Sigma [Sxx Syy Szz sqrt(2)*Syz sqrt(2)*Sxz sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy Ezz sqrt(2)*Eyz sqrt(2)*Exz sqrt(2)*Exy]

_Update() → None

Updates the constitutives laws by updating the C stiffness and S compliance matrices. in Kelvin Mandel notation

get_bulk()

Bulk modulus

get_lambda()

get_mu()

Shear coefficient

E: float

Young’s modulus

v: float

Poisson’s ratio (-1<v<0.5)

class EasyFEA.Models.Elastic.Orthotropic(

dim: int,

E1: float,

E2: float,

E3: float,

G23: float,

G13: float,

G12: float,

v23: float,

v13: float,

v12: float,

axis_1: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (1, 0, 0),

axis_2: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (0, 1, 0),

planeStress: bool = True,

thickness: float = 1.0,

)

Bases: _Elastic

Orthotropic Linearized Elastic material.

Walpole_Decomposition() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Walpole’s decomposition in Kelvin Mandel notation such that:

C = sum(ci * Ei).

returns ci, Ei

_Behavior(dim: int | None = None)

Updates the constitutives laws by updating the C stiffness and S compliance matrices in Kelvin Mandel notation.

In 2D:

C -> C : Epsilon = Sigma [Sxx Syy sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy sqrt(2)*Exy]

In 3D:

C -> C : Epsilon = Sigma [Sxx Syy Szz sqrt(2)*Syz sqrt(2)*Sxz sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy Ezz sqrt(2)*Eyz sqrt(2)*Exz sqrt(2)*Exy]

_Update() → None

Updates the constitutives laws by updating the C stiffness and S compliance matrices. in Kelvin Mandel notation

E1: float

Young’s modulus along axis_1.

E2: float

Young’s modulus along axis_2.

E3: float

Young’s modulus along axis_3.

G12: float

Shear modulus in the 1-2 plane.

G13: float

Shear modulus in the 1-3 plane.

G23: float

Shear modulus in the 2-3 plane.

property _c11: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c12: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c13: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c22: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c23: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c33: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c44: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c55: float | ndarray[tuple[Any, ...], dtype[floating]]

property _c66: float | ndarray[tuple[Any, ...], dtype[floating]]

property axis_1: ndarray[tuple[Any, ...], dtype[floating]]

Axis 1

property axis_2: ndarray[tuple[Any, ...], dtype[floating]]

Axis 2

v12: float

Poisson’s ratio for transverse strain along the axis_2 when stressed along the axis_1.

v13: float

Poisson’s ratio for transverse strain along the axis_3 when stressed along the axis_1.

v23: float

Poisson’s ratio for transverse strain along the axis_3 when stressed along the axis_2.

class EasyFEA.Models.Elastic.TransverselyIsotropic(

dim: int,

El: float,

Et: float,

Gl: float,

vl: float,

vt: float,

axis_l: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (1, 0, 0),

axis_t: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (0, 1, 0),

planeStress: bool = True,

thickness: float = 1.0,

)

Bases: _Elastic

Transversely Isotropic Linearized Elastic material.

Walpole_Decomposition() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Walpole’s decomposition in Kelvin Mandel notation such that:

C = sum(ci * Ei).

returns ci, Ei

_Behavior(dim: int | None = None)

Updates the constitutives laws by updating the C stiffness and S compliance matrices in Kelvin Mandel notation.

In 2D:

C -> C : Epsilon = Sigma [Sxx Syy sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy sqrt(2)*Exy]

In 3D:

C -> C : Epsilon = Sigma [Sxx Syy Szz sqrt(2)*Syz sqrt(2)*Sxz sqrt(2)*Sxy]

S -> S : Sigma = Epsilon [Exx Eyy Ezz sqrt(2)*Eyz sqrt(2)*Exz sqrt(2)*Exy]

_Update() → None

Updates the constitutives laws by updating the C stiffness and S compliance matrices. in Kelvin Mandel notation

El: float

Longitudinal Young’s modulus.

Et: float

Transverse Young’s modulus.

Gl: float

Longitudinal shear modulus.

property Gt: float | ndarray[tuple[Any, ...], dtype[floating]]

Transverse shear modulus.

property axis_l: ndarray[tuple[Any, ...], dtype[floating]]

Longitudinal axis

property axis_t: ndarray[tuple[Any, ...], dtype[floating]]

Transversal axis

property kt: float | ndarray[tuple[Any, ...], dtype[floating]]

vl: float

Longitudinal Poisson’s ratio (-1<vl<0.5).

vt: float

Transverse Poisson ratio (-1<vt<1)

class EasyFEA.Models.Elastic._Elastic(dim: int, thickness: float, planeStress: bool)

Bases: _IModel, ABC

Linearized Elasticity material.

ElasIsot, ElasIsotTrans and ElasAnisot inherit from _Elas class.

static Available_Laws()

Get_sqrt_C_S() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Returns the Matrix square root of C and S.

abstractmethod Walpole_Decomposition() → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Walpole’s decomposition in Kelvin Mandel notation such that:

C = sum(ci * Ei).

returns ci, Ei

_Apply_basis_transformation(

dim: int,

material_cM: ndarray[tuple[Any, ...], dtype[floating]],

material_sM: ndarray[tuple[Any, ...], dtype[floating]],

axis_1: ndarray[tuple[Any, ...], dtype[floating]],

axis_2: ndarray[tuple[Any, ...], dtype[floating]],

) → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Performs a basis transformation from the material’s (1,2,3) coordinate system to the (x,y,z) coordinate system to orient the material in space.

Parameters:

• dim (int) – dimension

• material_cM (_types.FloatArray) – stiffness matrix

• material_sM (_types.FloatArray) – compliance matrix

• axis_1 (_types.FloatArray) – Axis 1

• axis_2 (_types.FloatArray) – Axis 2

Returns:

global_cM, global_sM

Return type:

tuple[_types.FloatArray, _types.FloatArray]

abstractmethod _Update() → None

Updates the constitutives laws by updating the C stiffness and S compliance matrices. in Kelvin Mandel notation

property C: ndarray[tuple[Any, ...], dtype[floating]]

Stifness matrix in Kelvin Mandel notation such that:

In 2D: C -> C: Epsilon = Sigma [Sxx, Syy, sqrt(2)*Sxy]

In 3D: C -> C: Epsilon = Sigma [Sxx, Syy, Szz, sqrt(2)*Syz, sqrt(2)*Sxz, sqrt(2)*Sxy].

(Lame’s law)

property S: ndarray[tuple[Any, ...], dtype[floating]]

Compliance matrix in Kelvin Mandel notation such that:

In 2D: S -> S : Sigma = Epsilon [Exx, Eyy, sqrt(2)*Exy]

In 3D: S -> S: Sigma = Epsilon [Exx, Eyy, Ezz, sqrt(2)*Eyz, sqrt(2)*Exz, sqrt(2)*Exy].

(Hooke’s law)

property coef: float

kelvin mandel coef -> sqrt(2)

dim: int

property isHeterogeneous: bool

indicates whether the model has heterogeneous parameters

property modelType: ModelType

model type

planeStress: bool

the model uses plane stress simplification

property simplification: str

simplification used for the model

thickness: float

Hyper elastic module.

class EasyFEA.Models.HyperElastic.CiarletGeymonat(

dim: int,

K1: float | ndarray[tuple[Any, ...], dtype[floating]],

K2: float | ndarray[tuple[Any, ...], dtype[floating]],

K: float | ndarray[tuple[Any, ...], dtype[floating]] = 0.0,

thickness=1.0,

)

Bases: _HyperElastic

Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

_Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

K: float

Bulk modulus

K1: float

Kappa1

K2: float

Kappa2

class EasyFEA.Models.HyperElastic.HolzapfelOgden(

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: ndarray[tuple[Any, ...], dtype[floating]],

T2: ndarray[tuple[Any, ...], dtype[floating]],

ks: float = 100,

thickness=1.0,

)

Bases: _HyperElastic

Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

_Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

C0: float

C1: float

C2: float

C3: float

C4: float

C5: float

C6: float

C7: float

K: float

Bulk modulus

Mu1: float

Mu2: float

T1

direction(s) 1, used for the invariants I4 and I8

T2

direction(s) 2, used for the invariants I6 and I8

class EasyFEA.Models.HyperElastic.HyperElasticState(

groupElem: _GroupElem,

displacement: ndarray[tuple[Any, ...], dtype[floating]],

matrixType: MatrixType,

)

Bases: object

Hyperelastic state.

static _CheckFormat(

groupElem: _GroupElem,

u: ndarray[tuple[Any, ...], dtype[floating]],

) → None

Compute_C() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Right Cauchy-Green tensor C(u) = F(u)ᵀ · F(u) at the Gauss points, padded to 3×3.

Returns:

Layout depends on dim (unused off-diagonal terms are zero; missing diagonal terms stay at 1, so the result is always 3×3):

dim == 1: ` [ cxx   0    0  ] [  0    1    0  ] [  0    0    1  ] `

dim == 2: ` [ cxx  cxy   0  ] [ cyx  cyy   0  ] [  0    0    1  ] `

dim == 3: ` [ cxx  cxy  cxz ] [ cyx  cyy  cyz ] [ czx  czy  czz ] `

Shape: (Ne, nPg, 3, 3).

Return type:

FeArray

Compute_De() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Kinematic operator such that Ė_vec = De(u) · flat(∇v) where Ė_vec is the Green-Lagrange strain rate in Kelvin-Mandel vector form.

It is __Build_De with G = F = I + ∇u the deformation gradient (so De · flat(∇v) = sym(Fᵀ∇v) = Ė).

With c = 1/√2 (Kelvin-Mandel scaling on shear rows):

2D — 3 strain components (Ėxx, Ėyy, √2·Ėxy), columns index flat(∇v) = [dxux, dyux, dxuy, dyuy]: ` [ 1+dxux       0         dxuy         0      ] [    0       dyux         0        1+dyuy    ] [ c·dyux  c·(1+dxux)  c·(1+dyuy)   c·dxuy   ] `

3D — 6 strain components (Ėxx, Ėyy, Ėzz, √2·Ėyz, √2·Ėxz, √2·Ėxy), columns index flat(∇v) = [dxux, dyux, dzux, dxuy, dyuy, dzuy, dxuz, dyuz, dzuz]: ` [ 1+dxux      0          0         dxuy        0          0          dxuz        0          0      ] [    0      dyux         0           0       1+dyuy       0            0        dyuz        0      ] [    0        0        dzux          0          0        dzuy          0          0       1+dzuz   ] [    0    c·dzux     c·dyux          0       c·dzuy   c·(1+dyuy)       0     c·(1+dzuz)  c·dyuz   ] [ c·dzux     0      c·(1+dxux)   c·dzuy        0       c·dxuy     c·(1+dzuz)     0       c·dxuz   ] [ c·dyux c·(1+dxux)    0       c·(1+dyuy)   c·dxuy       0          c·dyuz     c·dxuz        0    ] `

Shape: (Ne, nPg, 3, 4) in 2D, (Ne, nPg, 6, 9) in 3D.

Compute_Deta(

velocity: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Configuration derivative of the Green-Lagrange strain rate: the operator Deta such that ∂Ė_vec = Deta · flat(∇δu) at fixed velocity.

Ė_vec = De(u) · flat(∇v) = sym(Fᵀ∇v) is bilinear in (∇u, ∇v), so ∂Ė/∂(∇u) = sym(∇vᵀ∇δu) is __Build_De with G = ∇v the velocity gradient (the same builder as Compute_De(), with ∇v in place of F). It feeds the material-like piece η · Bᵀ · (Deta · grad) of the viscous configuration tangent (Kgeo) returned by Operators.NonLinear.KelvinVoigtDamping().

Parameters:

• velocity (_types.FloatArray) – velocity field, same (xi,yi,zi,...) layout as the displacement.

• ``(Ne (Returns) –

• nPg

• 3

• 2D (4)`` in)

• ``(Ne –

• nPg

• 6

:param 9)`` in 3D — same layout as Compute_De().:

Compute_Edot_vec(

velocity: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Green–Lagrange strain rate Ė in Kelvin–Mandel vector form.

Uses the identity Ė = sym(Fᵀ · ∇v) = De(u) · flat(∇v) — the same kinematic operator Compute_De() that maps flat(∇u̇) to the small-strain rate also maps flat(∇v) to Ė_vec evaluated at the current u.

Parameters:

• velocity (_types.FloatArray) – velocity field, same (xi,yi,zi,...) layout as the displacement.

• ``(Ne (Returns) –

• nPg

• 2D (nstrain)`` — nstrain = 3 in)

• 3D. (6 in)

Compute_Epsilon() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Linearized strain ε = 1/2 (∇uᵀ + ∇u) in Kelvin-Mandel vector form at the Gauss points.

Returns:

Layout depends on dim (off-diagonal components carry the Kelvin-Mandel √2 factor):

dim == 2: ` [ εxx,  εyy,  √2·εxy ] ` Shape: (Ne, nPg, 3).

dim == 3: ` [ εxx,  εyy,  εzz,  √2·εyz,  √2·εxz,  √2·εxy ] ` Shape: (Ne, nPg, 6).

Return type:

FeArray

Compute_F() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Deformation gradient F(u) = I + ∇u at the Gauss points, padded to 3×3.

Returns:

Layout depends on dim (unused off-diagonal terms are zero; missing diagonal terms stay at 1, so the result is always 3×3):

dim == 1: ` [ 1+dxux    0       0    ] [   0       1       0    ] [   0       0       1    ] `

dim == 2: ` [ 1+dxux   dyux     0    ] [  dxuy   1+dyuy    0    ] [   0       0       1    ] `

dim == 3: ` [ 1+dxux   dyux    dzux  ] [  dxuy   1+dyuy   dzuy  ] [  dxuz    dyuz   1+dzuz ] `

Shape: (Ne, nPg, 3, 3).

Return type:

FeArray

Compute_GreenLagrange() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Green-Lagrange strain E = 1/2 (C - I) at the Gauss points, padded to 3×3.

Returns:

Layout depends on dim (unused components are zeroed so the result is always 3×3):

dim == 1: ` [ Exx   0    0  ] [  0    0    0  ] [  0    0    0  ] `

dim == 2: ` [ Exx  Exy   0  ] [ Eyx  Eyy   0  ] [  0    0    0  ] `

dim == 3: ` [ Exx  Exy  Exz ] [ Eyx  Eyy  Eyz ] [ Ezx  Ezy  Ezz ] `

Shape: (Ne, nPg, 3, 3).

Return type:

FeArray

Compute_I1() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I1(u)

Returns:

I1_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_I2() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I2(u)

Returns:

I2_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_I3() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I3(u)

Returns:

I3_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_I4(

T: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I4(u)

Parameters:

T (_types.FloatArray) – direction(s)

Returns:

I4_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_I6(

T: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I6(u)

Parameters:

T (_types.FloatArray) – direction(s)

Returns:

I6_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_I8(

T1: ndarray[tuple[Any, ...], dtype[floating]],

T2: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes I8(u)

Parameters:

• T1 (_types.FloatArray) – direction(s) 1

• T2 (_types.FloatArray) – direction(s) 2

Returns:

I8_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_J() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes the deformation gradient J = det(F)

Returns:

J_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_d2I1dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I1dC(u)

Returns:

d2I1dC of shape (d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_d2I2dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I2dC(u)

Returns:

d2I2dC of shape (d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_d2I3dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I3dC(u)

Returns:

d2I3dC_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_d2I4dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I4dC(u)

Returns:

d2I4dC of shape (d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_d2I6dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I6dC(u)

Returns:

d2I6dC of shape (d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_d2I8dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes d2I8dC(u)

Returns:

d2I8dC of shape (d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI1dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI1dC(u)

Returns:

dI1dC of shape (d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI2dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI2dC(u)

Returns:

dI2dC_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI3dC() → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI3dC(u)

Returns:

dI3dC_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI4dC(

T: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI4dC(u)

Parameters:

T (_types.FloatArray) – direction(s)

Returns:

dI4dC_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI6dC(

T: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI6dC(u)

Parameters:

T (_types.FloatArray) – direction(s)

Returns:

dI6dC_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

Compute_dI8dC(

T1: ndarray[tuple[Any, ...], dtype[floating]],

T2: ndarray[tuple[Any, ...], dtype[floating]],

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Computes dI8dC(u)

Parameters:

• T1 (_types.FloatArray) – direction(s) 1

• T2 (_types.FloatArray) – direction(s) 2

Returns:

dI8dC_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_Anisotropic_Invariants(

T1: ndarray[tuple[Any, ...], dtype[floating]],

T2: ndarray[tuple[Any, ...], dtype[floating]],

)

_Compute_Anisotropic_Invariants_First_Derivatives(

T1: ndarray[tuple[Any, ...], dtype[floating]],

T2: ndarray[tuple[Any, ...], dtype[floating]],

)

_Compute_C() → list[FeArray | ndarray[tuple[Any, ...], dtype[Any]]]

Computes the right Cauchy-Green tensor components C(u) = F(u)’.F(u)

returns [cxx, cxy, cxz, cyx, cyy, cyz, czx, czy, czz]

_GetDims() → tuple[int, int, int]

return Ne, nPg, dim

_Get_normalized_components(

T: ndarray[tuple[Any, ...], dtype[floating]],

) → tuple[FeArray | ndarray[tuple[Any, ...], dtype[Any]], FeArray | ndarray[tuple[Any, ...], dtype[Any]], FeArray | ndarray[tuple[Any, ...], dtype[Any]]]

_Slice_Matrix(

matrix: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

)

_Slice_Vector(

vector: FeArray | ndarray[tuple[Any, ...], dtype[Any]],

)

property displacement

displacement field in (xi,yi,zi,…,xn,yn,zn) format.

property groupElem

group of elements.

property matrixType

matrix type.

class EasyFEA.Models.HyperElastic.MooneyRivlin(

dim: int,

K1: float | ndarray[tuple[Any, ...], dtype[floating]],

K2: float | ndarray[tuple[Any, ...], dtype[floating]],

K: float | ndarray[tuple[Any, ...], dtype[floating]] = 0.0,

thickness=1.0,

)

Bases: _HyperElastic

Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

_Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

K: float

Bulk modulus

K1: float

Kappa1

K2: float

Kappa2

class EasyFEA.Models.HyperElastic.NeoHookean(

dim: int,

K: float | ndarray[tuple[Any, ...], dtype[floating]],

thickness=1.0,

)

Bases: _HyperElastic

Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

_Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

K: float

Bulk modulus

class EasyFEA.Models.HyperElastic.SaintVenantKirchhoff(

dim: int,

lmbda: float | ndarray[tuple[Any, ...], dtype[floating]],

mu: float | ndarray[tuple[Any, ...], dtype[floating]],

K: float | ndarray[tuple[Any, ...], dtype[floating]] = 0.0,

thickness=1.0,

)

Bases: _HyperElastic

Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

_Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

_Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

K: float

Bulk modulus

lmbda: float

Lame’s first parameter

mu: float

Shear modulus

class EasyFEA.Models.HyperElastic._HyperElastic(dim: int, thickness: float)

Bases: _IModel, ABC

HyperElastic material.

NeoHookean, MooneyRivlin, SaintVenantKirchhoff and HolzapfelOgden inherit from _HyperElas class.

static Available_Laws()

abstractmethod Compute_W(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the quadratic energy W(u).

Parameters:

hyperElasticState (HyperElasticState) – Hyperelastic state containing the mesh, the discretized field, and the matrix type.

Returns:

We_e_pg of shape (Ne, pg)

Return type:

FeArray

Compute_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Total consistent tangent dΣ/de in Kelvin–Mandel matrix form.

Pass-through to _Compute_elastic_d2Wde() today. Same rationale as Compute_dWde(): the wrapper exists so future contributions with non-zero ∂Σ/∂E can compose here. Kelvin–Voigt viscosity is independent of E and is delivered through the damping matrix Operators.NonLinear.KelvinVoigtDamping().

Compute_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Total PK2 in Kelvin-Mandel vector form — elastic + active stress.

Composition:

1. _Compute_elastic_dWde(state) — subclass-specific elastic part.

2. + active_stress · (T̂ ⊗ T̂) if active_stress > 0 and the direction T̂ has been registered via Set_active_stress_vec() (E-independent active PK2 addend).

The result feeds the residual ∫ Bᵀ · dWde dΩ and the geometric tangent via Sig = block(P(dWde)).

Kelvin–Voigt viscosity does not fold in here — it lives in the separate damping matrix Operators.NonLinear.KelvinVoigtDamping(), mirroring how Rayleigh damping works in Elastic. The simulation handles both the residual contribution (b -= C @ v_t) and the linear tangent piece (coefC · C in the global assembly) — folding η · Ė here would double-count the viscous force.

Set_active_stress_vec(T) → None

Registers the direction of the active PK2 contribution.

Normalises T per Gauss point and precomputes the Kelvin-Mandel vector form of T̂ ⊗ T̂, so that Compute_dWde() returns Σ = Σ_elastic + active_stress · (T̂ ⊗ T̂) whenever active_stress > 0. Only the scalar active_stress is updated between steps; the direction is constant across a run, so call this once during setup.

Note: with a time scheme, set active_stress at the scheme’s effective time (midpoint t + dt/2, HHT t + (1−α)·dt, Newmark/Euler t + dt), not the endpoint t.

Parameters:

T – Direction tensor at every Gauss point, shape (Ne, nPg, 3). Need not be unit-norm — divided by |T| per Gauss point internally.

abstractmethod _Compute_elastic_d2Wde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic consistent tangent dΣ_elastic/de.

Subclass hook — implement the elastic tangent here. The public Compute_d2Wde() wraps this; Kelvin–Voigt viscosity adds no constitutive ∂Σ/∂E term (Σ_visco = η · Ė is independent of E), so the base wrapper is a pass-through today.

Returns:

dΣde_e_pg of shape (Ne, pg, d, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

Return type:

FeArray

abstractmethod _Compute_elastic_dWde(

hyperElasticState: HyperElasticState,

) → FeArray

Computes the elastic second Piola-Kirchhoff tensor Σ_elastic(u).

Subclass hook — implement the elastic kernel here. The public Compute_dWde() wraps this with viscous / future active-stress contributions when applicable.

Returns:

• FeArray – Σ_e_pg of shape (Ne, pg, d), where d = 1, 3, 6 depending on whether the solution dimension is 1D, 2D, or 3D.

• Σxx, Σyy, Σzz, sqrt(2) Σyz, sqrt(2) Σxz, sqrt(2) Σxy

active_stress: float

Active stress magnitude (per-time-step scalar). 0 (default) → inactive. When > 0 and the direction tensor has been registered via Set_active_stress_vec(), the PK2 contribution active_stress · (T̂ ⊗ T̂) is folded into Compute_dWde(). Typical cardiac use: precompute the fiber direction tensor once with Set_active_stress_vec(), then update only this scalar between Solve() calls — material.active_stress = float(tau_values[i]).

property coef: float

kelvin mandel coef -> sqrt(2)

dim: float

eta: float

Kelvin–Voigt viscosity. 0 (default) → purely elastic When > 0 and a velocity field is available, the viscous contribution is delivered via the damping matrix Operators.NonLinear.KelvinVoigtDamping(). The simulation handles both the residual contribution (b -= C @ v_t) and the linear tangent piece (coefC · C in the global assembly) — same uniform pattern as Rayleigh damping in Elastic.

property isHeterogeneous: bool

indicates whether the model has heterogeneous parameters

property modelType: ModelType

model type

thickness: float

Beam module.

class EasyFEA.Models.Beam.BeamStructure(beams: list[_Beam])

Bases: _IModel

Beam structure class.

Calc_D_e_pg(

groupElem: _GroupElem,

matrixType: MatrixType = MatrixType.beam,

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Returns a matrix characterizing the beams’s stiffness behavior.

Calc_M_e_pg(

groupElem: _GroupElem,

) → FeArray | ndarray[tuple[Any, ...], dtype[Any]]

Returns a matrix characterizing the beams’s stiffness behavior.

Get_axis_e(

groupElem: _GroupElem,

) → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Returns the fiber and cross bar vertical axis on every elements.

return xAxis_e, yAxis_e

property areas: list[float]

beams areas

property beams: list[_Beam]

property dim: int

model dimensions

1D -> tension compression

2D -> tension compression + bending + flexion

3D -> all

property dof_n: int

Degrees of freedom per node.

1D -> [u1, … , un]

2D -> [u1, v1, rz1, …, un, vn, rzn]

3D -> [u1, v1, w1, rx1, ry1, rz1, …, u2, v2, w2, rx2, ry2, rz2]

property modelType: ModelType

model type

property nBeam: int

Number of beams in the structure

property thickness: float

The beam structure can have several beams and therefore different sections.

You need to look at the section of the beam you are interested in.

class EasyFEA.Models.Beam.Isotropic(

dim: int,

line: Line,

section: Mesh,

E: float,

v: float,

yAxis: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (0, 1, 0),

)

Bases: _Beam

Isotropic elastic beam.

Get_D(

useTimoshenko: bool = False,

) → ndarray[tuple[Any, ...], dtype[floating]]

Returns a matrix characterizing the beam’s stiffness behavior.

Get_M() → ndarray[tuple[Any, ...], dtype[floating]]

Returns a matrix characterizing the beam’s mass behavior.

E: float

Young’s modulus

property mu: float

shear modulus (G)

v: float

Poisson’s ratio

class EasyFEA.Models.Beam._Beam(

dim: int,

line: Line,

section: Mesh,

yAxis: ndarray[tuple[Any, ...], dtype[Any]] | Collection[int | float] = (0, 1, 0),

)

Bases: _IModel

Beam class model.

abstractmethod Get_D(

useTimoshenko: bool = False,

) → ndarray[tuple[Any, ...], dtype[floating]]

Returns a matrix characterizing the beam’s stiffness behavior.

abstractmethod Get_M() → ndarray[tuple[Any, ...], dtype[floating]]

Returns a matrix characterizing the beam’s mass behavior.

_Calc_P() → ndarray[tuple[Any, ...], dtype[floating]]

P matrix use to transform beam coordinates to global coordinates.

[ix, jx, kx

iy, jy, ky

iz, jz, kz]

coord(x,y,z) = P • coord(i,j,k)

_Get_shear_correction_factor(axis: str = 'y') → float

Cowper’s (1966) shear correction factor k for the cross-section, evaluated at Poisson’s ratio ν = 0.

Used to populate self._ky / self._kz in __init__ — those attributes are what Isotropic.Get_D actually reads, so users can override them with any value (Pure Jouravski, a textbook table value, Cowper-with-ν=v, …) without touching this helper.

Computes k by solving one Saint-Venant Poisson problem on the 2-D section mesh:

∇²φ = -s in the section S ∂φ/∂n = 0 on the boundary ∂S

where s is the section coordinate along the shear direction: axis="y" → s = y (gives k_y, pairs with Iz) axis="z" → s = x (gives k_z, pairs with Iy). The rigid-body (constant) mode of the Neumann problem is fixed by clamping one node to 0; k below is translation-invariant.

Using the Neumann BC, the energy identity gives

uᵀ · f = ∫_S |∇φ|² dS = ∫_S s · φ dS

and Cowper at ν = 0 reduces to

k = I² / (A · uᵀ · f)

with I = Iz for axis=”y”, I = Iy for axis=”z”.

Reference values (any mesh size, any aspect ratio):

rectangle → 5/6 ≈ 0.8333 circle → 6/7 ≈ 0.8571 I-beam, hollow tube, channel, … → whatever the geometry yields

Note vs. the textbook “Pure Jouravski” value:

• The two formulas agree on the rectangle (both 5/6).

• For curved boundaries they differ: Pure Jouravski’s 1-D τ_yz = V·Q/(I·b), τ_xz = 0 assumption gives 9/10 for a circle; this Poisson PDE returns 6/7 because it accounts for the non-zero τ_xz that the curved boundary forces. 6/7 is the physically correct ν = 0 value (matches 3-D elasticity).

Cowper-with-ν ≠ 0 needs a different PDE (non-zero Neumann boundary term involving ν); not implemented here. If you need a specific k value other than Cowper-ν=0 — Pure Jouravski, Cowper-with-ν, a steel-table value — set beam._ky / beam._kz directly.

property Iy: float

squared moment of the cross-section around y-axis.

int_S z^2 dS

property Iz: float

squared moment of the cross-section around z-axis.

int_S y^2 dS

property J: float

polar area moment of inertia (Iy + Iz).

_ky

Shear correction factor for transverse shear in the beam’s y direction (paired with Iz bending in Isotropic.Get_D).

Populated in _Beam.__init__ with the Cowper-ν=0 value returned by _Get_shear_correction_factor("y"). Override directly to use a different k:

# Pure Jouravski (textbook): same as default for rectangle, 9/10 for circle beam._ky = 9 / 10

# Cowper-with-ν for a rectangular section (k = 10(1+ν)/(12+11ν)) beam._ky = 10 * (1 + beam.v) / (12 + 11 * beam.v)

# Cowper-with-ν for a circular section (k = 6(1+ν)/(7+6ν)) beam._ky = 6 * (1 + beam.v) / (7 + 6 * beam.v)

# Value from a steel-section table, etc. beam._ky = 0.49

_kz

Shear correction factor for transverse shear in the beam’s z direction (paired with Iy bending in Isotropic.Get_D). 3-D beams only — unused in 1-D / 2-D simulations.

Populated in _Beam.__init__ with the Cowper-ν=0 value returned by _Get_shear_correction_factor("z"). Override directly to use a different k — same Cowper-with-ν / Pure Jouravski / table-value patterns as _ky, just applied to the z direction (so e.g. for a rectangle bxh, the paired bending inertia is Iy = h·b³/12 and the formula factor refers to the b direction).

property area: float

cross-section area.

property dim

model dimension

property dof_n: int

Degrees of freedom per node 1D -> [u1, … , un]

2D -> [u1, v1, rz1, …, un, vn, rzn]

3D -> [u1, v1, w1, rx1, ry1, rz1, …, u2, v2, w2, rx2, ry2, rz2]

property line: Line

average fiber line of the beam.

property modelType: ModelType

model type

property name: str

beam name/tag

property section: Mesh

beam cross-section in (x,y) plane

property thickness: float

thickness used in the model

property xAxis: ndarray[tuple[Any, ...], dtype[floating]]

perpendicular cross-beam axis (fiber)

property yAxis: ndarray[tuple[Any, ...], dtype[floating]]

vertical cross-beam axis