Simulations

The EasyFEA.Simulations module provides essential tools for creating and managing simulations. These simulations are built using a Mesh and a _IModel (see Models).

In the simulation workflow, Simulations is the central step: it takes a mesh and a model, exposes boundary-condition methods (add_dirichlet, add_surfLoad, …), drives the linear solver, and stores the solution history.

With this module, you can construct:

• Linear elastic simulations with Elastic.

• Nonlinear hyperelastic simulations with HyperElastic.

• Euler-Bernoulli beam simulations with Beam.

• PhaseField damage simulations for quasi-static brittle fracture with PhaseField.

• Thermal simulations with Thermal.

• Weak form simulations with WeakForms.

See also

• Apply boundary conditions

• Understand the solve pipeline

Matrix System Solvers

EasyFEA automatically manages the resolution of elliptic, parabolic, and hyperbolic matrix systems, allowing developers to focus exclusively on constructing local matrices via the Construct_local_matrix_system method.

Elliptic

(1)\[ \Krm \, \urm = \Frm \]

Parabolic

(2)\[ \Krm \, \urm^{n+\alpha} + \Crm \, \vrm^{n+\alpha} = \Frm^{n+\alpha} \]

Set with simu.Solver_Set_Parabolic_Algorithm(dt, alpha=0.5).

Method	α	Order	Stability
Forward Euler	0	1st	Conditionally stable
Crank–Nicolson	0.5	2nd	Unconditionally stable
Backward Euler	1	1st	Unconditionally stable

Hyperbolic

(3)\[ \Krm \, \urm + \Crm \, \vrm + \Mrm \, \arm = \Frm \]

Set with simu.Solver_Set_Hyperbolic_Algorithm(dt, algo=AlgoType.newmark).

Method	AlgoType	Order	Stability	Notes
Newmark β	newmark	2nd	Unconditionally stable	Default; energy-conserving (β=1/4, γ=1/2)
Midpoint	midpoint	2nd	Unconditionally stable	Energy-conserving
HHT-α	hht	2nd	Unconditionally stable	Numerical damping (α ∈ [0, 1[)
Euler implicit	euler_implicit	1st	Unconditionally stable	Dissipative
Euler explicit	euler_explicit	1st	Conditionally stable (dt < h_e/c)	Linear only

How to Create New Simulations in EasyFEA?

To create new simulation classes, you can take inspiration from existing implementations. Make sure to follow the _Simu interface. The Thermal class is relatively simple and can serve as a good starting point.

See also

• Create a custom simulation

• EasyFEA/Simulations/_thermal.py source code

• Understand the solve pipeline

Simulations API

This module contains all available simulation classes.

class EasyFEA.Simulations.Beam(

mesh: Mesh,

model: BeamStructure,

folder: str = '',

verbosity=False,

useTimoshenko: bool = False,

)

Bases: _Simu

Beam simulation (Euler-Bernoulli or Timoshenko).

Timoshenko elements use pure Lagrange interpolation for every kinematic field with selective reduced integration on the shear term — the standard formulation used by Abaqus B31/B32, ANSYS BEAM188/189, Calculix B31 and FEniCS dolfinx. As in those codes, accuracy depends on the polynomial degree of the element:

elem | v interpolation | tip-load convergence | accurate at coarse mesh? ——-+—————–+———————-+————————- SEG2 | linear | O(h²) | no — needs fine mesh SEG3 | quadratic | exact (cubic v) | yes SEG4 | cubic | exact | yes SEG5 | quartic | exact | yes

SEG2+Timoshenko remains a valid choice (it matches industry-standard 2-node Timoshenko elements), but for results comparable to the Euler- Bernoulli SEG2 element (which is exact for cubic v), prefer SEG3 or higher when activating Timoshenko.

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Get_Iteration_Summary() → str

Returns the iteration’s summary.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(details=False) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'displacement': array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.23099308e-01, 1.25163294e-03, -5.45209754e-01, 5.61452028e-03, -8.85989704e-03, 2.54610088e-03, 1.23099308e-01, 1.25163294e-03, -5.45209754e-01, 5.61452028e-03, -8.85989704e-03, 2.54610088e-03, 1.23047581e-01, 1.90036372e-01, 4.62082368e-02, 5.61452028e-03, -1.02319399e-02, 3.37444122e-03, 1.11123892e-03, 1.45329487e-04, -1.31196704e-02, 1.89978163e-03, -9.84433005e-04, 1.68127259e-04, 4.52854098e-03, 2.91302287e-04, -4.88544993e-02, 3.39351094e-03, -1.96886601e-03, 3.46609068e-04, 1.04303616e-02, 4.38136701e-04, -1.01867041e-01, 4.49477685e-03, -2.95329901e-03, 5.42363904e-04, 1.91149957e-02, 5.85279934e-04, -1.67088408e-01, 5.23052629e-03, -3.93773202e-03, 7.65888380e-04, 3.10254663e-02, 7.31222884e-04, -2.39910485e-01, 5.64358034e-03, -4.92216502e-03, 1.02790675e-03, 4.67349035e-02, 8.73558778e-04, -3.16426059e-01, 5.79683933e-03, -5.90659803e-03, 1.33678575e-03, 6.68963168e-02, 1.00935740e-03, -3.93707308e-01, 5.77389137e-03, -6.89103103e-03, 1.69564395e-03, 9.21591406e-02, 1.13588759e-03, -4.70056980e-01, 5.67501394e-03, -7.87546404e-03, 2.10098540e-03, 2.75455672e-04, 7.26194756e-05, -3.39361161e-03, 1.00102142e-03, -4.92216502e-04, 8.29849030e-05, 2.52266933e-03, 2.18207541e-04, -2.84980014e-02, 2.69685194e-03, -1.47664951e-03, 2.55765353e-04, 7.15357189e-03, 3.64625556e-04, -7.35262552e-02, 3.99184616e-03, -2.46108251e-03, 4.41712430e-04, 1.44005182e-02, 5.11741398e-04, -1.33252942e-01, 4.90608730e-03, -3.44551552e-03, 6.49957645e-04, 2.46338360e-02, 6.58533896e-04, -2.02813204e-01, 5.47396872e-03, -4.42994852e-03, 8.91477319e-04, 3.83655034e-02, 8.03019435e-04, -2.77915763e-01, 5.74790714e-03, -5.41438153e-03, 1.17613487e-03, 5.62173753e-02, 9.42461013e-04, -3.55108421e-01, 5.80146523e-03, -6.39881453e-03, 1.51004825e-03, 7.88519112e-02, 1.07391662e-03, -4.32050913e-01, 5.72711909e-03, -7.38324754e-03, 1.89291040e-03, 1.06887123e-01, 1.19512334e-03, -5.07738593e-01, 5.63230832e-03, -8.36768054e-03, 2.31899592e-03, 1.23079260e-01, 3.40813982e-02, -4.35107134e-01, 5.61452028e-03, -9.46243345e-03, 2.90725625e-03, 1.23064768e-01, 7.05906867e-02, -3.18855017e-01, 5.61452028e-03, -9.88122172e-03, 3.15947868e-03, 1.23055020e-01, 1.09481470e-01, -1.98678990e-01, 5.61452028e-03, -1.01206730e-02, 3.30550829e-03, 1.23049411e-01, 1.49573139e-01, -7.65319486e-02, 5.61452028e-03, -1.02175422e-02, 3.36476293e-03, 1.23088534e-01, 1.71249178e-02, -4.91064004e-01, 5.61452028e-03, -9.18228496e-03, 2.74037562e-03, 1.23071374e-01, 5.19577064e-02, -3.77608992e-01, 5.61452028e-03, -9.69562556e-03, 3.04698087e-03, 1.23059344e-01, 8.98176403e-02, -2.59125004e-01, 5.61452028e-03, -1.00215381e-02, 3.24511553e-03, 1.23051726e-01, 1.29439638e-01, -1.37748938e-01, 5.61452028e-03, -1.01838427e-02, 3.34389438e-03, 1.23048036e-01, 1.69791753e-01, -1.51807437e-02, 5.61452028e-03, -1.02301569e-02, 3.37311816e-03])},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

_Calc_Epsilon_e_pg(

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

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

Construct deformations for each element and each Gauss point.

a’ denotes here da/dx

1D -> [ux’]

2D -> [ux’, rz’]

3D -> [ux’, rx’, ry’, rz’]

_Calc_InternalForces_e_pg(

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

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

Calculation of internal forces.

1D -> [N]

2D -> [N, Mz]

3D -> [N, Mx, My, Mz]

_Calc_Sigma_e_pg(

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

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

Calculates stresses from strains.

1D -> [Sxx]

2D -> [Sxx, Syy, Sxy]

3D -> [Sxx, Syy, Szz, Syz, Sxz, Sxy]

_Check_dim_mesh_material() → None

Checks that the material dim matches the mesh dim.

_indexResult(result: str) → int

add_connection(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list[str],

description: str,

)

Connects beams together in the specified unknowns.

Parameters:

• nodes (_types.IntArray) – nodes

• unknowns (list[str]) – unknowns

• description (str) – description

add_connection_fixed(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

description='Fixed',

)

Adds a fixed connection.

Parameters:

• nodes (_types.IntArray) – nodes

• description (str, optional) – description, by default “Fixed”

add_connection_hinged(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

unknowns=[''],

description='Hinged',

)

Adds a hinged connection.

Parameters:

• nodes (_types.IntArray) – nodes

• unknowns (list, optional) – unknowns, by default [‘’]

• description (str, optional) – description, by default “Hinged”

add_lineLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds a linear load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Distributed load of -1000 N/m along y (elastic simulation):

>>> simu.add_lineLoad(nodes, [-1000], ["y"])

Spatially varying load (uses integration-point coordinates):

>>> simu.add_lineLoad(nodes, [lambda x, y, z: -500 * (1 + x)], ["y"])

add_surfLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list,

problemType=None,

description='',

)

Adds a surface load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Uniform pressure of -800 Pa along x:

>>> simu.add_surfLoad(nodes, [-800], ["x"])

add_volumeLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list,

problemType=None,

description='',

)

Adds a volumetric load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Gravity in -y direction, ρ = 7800 kg/m³ (elastic simulation):

>>> simu.add_volumeLoad(mesh.nodes, [-7800 * 9.81], ["y"])

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

Center of mass / barycenter / inertia center

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

Displacement vector field.

1D [uxi, …]

2D [uxi, uyi, rzi, …]

3D [uxi, uyi, uzi, rxi, ryi, rzi, …]

property inDim: int

property mass: float

property structure: BeamStructure

Beam structure.

useTimoshenko: bool

class EasyFEA.Simulations.DIC(

mesh: Mesh,

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

lr: float = 0.0,

verbosity: bool = False,

)

Bases: _IObserver

Digital Image Correlation (DIC) analysis

static Get_Circle(

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

threshold: float,

boundary=None,

radiusCoef=1.0,

) → tuple[float, float, float]

Returns the circle properties in the image.

Parameters:

• img (_types.FloatArray) – image used

• threshold (float) – threshold for pixel color

• boundary (tuple[tuple[float, float], tuple[float, float]], optional) – ((xMin, xMax),(yMin, yMax)), by default None

• radiusCoef (float, optional) – multiplier coef for radius, by default 1.0

Returns:

circle coordinates and radius

Return type:

XC, YC, radius

static _Get_coords(

img: ndarray,

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

Calc_pixelDisplacement(

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

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

Computes pixel displacements based on the finite element displacement field.

Calc_r_dic(

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

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

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

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

Computes the dic residual between img and imgRef (as a Np x Np matrix).

r_dic = f(x) - g(x + u(x))

Parameters:

• u (_types.FloatArray) – finite element displacement field (Ndof)

• img (_types.FloatArray) – deformed image

• imgRef (_types.FloatArray, optional) – reference image, by default None

Returns:

the dic residual between images (as a Np x Np matrix)

Return type:

_types.FloatArray

Get_scaled_mesh(scale: float = 1.0) → Mesh

Save(folder: str, filename: str = 'dic') → None

Saves the dic analysis in folder as ‘filename.pickle’.

Solve(

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

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

iterMax: int = 1000,

tolConv: float = 1e-06,

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

verbosity=True,

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

Computes the displacement field between the two images.

Parameters:

• img (_types.FloatArray) – deformed image (g)

• u0 (_types.FloatArray, optional) – initial displacement field, by default None If None, the field is initialized using the _Get_u_from_images(imgRef, img) function if opencv-python is installed; otherwise, it is initialized to zeros.

• iterMax (int, optional) – maximum number of iterations, by default 1000

• tolConv (float, optional) – convergence tolerance (converged once ||b|| <= tolConv), by default 1e-6

• imgRef (_types.FloatArray, optional) – reference image (f), by default None

• verbosity (bool, optional) – display iterations, by default True

Returns:

computed displacement field (Ndof)

Return type:

_types.FloatArray

_Compute_L_M(

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

) → None

Computes DIC matrices.

_Get_ldic(coef: float = 8.0) → float

Get the characteristic length of the mesh, i.e. coef * the average mesh size.

_Get_u_from_images(

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

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

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

Use open cv to calculate displacements between images.

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

Returns the 2D periodic vector field.

_Update(observable: Observable, event: str) → None

Receive an update/event from an observable object (observer pattern).

property ROI: ndarray[tuple[Any, ...], dtype[int64]]

roi as a matrix.

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

reference image (f)

Warning

Updating _imgRef will automatically update the roi, initialize the shape functions and construct the laplacian operator!

property _lr: float

regularization length

property alpha: float

property ldic: float

8 * mean(meshSize) (see Get_ldic())

property mesh: Mesh

pixel-based mesh used for dic.

property roi: ndarray[tuple[Any, ...], dtype[int64]]

roi as a vector.

property shape: tuple[int, int]

reference image shape

class EasyFEA.Simulations.Elastic(

mesh: Mesh,

model: _Elastic,

folder: str = '',

verbosity=False,

)

Bases: _Simu

Linearized elasticity.

Strong form:

\[\begin{split}\diver{\Sig(\ub)} + \fb &= \rho \, \ddot{\ub} && \quad \text{in } \Omega, \\ % \Sig(\ub) \cdot \nb &= \tb && \quad \text{on } \partial\Omega_t, \\ % \Sig(\ub) &= \Cbb : \Eps(\ub) && \quad \text{in } \Omega, \\ % \ub &= \ub_D && \quad \text{on } \partial\Omega_u,\end{split}\]

Weak form:

\[\int_\Omega \Sig(\ub) : \Eps(\vb) \, \dO + \int_\Omega \rho \, \ddot{\ub} \cdot \vb \, \dO = \int _{\partial\Omega_t} \tb\cdot\vb \, \dS + \int _{\Omega} \fb\cdot\vb \, \dO \quad \forall \, \vb \in V\]

The implemented elastic laws are available here.

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Get_Iteration_Summary() → str

Returns the iteration’s summary.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(details=False) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'accel': array([0., 0., 0., ..., 0., 0., 0.], shape=(6615,)), 'displacement': array([2.06877508e-04, -2.05036557e-04, 5.26255295e-08, ..., -2.10638322e-04, -9.99171321e-05, 4.27128024e-07], shape=(8049,)), 'speed': array([0., 0., 0., ..., 0., 0., 0.], shape=(6615,))},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

Set_Rayleigh_Damping_Coefs(coefM=0.0, coefK=0.0)

Sets damping coefficients ( C = coefK * K + coefM * M ).

_Calc_Epsilon_e_pg(

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

matrixType=MatrixType.rigi,

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

Computes strain field from the displacement vector field.

2D : [Exx Eyy sqrt(2)*Exy]

3D : [Exx Eyy Ezz sqrt(2)*Eyz sqrt(2)*Exz sqrt(2)*Exy]

Parameters:

u (_types.FloatArray) – displacement vector (Ndof)

Returns:

Computed strain field (Ne,pg,(3 or 6))

Return type:

FeArray

_Calc_Psi_Elas(

returnScalar=True,

smoothedStress=False,

matrixType=MatrixType.rigi,

)

Computes the kinematically admissible deformation energy. Wdef = 1/2 int_Ω Sig : Eps dΩ

_Calc_Sigma_e_pg(

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

matrixType=MatrixType.rigi,

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

Computes stress field from strain field.

2D : [Sxx Syy sqrt(2)*Sxy]

3D : [Sxx Syy Szz sqrt(2)*Syz sqrt(2)*Sxz sqrt(2)*Sxy]

Parameters:

Epsilon_e_pg (FeArray.FeArrayALike) – Strain field (Ne,pg,(3 or 6))

Returns:

Computed stress field (Ne,pg,(3 or 6))

Return type:

FeArray

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

Computes the ZZ1 error.

For more details, [F.Pled, Vers une stratégie robuste … ingénierie mécanique] page 20/21

Returns the global error and the error on each element.

Return type:

error, error_e

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

Acceleration vector field.

2D [axi, ayi, …]

3D [axi, ayi, azi, …]

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

Displacement vector field.

2D [uxi, uyi, …]

3D [uxi, uyi, uzi, …]

property material: _Elastic

elastic material

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

Velocity vector field.

2D [vxi, vyi, …]

3D [vxi, vyi, vzi, …]

class EasyFEA.Simulations.HyperElastic(

mesh: Mesh,

model: _HyperElastic,

folder: str = '',

tolConv=1e-05,

maxIter=20,

verbosity=False,

)

Bases: _Simu

Hyperelastic simulation.

Weak form:

\[R(\ub; \vb) = \int_{\Omega_0} \boldsymbol{\Sigma}(\ub) : \Drm_\ub \eb(\ub) \cdot \vb \, \dO + \int_{\Omega_0} \rho \, \ddot{\ub} \cdot \vb \, \dO - \int_{\partial\Omega_0^t} \tb\cdot\vb \, \dS - \int_{\Omega_0} \fb\cdot\vb \, \dO \quad \forall \, \vb \in V\]

where \(\boldsymbol{\Sigma} := J \, \Fb^{-1} \cdot \Sig \cdot \Fb^{-T}\), is the second Piola Kirchhoff stress tensor (PK2), \(\eb := \frac{1}{2} \left( \Cb - \boldsymbol{1} \right) = \frac{1}{2} \left( \Fb^T \cdot \Fb - \boldsymbol{1} \right)\) is the Green-Lagrange strain tensor and \(\Fb := \boldsymbol{1} + \grad \ub\) the deformation gradient.

This non linear problem is solve using the newton rapshon algorithm:

\[A(\ub; \vb, \wb) \, \Delta \ub = - R(\ub; \vb) \quad \forall \, (\vb, \wb) \in \Vc \times \Wc,\]

where the tangent \(A(\ub; \vb, \wb)\) is defined \(\forall \, (\vb, \wb) \in \Vc \times \Wc\) as:

\[\begin{split}A(\ub; \vb, \wb) &= \dpartial{R(\ub; \vb)}{\ub} \cdot \wb \\ &= \int_{\Omega_0} \Drm_\ub \eb(\ub) \cdot \wb : \dNpartial{2}{W}{\eb}(\ub) : \Drm_\ub \eb(\vb) \, \dO + \int_{\Omega_0} \dpartial{W}{\eb}(\ub) : \Drm_\ub^2 \eb(\vb, \wb) \, \dO\end{split}\]

The implemented hyperelastic laws are available here and were constructed by the ComputeHyperelasticLaws script.

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_problemTypes()

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_dict_Energy()

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(details=False)

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'accel': array([0., 0., 0., ..., -0.06848857, -1.61017255, 0.08481901], shape=(1344,)), 'displacement': array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ..., -8.85908008e-01, 1.09300780e+01, -4.97057301e-06], shape=(1344,)), 'speed': array([0., 0., 0., ..., -0.18831023, 1.65846028, 0.00308314], shape=(1344,))},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter=-1, resetAll=False)

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

_Calc_GreenLagrange(matrixType=MatrixType.rigi)

_Calc_SecondPiolaKirchhoff(matrixType=MatrixType.rigi)

_Calc_W(returnScalar=True, matrixType=MatrixType.rigi)

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

Displacement vector field.

[uxi, uyi, uzi, …]

property material: _HyperElastic

hyperelastic material

class EasyFEA.Simulations.PhaseField(

mesh: Mesh,

model: PhaseField,

folder: str = '',

verbosity=False,

)

Bases: _Simu

PhaseField damage simulations for quasi-static brittle fracture.

Strong form:

Damaged linear elastic problem

\[\begin{split}-\diver{\Sig(\ub, \phi)} &= \fb && \quad \text{in } \Omega, \\ % \Sig(\ub, \phi) \cdot \nb &= \tb && \quad \text{on } \partial\Omega_t, \\ % \Sig(\ub, \phi) &= \Cbb(\phi) : \Eps(\ub) && \quad \text{in } \Omega, \\ % \ub &= \ub && \quad \text{on } \partial\Omega_u,\end{split}\]

Damage problem

\[\begin{split}- \nabla \cdot \left( \dfrac{2 \, G_c \, \ell}{c_w} \, \nabla\phi \right) + \dfrac{G_c}{c_w \, \ell} \, w'(\phi) &= Y(\Eps, \phi) && \quad \text{in } \Omega, \\ % \nabla \phi \cdot \nb &= 0 && \quad \text{on } \partial\Omega,\end{split}\]

Weak form:

Damaged linear elastic problem

\[\int_\Omega \Sig(\ub, \phi) : \Eps(\vb) \, \dO = \int _{\partial\Omega_t} \tb\cdot\vb \, \dS + \int _{\Omega} \fb\cdot\vb \, \dO \quad \forall \, \vb \in V\]

Damage problem

\[\int_\Omega k_w \, \nabla \phi \cdot \nabla \delta \phi + r_w \, \phi \, \delta \phi \, \dO = \int_\Omega f_w \, \delta \phi \, \dO \quad \forall \, \delta \phi \in V,\]

Further Reading

See section 3.1. of https://univ-eiffel.hal.science/hal-05115523 for additional mathematical details.

static Folder(

folder: str,

material: str,

split: str,

regu: str,

simpli2D: str,

tolConv: float,

solver: str,

test: bool,

optimMesh=False,

closeCrack=False,

nL=0,

theta=0.0,

) → str

Creates a phase field folder based on the specified arguments.

Bc_dofs_nodes(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list[str],

problemType=ModelType.elastic,

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

Returns degrees of freedom associated with the nodes, based on the problem type and unknowns.

Parameters:

• nodes (_types.IntArray) – nodes.

• unknowns (list) – unknowns (e.g [“x”,”y”,”rz”])

• problemType (str) – Problem type.

Returns:

Degrees of freedom.

Return type:

_types.IntArray

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Detect_Damage(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

threshold: float = 1.0,

) → bool

Returns True if damage >= threshold on any node in nodes, across all MPI ranks.

Parameters:

• nodes (IntArray) – Global node indices to check. May be empty.

• threshold (float, optional) – Damage threshold, by default 1.0.

Returns:

True if any node has damage >= threshold (reduced across all MPI ranks).

Return type:

bool

Examples

Stop the load loop when the crack reaches the boundary 10 consecutive times:

>>> nDetect = 0
>>> if simu.Detect_Damage(nodes):
...     nDetect += 1
...     if nDetect == 10:
...         break

Get_K_C_M_F(

problemType=None,

) → tuple[csr_matrix, csr_matrix, csr_matrix, csr_matrix]

Returns the assembled matrices of \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_lb_ub(

problemType: ModelType,

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

Returns the lower bound and upper bound.

Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Need_Update(value=True) → None

Sets whether the simulation needs to reconstruct matrices K, C, M and F.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

) → ndarray[tuple[Any, ...], dtype[floating]] | float | None

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Get_Bc_Summary() → str

Returns the simulation loading summary.

Results_Get_Iteration_Summary() → str

Returns the iteration’s summary.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_Set_Bc_Summary(config: str)

Sets the simulation loading summary.

Results_Set_Iteration_Summary(

iter: int,

load: float,

unitLoad: str,

percentage=0.0,

remove=False,

) → str

Creates the iteration summary for the damage problem.

Parameters:

• iter (int) – iteration

• load (float) – loading

• unitLoad (str) – loading unit

• percentage (float, optional) – percentage of simualtion performed, by default 0.0

• remove (bool, optional) – removes line from terminal after display, by default False

Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(

details=False,

) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'Niter': 1, 'convIter': np.float64(0.0014931281647901395), 'damage': array([0.02516667, 0.99880179, 0.01021515, ..., 0.15390196, 0.23204452, 0.41399389], shape=(1100,)), 'displacement': array([6.94347839e-08, -3.59175015e-10, -1.96950975e-09, ..., 8.02869807e-06, -1.40620306e-09, 8.02838724e-06], shape=(2200,)), 'timeIter': np.float64(0.019930601119995117)},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

Solve(

tolConv=1.0,

maxIter=500,

convOption=2,

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

Solves the iterative damage problem using the staggered scheme.

Parameters:

• tolConv (float, optional) – threshold used to check convergence (𝜖), by default 1.0

• maxIter (int, optional) – Maximum iterations for convergence, by default 500

• convOption (int, optional) –

  • 0 -> convergence on damage np.max(np.abs(d_np1-dk)) equivalent normInf(d_np1-dk)

  • 1 -> convergence on crack energy np.abs(psi_crack_n - psi_crack_np1)/psi_crack_np1

  • 2 -> convergence on total energy np.abs(psi_tot_n - psi_tot_np1)/psi_tot_np1

  eq (39) Ambati 2015 10.1007/s00466-014-1109-y

  • 3 -> (convD <= tolConv) and (convU <= tolConv*0.999)

  eq (25) Pech 2022 10.1016/j.engfracmech.2022.108591

Returns:

u_np1, d_np1, converged

such that:

• u_np1: displacement vector field

• d_np1: damage scalar field

• converged: the solution has converged

Return type:

_types.FloatArray, _types.FloatArray, csr_matrix, bool

_Calc_Epsilon_e_pg(

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

matrixType=MatrixType.rigi,

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

Computes strain field (Ne,pg,(3 or 6)).

2D : [Exx Eyy sqrt(2)*Exy]

3D : [Exx Eyy Ezz sqrt(2)*Eyz sqrt(2)*Exz sqrt(2)*Exy]

Parameters:

sol (_types.FloatArray) – Displacement vector

Returns:

Computed strain field (Ne,pg,(3 or 6))

Return type:

FeArray

_Calc_Psi_Crack() → float

Computes crack’s energy.

_Calc_Psi_Elas() → float

Computes of the kinematically admissible damaged deformation energy.

Psi_Elas = 1/2 int_Ω Sig : Eps dΩ

_Calc_Psi_Ext(

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

) → float

Computes external’s energy.

_Calc_Sigma_e_pg(

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

matrixType=MatrixType.rigi,

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

Computes stress field from strain field.

2D : [Sxx Syy sqrt(2)*Sxy]

3D : [Sxx Syy Szz sqrt(2)*Syz sqrt(2)*Sxz sqrt(2)*Sxy]

Parameters:

Epsilon_e_pg (FeArray.FeArrayALike) – Strain field (Ne,pg,(3 or 6))

Returns:

Computed damaged stress field.

Return type:

FeArray

_Update(

observable: Observable,

event: str,

) → None

Receive an update/event from an observable object (observer pattern).

add_dirichlet(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=ModelType.elastic,

description='',

)

Adds Dirichlet’s boundary conditions.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contains floats, arrays or functions or functions.

  e.g [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  The functions use the x, y and z nodes coordinates.

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Fix all displacement DOFs to zero on the left wall (elastic simulation):

>>> simu.add_dirichlet(nodes, [0, 0], ["x", "y"])

Impose temperature 100°C on the top boundary (thermal simulation):

>>> simu.add_dirichlet(nodes, [100.0], ["t"])

Spatially varying temperature (uses node coordinates):

>>> simu.add_dirichlet(nodes, [lambda x, y, z: 100 * x], ["t"])

add_lineLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=ModelType.elastic,

description='',

)

Adds a linear load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Distributed load of -1000 N/m along y (elastic simulation):

>>> simu.add_lineLoad(nodes, [-1000], ["y"])

Spatially varying load (uses integration-point coordinates):

>>> simu.add_lineLoad(nodes, [lambda x, y, z: -500 * (1 + x)], ["y"])

add_neumann(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=ModelType.elastic,

description='',

)

Adds Neumann’s boundary conditions.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contains floats, arrays or functions or functions.

  e.g [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  The functions use the x, y and z nodes coordinates.

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Concentrated force of -500 N in y at the beam tip (beam simulation):

>>> nodes = mesh.Nodes_Point((L, h/2))
>>> simu.add_neumann(nodes, [-500], ["y"])

add_pressureLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

magnitude: float,

problemType=ModelType.elastic,

description='',

) → None

Adds a pressure.

Parameters:

• nodes (_types.IntArray) – nodes. (must belong to the edge of the mesh.)

• magnitude (float) – pressure magnitude

• problemType (str, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Normal inward/outward pressure of 1 MPa on a cylinder surface (elastic simulation):

>>> simu.add_pressureLoad(cylinder_nodes, magnitude=1e6)

add_surfLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=ModelType.elastic,

description='',

)

Adds a surface load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Uniform pressure of -800 Pa along x:

>>> simu.add_surfLoad(nodes, [-800], ["x"])

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

Damage scalar field.

[di, …]

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

Displacement vector field.

2D [uxi, uyi, …]

3D [uxi, uyi, uzi, …]

property needUpdate: bool

The object needs to be updated.

property phaseFieldModel: PhaseField

damage model

class EasyFEA.Simulations.Thermal(

mesh: Mesh,

model: Thermal,

folder: str = '',

verbosity=False,

)

Bases: _Simu

Thermal simulation.

Strong form:

\[\begin{split}\diver q + \rho \, c \, \dot{t} &= r && \quad \text{in } \Omega, \\ % q &= -k \, \nabla t && \quad \text{in } \Omega, \\ % q \cdot \nb &= -q_n && \quad \text{on } \partial\Omega_q, \\ % t &= t_D && \quad \text{on } \partial\Omega_t,\end{split}\]

Weak form:

\[\int_\Omega k \, \nabla t \cdot \nabla \delta t \, \dO + \int_\Omega \rho \, c \, \dot{t} \, \delta t \, \dO = \int_{\partial\Omega_q} q_n \, \delta t \, \dS + \int_{\partial\Omega_t} r \, \delta t \, \dO \quad \forall \, \delta t \in V\]

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(details=False) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'thermal': array([2.75928312, 2.75928312, 2.75928312, ..., -1.61828253, -1.61828253, -1.61828253], shape=(1584,)), 'thermalDot': array([1.26156081, 1.26156081, 1.26156081, ..., 1.21148448, 1.21148448, 1.21148448], shape=(1584,))},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

_Check_dim_mesh_material() → None

Checks that the material dim matches the mesh dim.

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

Scalar temperature field.

[ti, ….]

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

Time derivative of the scalar temperature field.

[d(ti)/dt, ….]

property thermalModel: Thermal

Thermal simulation model.

class EasyFEA.Simulations.WeakForms(

mesh: Mesh,

model: WeakForms,

folder: str = '',

isNonLinear=False,

tolConv=1e-05,

maxIter=20,

verbosity=False,

)

Bases: _Simu

Construct_local_matrix_system(problemType)

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

Get_x0(problemType=None)

Returns the solution from the previous iteration.

Result(

result: str,

nodeValues: bool = True,

iter: int | None = None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

Results_nodeFields_elementFields(details=False) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save_Iter(

iter={'a': array([0., 0., 0., ..., -3568.98349987, 698261.93869521, -686178.95162937], shape=(1344,)), 'u': array([0., 0., -0.63468269, ..., -2.57176296, 0.58676121, -2.57181197], shape=(3782,)), 'v': array([0., 0., 0., ..., 1.11176005, -22.89989076, -22.73918473], shape=(1344,))},

)

Saves iteration results in _results or in the specified simu.folder.

Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

_Check_dim_mesh_material() → None

Checks that the material dim matches the mesh dim.

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

node field a = d2udt2

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

node field u.

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

node field v = dudt

property weakForms: WeakForms

Weak form manager.

class EasyFEA.Simulations._Simu(

mesh: Mesh,

model: _IModel,

folder: str = '',

verbosity: bool = True,

)

Bases: _IObserver, Updatable, ABC

The following classes inherit from the parent class _Simu:

• Elastic

• HyperElastic

• PhaseField

• Beam

• Thermal

• WeakForm

To create new simulation classes, you can take inspiration from existing implementations.n Make sure to follow this interface.n The ThermalSimuclass is relatively simple and can serve as a good starting point.n See simulations/_thermal.py for more details.

To use the interface/inheritance, 13 methods need to be defined.

General:

• def Get_problemTypes(self) -> list[ModelType]:

• def Get_unknowns(self, problemType=None) -> list[str]:

• def Get_dof_n(self, problemType=None) -> int:

These functions provides access to the available unknowns and degrees of freedom (dofs).

Solve:

• def Get_x0(self, problemType=None):

• def Construct_local_matrix_system(self, problemType):

These functions assemble the matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\).

Iterations:

• def Save_Iter(self, iter: dict[str, Any]={}) -> None:

• def Set_Iter(self, index=-1) -> None:

These functions are used to save or load iterations.

Results:

• def Results_Available(self) -> list[str]:

• def Result(self, result: str, nodeValues=True, iter=None) -> float | _types.FloatArray:

• def Results_Iter_Summary(self) -> tuple[list[int], list[tuple[str, _types.FloatArray]]]:

• def Results_dict_Energy(self) -> dict[str, float]:

• def Results_displacement_matrix(self) -> _types.FloatArray:

• def Results_nodesField_elementsField(self, details=False) -> tuple[list[str], list[str]]:

These functions are used to process the results.

Assembly(

problemType: ModelType,

) → tuple[csr_matrix, csr_matrix, csr_matrix, csr_matrix]

Assemble the matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problemType.

Bc_Init() → None

Initializes Dirichlet, Neumann and Lagrange boundary conditions

Bc_dofs_Dirichlet(

problemType=None,

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

Returns dofs related to Dirichlet conditions.

Bc_dofs_Neumann(

problemType=None,

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

Returns dofs related to Neumann conditions.

Bc_dofs_known_unknown(

problemType: ModelType,

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

Returns known and unknown dofs.

Bc_dofs_nodes(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list[str],

problemType=None,

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

Returns degrees of freedom associated with the nodes, based on the problem type and unknowns.

Parameters:

• nodes (_types.IntArray) – nodes.

• unknowns (list) – unknowns (e.g [“x”,”y”,”rz”])

• problemType (str) – Problem type.

Returns:

Degrees of freedom.

Return type:

_types.IntArray

Bc_values_Dirichlet(

problemType=None,

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

Returns dofs values related to Dirichlet conditions.

Bc_values_Neumann(

problemType=None,

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

Returns dofs values related to Neumann conditions.

Bc_vector_Dirichlet(

problemType=None,

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

Returns a vector filled with Dirichlet boundary conditions values.

Bc_vector_Neumann(

problemType=None,

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

Returns a vector filled with Neuman boundary conditions values.

Calc_Reaction(

dofs: ndarray[tuple[Any, ...], dtype[int64]] = None,

problemType=None,

) → ndarray

Resultant reaction at nodes: K[dofs, :] @ u, MPI-safe.

abstractmethod Construct_local_matrix_system(

problemType,

) → tuple[FeArray | None, FeArray | None, FeArray | None, FeArray | None]

Construct the local matrix system \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_K_C_M_F(

problemType=None,

) → tuple[csr_matrix, csr_matrix, csr_matrix, csr_matrix]

Returns the assembled matrices of \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\) for the given problem.

Get_contact(

masterMesh: Mesh,

slaveNodes: ndarray[tuple[Any, ...], dtype[int64]] | None = None,

masterNodes: ndarray[tuple[Any, ...], dtype[int64]] | None = None,

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

Returns the simulation nodes detected in the master mesh with the associated displacement matrix to the interface.

Parameters:

• masterMesh (Mesh) – master mesh

• slavenodes (_types.IntArray, optional) – slave nodes, by default None

• masternodes (_types.IntArray, optional) – master nodes, by default None

Returns:

nodes, displacementMatrix

Return type:

tuple[_types.IntArray, _types.FloatArray]

abstractmethod Get_dof_n(problemType=None) → int

Returns the number of degrees of freedom per node.

Get_dofs(problemType=None)

Returns (owned) dofs associated with the problem type.

Get_lb_ub(

problemType: ModelType,

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

Returns the lower bound and upper bound.

abstractmethod Get_problemTypes() → list[ModelType]

Returns the problem types available through the simulation.

abstractmethod Get_unknowns(problemType=None) → list[str]

Returns a list of unknowns available in the simulation.

abstractmethod Get_x0(

problemType: ModelType | None = None,

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

Returns the solution from the previous iteration.

Need_Update(value=True)

Sets whether the simulation needs to reconstruct matrices K, C, M and F.

abstractmethod Result(

option: str,

nodeValues=True,

iter=None,

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

Returns the result. Use Results_Available() to know the available results.

Examples

For elastic simulation.

>>> simu.Results_Available()   # list valid option strings

Von Mises stress ():

>>> svm = simu.Result("Svm", nodeValues=False)  # shape (Ne,)

Node-interpolated von Mises stress:

>>> svm = simu.Result("Svm")  # shape (Nn,)

Result at a specific past iteration:

>>> u5 = simu.Result("displacement_norm", iter=5)

abstractmethod Results_Available() → list[str]

Returns a list of available results in the simulation.

Results_Get_Bc_Summary() → str

Returns the simulation loading summary.

Results_Get_Iteration_Summary() → str

Returns the iteration’s summary.

abstractmethod Results_Iter_Summary() → tuple[list[int], list[tuple[str, ndarray[tuple[Any, ...], dtype[floating]]]]]

Returns the values to be displayed in Plot_Iter_Summary.

Results_Reshape_values(

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

nodeValues: bool,

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

Reshapes input values based on whether they are stored at nodes or elements.

Parameters:

• values (_types.FloatArray) – Input values to reshape.

• nodeValues (bool) – If True, the output will represent values at nodes; if False, values on elements will be derived.

Returns:

Reshaped values on nodes or elements.

Return type:

_types.FloatArray

Raises:

Exception – Raised if unexpected conditions occur during the calculation.

Results_Set_Bc_Summary() → None

Sets the simulation loading summary.

Results_Set_Iteration_Summary(text: str, remove: bool = False) → None

Sets and print the iteration’s summary.

abstractmethod Results_dict_Energy() → dict[str, float]

Returns a dictionary containing the names and values of the calculated energies.

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

Returns displacements as a matrix [dx, dy, dz] (Nn,3).

abstractmethod Results_nodeFields_elementFields(details=False) → tuple[list[str], list[str]]

Returns lists of nodesFields and elementsFields displayed in paraview.

Save(

folder: str,

filename: str = 'simulation',

gather=False,

additionalInfos: str = '',

) → str

Saves the simulation and its summary in the folder. Saves the simulation as ‘filename.pickle’.

abstractmethod Save_Iter(iter: dict[str, Any] = {})

Saves iteration results in _results or in the specified simu.folder.

abstractmethod Set_Iter(iter: int = -1, resetAll=False) → dict

Sets the simulation to the specified iteration (the last one by default) and then reset (or not) the internal variables.

Examples

Returns the results dictionary.

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

Computes the solution field for the current boundary conditions.

Returns:

The solution of the simulation.

Return type:

_types.FloatArray

Solver_Set_Elliptic_Algorithm() → None

Sets the algorithm’s resolution properties for an elliptic problem.

Used to solve \(\Krm \, \mathrm{u} = \Frm\).

Examples

>>> simu.Solver_Set_Elliptic_Algorithm()
>>> u = simu.Solve()

Solver_Set_Hyperbolic_Algorithm(

dt: float,

algo=AlgoType.newmark,

beta=0.25,

gamma=0.5,

alpha=0.5,

) → None

Sets the algorithm’s resolution properties for a Hyperbolic problem.

Used to solve \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\).

Parameters:

• dt (float) – The time increment.

• algo (AlgoType, optional) – Algo used to solve \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\), by default AlgoType.newmark \(\Krm \, \mathrm{u}^{n+1} + \Crm \, \vrm^{n+1} + \Mrm \, \arm^{n+1} = \Frm^{n+1}\).

• beta (float, optional) – The coefficient beta, by default 1/4.

• gamma (float, optional) – The coefficient gamma, by default 1/2.

• alpha (float, optional) – The coefficient alpha, by default 1/2.

Examples

Newmark average acceleration (unconditionally stable, default):

>>> simu.Solver_Set_Hyperbolic_Algorithm(dt=0.001)  # beta=0.25, gamma=0.5
>>> simu.Solve()

HHT-α (Hilber-Hughes-Taylor, numerical damping with α=0.1):

>>> from EasyFEA import AlgoType
>>> simu.Solver_Set_Hyperbolic_Algorithm(dt=0.001, algo=AlgoType.hht, alpha=0.1)

Midpoint rule (second-order, unconditionally stable):

>>> from EasyFEA import AlgoType
>>> simu.Solver_Set_Hyperbolic_Algorithm(dt=0.001, algo=AlgoType.midpoint)

Solver_Set_Parabolic_Algorithm(dt: float, alpha=0.5) → None

Sets the algorithm’s resolution properties for a parabolic problem.

Used to solve \(\Krm \, u^{n+1} + \Crm \, \vrm^{n+1} = F^{n+1}\) with:

\(\mathrm{u}^{n+1} = \mathrm{u}^n + \dt \, \vrm^{n+\alpha}\)

Parameters:

• dt (float) – The time increment.

• alpha (float, optional) – The alpha criterion, by default 1/2n - 0 -> Forward Euler - 1 -> Backward Euler - 1/2 -> midpoint

Examples

Transient heat equation with Crank-Nicolson (α=0.5):

>>> simu.Solver_Set_Parabolic_Algorithm(dt=0.01, alpha=0.5)
>>> simu.add_dirichlet(nodes, [100.0], ["t"])
>>> for t in times:
...     simu.Solve()
...     simu.Save_Iter()

_Bc_Add_Dirichlet(

problemType: ModelType,

nodes: ndarray[tuple[Any, ...], dtype[int64]],

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

dofs: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list,

description='',

) → None

Adds Dirichlet’s boundary conditions.

If a Dirichlet’s dof is entered more than once, the conditions are added together.

_Bc_Add_Display(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list[str],

description: str,

problemType=None,

) → None

Adds a display condition.

_Bc_Add_Lagrange(newBc: LagrangeCondition)

Adds Lagrange conditions.

_Bc_Add_Neumann(

problemType: ModelType,

nodes: ndarray[tuple[Any, ...], dtype[int64]],

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

dofs: ndarray[tuple[Any, ...], dtype[int64]],

unknowns: list[str],

description='',

) → None

Adds Neumann’s boundary conditions.

If a neumann condition is already applied to the dof, the condition will not be taken into account for the dof.

_Bc_Lagrange_dim(problemType=None) → int

Calculates the dimension required to resize the system to use Lagrange multipliers.

_Check_dim_mesh_material() → None

Checks that the material dim matches the mesh dim.

_Check_dofs(problemType: ModelType, unknowns: list) → None

Checks whether the specified unknowns are available for the problem.

_Gather()

Gathers mesh partitions from all MPI ranks onto root rank.

Each mesh in the simulation’s mesh history is gathered so that rank root holds the complete global mesh. Simulation results are already fully available on every rank (allgathered by the solver), so they are left untouched.

_Get_a_n(

problemType: ModelType,

asCsrMatrix: bool = False,

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

Returns the acceleration solution associated with the given problem.

_Get_u_n(

problemType: ModelType,

asCsrMatrix: bool = False,

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

Returns the solution associated with the given problem.

_Get_v_n(

problemType: ModelType,

asCsrMatrix: bool = False,

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

Returns the speed solution associated with the given problem.

_Results_Check_Available(result: str) → bool

Check that the result is available

_Set_solutions(

problemType: ModelType,

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

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

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

) → None

_Solver_Apply_Dirichlet(

problemType: ModelType,

b: csr_matrix,

resolution: ResolType,

) → tuple[csr_matrix, csr_matrix]

Fill in the Dirichlet conditions by constructing A and x from A x = b.

Parameters:

• problemType (ModelType) – The type of problem.

• b (sparse.csr_matrix) – The b matrix.

• resolution (ResolutionType) – The resolution type.

Returns:

The A and x matrices.

Return type:

tuple[sparse.csr_matrix, sparse.csr_matrix]

_Solver_Apply_Neumann(

problemType: ModelType,

) → csr_matrix

Fill in the Neumann boundary conditions by constructing b from A x = b.

Parameters:

problemType (ModelType) – problem type

Returns:

The b vector as a csr_matrix.

Return type:

sparse.csr_matrix

_Solver_Evaluate_u_v_a_for_time_scheme(

problemType: ModelType,

u_np1: ndarray,

) → tuple[ndarray, ndarray, ndarray]

Returns u_t, v_t, and a_t vectors according to the time scheme.

Parameters:

• problemType (ModelType) – problem type

• u_np1 (np.ndarray) – the u_np1 vector

Returns:

the evaluated solutions u_t, v_t, a_t.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

_Solver_Get_Newton_Raphson_current_solution() → ndarray

Returns the current newton raphson solution.

_Solver_Get_PETSc4Py_Options(

problemType: ModelType | None = None,

) → tuple[str, str, str]

Returns (kspType, pcType, solverType) petsc4py options for the given (or default) problem type.

_Solver_Set_Newton_Raphson_Algorithm(tolConv=1e-05, maxIter=20) → None

Sets the algorithm’s resolution properties for an non linear problem.

Used to solve A(u) Δu = - R(u).

Parameters:

• tolConv (float, optional) – threshold used to check convergence, by default 1e-5

• maxIter (int, optional) – Maximum iterations for convergence, by default 20

_Solver_Set_PETSc4Py_Options(

kspType: str = 'cg',

pcType: str = 'gamg',

solverType: str = 'petsc',

problemType: ModelType | None = None,

) → None

Configure PETSc KSP solver options for the linear system Ax = b.

The constructor auto-selects defaults based on problem type and mesh dimension (see __init__). Call this method to override those defaults.

Parameters:

• kspType (str, optional) –

  Krylov subspace method, by default "cg". Choose based on matrix symmetry:

  • SPD (symmetric positive definite): "cg", "groppcg"

  • Symmetric indefinite (Lagrange/saddle-point): "minres", "symmlq"

  • General / non-symmetric: "gmres", "fgmres", "lgmres", "bcgs"

  • Direct solve (no KSP iterations): "preonly"

  • Smoother / specialised: "chebyshev" (requires spectral bounds, best as MG smoother)

  https://petsc.org/release/manualpages/KSP/KSPType/

• pcType (str, optional) –

  Preconditioner type, by default "gamg".

  pcType	Context	Notes
  "gamg" "hypre" "asm" "bjacobi" "ilu" "icc" "lu" "cholesky" "sor" "jacobi" "none"	serial + MPI serial + MPI serial + MPI serial + MPI serial only serial only see solverType see solverType serial only serial + MPI serial + MPI	PETSc AMG, no extra pkg BoomerAMG, requires hypre overlapping Schwarz block Jacobi (ILU/block) incomplete LU; use bjacobi/asm in MPI incomplete Cholesky, SPD; use bjacobi/asm in MPI direct LU factorization direct Cholesky, SPD only SOR / Gauss-Seidel diagonal scaling no preconditioning

  https://petsc.org/release/manualpages/PC/PCType/

• solverType (str, optional) –

  Backend for direct-factorization PC ("lu" / "cholesky"), by default "petsc". Ignored when pcType is not a factorization type.

  solverType	Context	Notes
  "petsc" "superlu_dist" "mumps" "mkl_pardiso" "superlu" "umfpack" "cholmod"	serial only serial + MPI serial + MPI serial + MPI serial only serial only serial only	built-in, seqaij only most common parallel direct alternative to superlu_dist Intel, fast on Intel CPUs SPD only

  https://petsc.org/release/manual/ksp/#using-external-linear-solvers

• problemType (ModelType, optional) – Problem type to configure. None applies to all problem types.

Raises:

ValueError –

• Unknown kspType, pcType, or solverType. - pcType="lu"/"cholesky" with kspType other than "preonly" (direct factorization PCs require no KSP iterations). - External solverType paired with a non-factorization pcType. - pcType="sor" in MPI context (serial-only PC). - pcType="lu"/"cholesky" with solverType="petsc" in MPI context (seqaij only — use solverType="mumps" or "superlu_dist"). - Requested external package not available in the current PETSc build.

_Solver_Solve_Newton_Raphson(

problemType=None,

tolConv=1e-05,

maxIter=20,

) → tuple[ndarray[tuple[Any, ...], dtype[floating]], int, float, list[float]]

Solves the non-linear problem using the newton raphson algorithm.

Used to solve A(u) Δu = - R(u).

Parameters:

• problemType (ModelType, optional) – The problem type, by default self.problemType

• tolConv (float, optional) – threshold used to check convergence, by default 1e-5

• maxIter (int, optional) – Maximum iterations for convergence, by default 20

Returns:

return u, Niter, timeIter, list_norm_r

Return type:

tuple[_types.FloatArray, int, float, list[float]]

Warning

The Construct_local_matrix_system function must return K and F, where K contains the tangent matrix and F contains the residual.

_Solver_Solve_problemType(

problemType: ModelType,

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

Solves the problem.

It is recommended to call the resolution via the Solve() function.

_Solver_Update_solutions(

problemType: ModelType,

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

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

Update solutions u, v and a according to x array.

Parameters:

• problemType (ModelType) – The type of problem.

• u_np1 (_types.FloatArray) – computed array in _Solver_Solve()

Returns:

returns u_np1, v_np1, a_np1

Return type:

tuple[ _types.FloatArray, Optional[_types.FloatArray], Optional[_types.FloatArray] ]

_Update(observable: Observable, event: str) → None

Receive an update/event from an observable object (observer pattern).

add_dirichlet(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds Dirichlet’s boundary conditions.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contains floats, arrays or functions or functions.

  e.g [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  The functions use the x, y and z nodes coordinates.

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Fix all displacement DOFs to zero on the left wall (elastic simulation):

>>> simu.add_dirichlet(nodes, [0, 0], ["x", "y"])

Impose temperature 100°C on the top boundary (thermal simulation):

>>> simu.add_dirichlet(nodes, [100.0], ["t"])

Spatially varying temperature (uses node coordinates):

>>> simu.add_dirichlet(nodes, [lambda x, y, z: 100 * x], ["t"])

add_lineLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds a linear load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Distributed load of -1000 N/m along y (elastic simulation):

>>> simu.add_lineLoad(nodes, [-1000], ["y"])

Spatially varying load (uses integration-point coordinates):

>>> simu.add_lineLoad(nodes, [lambda x, y, z: -500 * (1 + x)], ["y"])

add_neumann(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds Neumann’s boundary conditions.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contains floats, arrays or functions or functions.

  e.g [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  The functions use the x, y and z nodes coordinates.

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Concentrated force of -500 N in y at the beam tip (beam simulation):

>>> nodes = mesh.Nodes_Point((L, h/2))
>>> simu.add_neumann(nodes, [-500], ["y"])

add_pressureLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

magnitude: float,

problemType=None,

description='',

) → None

Adds a pressure.

Parameters:

• nodes (_types.IntArray) – nodes. (must belong to the edge of the mesh.)

• magnitude (float) – pressure magnitude

• problemType (str, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Normal inward/outward pressure of 1 MPa on a cylinder surface (elastic simulation):

>>> simu.add_pressureLoad(cylinder_nodes, magnitude=1e6)

add_surfLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds a surface load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Uniform pressure of -800 Pa along x:

>>> simu.add_surfLoad(nodes, [-800], ["x"])

add_volumeLoad(

nodes: ndarray[tuple[Any, ...], dtype[int64]],

values: list,

unknowns: list[str],

problemType=None,

description='',

) → None

Adds a volumetric load.

Parameters:

• nodes (_types.IntArray) – nodes

• values (list) –

  list of values that can contain floats, arrays or functions or lambda functions.

  e.g = [10, lambda x,y,z: 10*x - 20*y + x*z, _types.FloatArray]

  functions use x, y and z integration points coordinates (x,y,z are in this case arrays of dim (e,p))

  Please note that the functions must take 3 input parameters in the order x, y, z, whether the problem is 1D, 2D or 3D.

• unknowns (list[str]) – unknowns where values will be applied (e.g [‘y’, ‘x’])

• problemType (ModelType, optional) – problem type, if not specified, we take the basic problem of the problem

• description (str, optional) – Description of the condition, by default “”.

Examples

Gravity in -y direction, ρ = 7800 kg/m³ (elastic simulation):

>>> simu.add_volumeLoad(mesh.nodes, [-7800 * 9.81], ["y"])

property Bc_Dirichlet: list[BoundaryCondition]

Returns a copy of the Dirichlet conditions.

property Bc_Display: list[BoundaryCondition | LagrangeCondition]

Returns a copy of the boundary conditions for display.

property Bc_Lagrange: list[LagrangeCondition]

Returns a copy of the Lagrange conditions.

property Bc_Neuman: list[BoundaryCondition]

Returns a copy of the Neumann conditions.

property Niter: int

Number of iterations

property Nmesh: int

property algo: AlgoType

The algorithm used to solve the problem. (elliptic, parabolic, hyperbolic) see:

• Solver_Set_Elliptic_Algorithm() \(\Krm \, \mathrm{u} = \Frm\)

• Solver_Set_Parabolic_Algorithm() \(\Krm \, \mathrm{u} + \Crm \, \vrm = \Frm\)

• Solver_Set_Hyperbolic_Algorithm() \(\Krm \, \mathrm{u} + \Crm \, \vrm + \Mrm \, \arm = \Frm\)

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

Center of mass / barycenter / inertia center

property dim: int

simulation’s dimension

property folder: str

property inDim: int

property isGathered: bool

Whether the simulation holds the complete global mesh.

property isNonLinear

Returns whether the simulation is non linear.

property mass: float

property mesh: Mesh

simulation’s mesh.

property model: _IModel

model used

property problemType: ModelType

Get the simulation problem type.

property results: list[dict] | list[str]

Returns a copy of the list of dictionary containing the results from each iteration.

rho

mass density

property solver: str

Solver used to solve the simulation.

EasyFEA.Simulations.Load_DIC(folder: str, filename: str = 'dic') → DIC

Loads the dic analysis from the specified folder.

Parameters:

• folder (str) – The name of the folder where the dic analysis is saved.

• filename (str, optional) – The dic’s file name, by default “dic”.

Returns:

The loaded dic analysis.

Return type:

DIC

EasyFEA.Simulations.Load_Simu(

folder: str,

filename: str = 'simulation',

gather=False,

) → _Simu

Loads the simulation from the specified folder.

Parameters:

• folder (str) – simulation’s folder.

• filename (str, optional) – The simualtion’s name, by default “simulation”.

• gather (bool, optional) – Gather simu on root rank, by default False.

Returns:

The loaded simulation.

Return type:

_Simu

EasyFEA.Simulations.Load_pickle(folder: str, filename: str)

Returns folder/filename.pickle object.

EasyFEA.Simulations.Mesh_Optim_ZZ1(

DoSimu: Callable[[str], Elastic],

folder: str,

threshold: float = 0.01,

iterMax: int = 20,

coef: float = 0.5,

) → Elastic

Optimizes the mesh using ZZ1 error criterion.

Parameters:

• DoSimu (Callable[[str], Displacement]) – Function that runs a simulation and takes a .pos file as argument for mesh optimization. The function must return a Displacement simulation.

• folder (str) – Folder in which .pos files are created and then deleted.

• threshold (float, optional) – targeted error, by default 1e-2

• iterMax (int, optional) – Maximum number of iterations, by default 20

• coef (float, optional) – mesh size division ratio, by default 1/2

Returns:

Displacement simulation

Return type:

Displacement

EasyFEA.Simulations.Save_pickle(obj, folder: str, filename: str) → None

Saves the object in folder/filename.pickle.

Solvers API

Interface module to various solvers available in Python for solving linear systems (A x = b).

class EasyFEA.Simulations.Solvers.AlgoType(value)

Bases: str, Enum

Temporal integration algorithm.

All time-dependent schemes assemble and solve:

\[\Krm \, \urm^t + \Crm \, \vrm^t + \Mrm \, \arm^t = \Frm^t\]

static Get_Hyperbolic_Types() → list[str]

static Get_Hyperbolic_and_Parabolic_Types() → list[str]

elliptic = 'elliptic'

Solve \(\Krm \, \urm = \Frm\).

euler_explicit = 'euler_explicit'

Forward-Euler (first-order, conditionally stable, linear only).

Solves \(\Mrm\,\arm^n = \Frm^n - \Crm\vrm^n - \Krm\urm^n\), then updates

\(\urm^{n+1} = \urm^n + \dt\vrm^n\),

\(\vrm^{n+1} = \vrm^n + \dt\arm^n\).

Stability requires \(\dt < h_e / c\) where \(c = \sqrt{E/\rho}\).

euler_implicit = 'euler_implicit'

Backward-Euler (first-order, unconditionally stable, dissipative).

\(\urm^t = \urm^{n+1}\),

\(\vrm^t = (\urm^{n+1}-\urm^n)/\dt\),

\(\arm^t = (\vrm^t - \vrm^n)/\dt\).

hht = 'hht'

Hilber–Hughes–Taylor-\(\alpha\) method (second-order, unconditionally stable).

Generalizes newmark by shifting evaluation points by \(\alpha \in [0, 1[\), introducing controllable numerical damping (\(\alpha = 0\) recovers newmark, \(\alpha = 1/2\) recovers midpoint).

Evaluation points:

\(\urm^t = (1-\alpha)\urm^{n+1} + \alpha\urm^n\),

\(\vrm^t = (1-\alpha)\vrm^{n+1} + \alpha\vrm^n\),

\(\arm^t = (1-\alpha)\arm^{n+1} + \alpha\arm^n\).

Update: identical to newmark (default \(\beta=1/4\), \(\gamma=1/2\)).

midpoint = 'midpoint'

Midpoint rule (hht with \(\alpha=1/2\), \(\beta=1/4\), \(\gamma=1/2\)).

Update:

\(\vrm^{n+1} = 2/\dt\,(\urm^{n+1}-\urm^n) - \vrm^n\),

\(\arm^{n+1} = 2/\dt\,(\vrm^{n+1}-\vrm^n) - \arm^n\).

newmark = 'newmark'

Newmark method (default \(\beta=1/4\), \(\gamma=1/2\)).

Predictor:

\(\tilde{\urm} = \urm^n + \dt\vrm^n + \dt^2/2\,(1-2\beta)\arm^n\).

Update:

\(\arm^{n+1} = (\urm^{n+1} - \tilde{\urm})/(\beta\dt^2)\),

\(\vrm^{n+1} = \vrm^n + \dt[(1-\gamma)\arm^n + \gamma\arm^{n+1}]\).

parabolic = 'parabolic'

Solve \(\Krm \, \urm^{n+\alpha} + \Crm \, \vrm^{n+\alpha} = \Frm^{n+\alpha}\).

Set with Solver_Set_Parabolic_Algorithm().

class EasyFEA.Simulations.Solvers.ResolType(value)

Bases: str, Enum

Resolution type.

r1 = 'r1'

\(\mathrm{x}_\irm = {\Arm_{\irm\irm}}^{-1} \cdot (b_\irm - \Arm_{\irm\crm} * \mathrm{x}_\crm)\), where \(\irm\) and \(\crm\) are unknown and known degrees of freedom.

r2 = 'r2'

Lagrange multipliers

r3 = 'r3'

Penality

class EasyFEA.Simulations.Solvers.SolverType(value)

Bases: str, Enum

Solver type used to solve the linear system \(\Arm \, \mathrm{x} = \brm\)

bicg = 'bicg'

scipy.sparse.linalg.bicg

cg = 'cg'

scipy.sparse.linalg.cg

gmres = 'gmres'

scipy.sparse.linalg.gmres

lgmres = 'lgmres'

scipy.sparse.linalg.lgmres

lsq_linear = 'lsq_linear'

scipy.optimize.lsq_linear

petsc = 'petsc'

pypardiso = 'pypardiso'

pypardiso.spsolve

scipy = 'scipy'

scipy.sparse.linalg.spsolve