Understand the solve pipeline#

The solve pipeline is the internal call chain executed each time simu.Solve() is called. Every _Simu in the EasyFEA.Simulations namespace runs the same pipeline, which lives in EasyFEA/Simulations/_simu.py. This guide traces it from the moment you call Solve() to the moment the solution is stored, aimed at advanced or curious users who want to understand the internals without reading the full source.

Note

The internals described here live in EasyFEA/Simulations/_simu.py, EasyFEA/Simulations/Solvers.py, and EasyFEA/FEM/_forms.py. The single-underscore methods (e.g. _Solver_Apply_Dirichlet) are advanced API; the double-underscore ones are private and should never be called directly.

Overview#

Every call to simu.Solve() performs the same three high-level operations:

  1. Build — assemble the global sparse matrices \((K, C, M, F)\) from element integrals (skipped if nothing changed since the last solve).

  2. Apply BCs — add Neumann contributions to the right-hand side, then enforce Dirichlet constraints to form the system \(A \, x = b\).

  3. Solve — pass \(A \, x = b\) to the linear algebra backend (scipy, PETSc, pypardiso) and store the solution.

simu.Solve()
    │
    ├── 1. BUILD ──────────────────────────────────────────────────────────
    │       Get_K_C_M_F()              (cached — rebuilds only if needUpdate)
    │           └─ Assembly()
    │                └─ Construct_local_matrix_system()   ← your code
    │                       element integrals → K_e, C_e, M_e, F_e
    │                   COO scatter → global sparse K, C, M, F
    │
    ├── 2. APPLY BCs ──────────────────────────────────────────────────────
    │       _Solver_Apply_Neumann()    add flux/load contributions to F
    │       _Solver_Apply_Dirichlet()  form A from K/C/M, enforce prescribed DOFs
    │
    └── 3. SOLVE ──────────────────────────────────────────────────────────
            _Solve_Axb()              scipy / PETSc / pypardiso
            _Solver_Update_solutions() recover u, v, a from time scheme
            _Set_solutions()          store solution on the simulation object

For non-linear problems (simu.isNonLinear = True), steps 1–3 are wrapped in a Newton–Raphson loop that repeats until the residual converges — each iteration rebuilds the tangent stiffness and solves for the increment \(\Delta u\).

Step 1 — Build: matrix assembly#

Caching (Get_K_C_M_F)#

Get_K_C_M_F(problemType) returns the global matrices \((K, C, M, F)\). It maintains a dirty flag (needUpdate) and only calls Assembly() when the flag is set. The flag is raised automatically when:

  • the mesh geometry changes,

  • model material properties change,

  • a non-linear iteration updates the solution state.

When nothing has changed, the cached matrices are returned immediately — so repeated solves with the same system cost almost nothing extra.

# force a manual rebuild on the next Solve() (rarely needed)
simu.Need_Update()

Global assembly (Assembly)#

Assembly(problemType) turns element arrays into global sparse matrices:

  1. Calls Construct_local_matrix_system() to get element arrays (K_e, C_e, M_e, F_e), each of shape (Ne, nPe·dof, nPe·dof).

  2. Retrieves the global DOF index maps (rows_e, columns_e) from the element group.

  3. Scatters all element contributions into global sparse CSR matrices.

Local integration (Construct_local_matrix_system)#

This is the method you implement when subclassing _Simu. It returns the element-level arrays for one physics. When using WeakForms, EasyFEA generates it automatically from @BiLinearForm / @LinearForm decorators — but you can always implement it directly for full control.

Each BiLinearForm evaluates the weak-form integrand at all Gauss points simultaneously (vectorized over all Ne elements) — no Python loops.

Step 2 — Apply boundary conditions#

Neumann (natural) BCs#

_Solver_Apply_Neumann(problemType) adds all Neumann loads (surface tractions, line loads, point forces) to the right-hand side vector \(F\). These were registered with add_neumann, add_surfLoad, etc.

Dirichlet (essential) BCs#

_Solver_Apply_Dirichlet(problemType, b, resolution) constructs the effective system matrix \(A\) from the assembled \(K\), \(C\), \(M\) weighted by the time-integration scheme, then eliminates prescribed DOFs:

Problem type

System matrix \(A\)

Elliptic (static)

\(A = K\)

Parabolic (heat, diffusion)

\(A = K + \frac{1}{\alpha \Delta t} C\)

Newmark / Midpoint / HHT

\(A = c_K K + c_C C + c_M M\)

For non-linear problems, \(A\) is built from the tangent stiffness \(K_t\), and prescribed values are converted to incremental form \(\Delta u_D = u_D - u_D^{\text{current}}\) so that the Newton correction is consistent.

Step 3 — Solve#

Linear solve (_Solve_Axb)#

_Solve_Axb selects a resolution strategy based on the BCs present, then dispatches to the active backend:

Strategy

When

Approach

r1 (default)

Standard Dirichlet

Partition into known/unknown DOFs; solve reduced system

r2

Lagrange multiplier BCs

Augment system with Lagrange multipliers

Backend

When active

scipy (default)

Serial, no MPI

pypardiso

Serial, Intel MKL available

petsc4py

MPI (MPI_SIZE > 1) or set manually

In parallel, PETSc solves the distributed system and EasyFEA performs an Allreduce so that every rank holds the same complete DOF solution — no manual gather needed inside the solve loop.

Recover and store (_Solver_Update_solutions + _Set_solutions)#

_Solver_Update_solutions applies the active time-integration scheme to recover all solution fields from the solved \(u\):

Scheme

Recovers

Elliptic

\(u\) only

Parabolic

\(u\), \(\dot{u}\) (via backward Euler / Crank–Nicolson)

Newmark / HHT / Midpoint

\(u\), \(\dot{u}\), \(\ddot{u}\)

_Set_solutions then stores all fields on the simulation object (simu.u, simu.v, simu.a).

Non-linear problems: Newton–Raphson loop#

When isNonLinear = True, steps 1–3 are wrapped in a loop. The loop solves for an increment \(\Delta u\) at each iteration and accumulates it:

u ← Get_x0()           start from the current solution state

REPEAT:
    Need_Update()           mark K as dirty so it is rebuilt with current u
    Δu = Solve_simu(...)    solve  Kt Δu = −R(u)
    u += Δu                 accumulate the Newton correction
    check ‖R‖ < tolConv     convergence on residual norm
      OR  ‖Δu‖ < 1e-11      convergence on increment norm
UNTIL converged OR maxIter reached

Construct_local_matrix_system must return:

  • K_e — the tangent stiffness \(K_t = \partial R / \partial u\)

  • F_e — the residual \(R(u)\)

Iteration management (Save_Iter / Set_Iter)#

After each time step or load increment, call Save_Iter() to snapshot the current solution:

simu.Save_Iter()   # stores u (and v, a if present) for this step

Derived quantities (stress, strain, …) are not stored — they are recomputed on demand by Result().

To restore a previous state (e.g. for post-processing or animation):

simu.Set_Iter(i)   # restores u, v, a to the values saved at step i

Important

If the simulation uses different meshes across iterations (adaptive remeshing, crack propagation, phase-field with refinement), Set_Iter also restores the mesh that was active at step i. No manual mesh switching is needed when replaying history or generating animations.

EasyFEA beyond forward solves#

The FEM infrastructure — mesh.Get_* functions, FeArray arithmetic, Gauss-point integration — is not restricted to _Simu subclasses. You can use it directly to evaluate arbitrary integrals or construct custom operators over a mesh, without ever calling Solve().

DIC (Digital Image Correlation) is the canonical example: it is a full analysis class built on the same mesh and integration machinery, but it never solves a linear system in the traditional sense. Instead it assembles correlation operators directly from mesh.Get_* functions and minimises a correlation functional. See Digital Image Correlation (DIC) analyses for worked examples.

This means EasyFEA can serve as a general-purpose FEM toolkit for any computation that benefits from structured Gauss-point integration over a mesh, not only for forward PDE solves.