femlabpy

Python FEM library for teaching. Based on MATLAB/Scilab FemLab. Visit: https://adzetto.github.io/femlabpy

Install

pip install femlabpy

Optional:

pip install "femlabpy[mesh]"  # Gmsh 4.x support
pip install "femlabpy[gui]"   # GUI tools

Workflow Example

Step 1: Create mesh with Gmsh

import gmsh

gmsh.initialize()
gmsh.model.add("plate_with_hole")

# Create geometry
gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
gmsh.model.occ.addDisk(0.5, 0.5, 0, 0.1, 0.1)
gmsh.model.occ.cut([(2, 1)], [(2, 2)])
gmsh.model.occ.synchronize()

# Generate mesh
gmsh.model.mesh.generate(2)
gmsh.write("plate_hole.msh")
gmsh.finalize()

Step 2: Load mesh into femlabpy

import femlabpy as fp

mesh = fp.load_gmsh2("plate_hole.msh")
T = mesh.triangles  # element connectivity
X = mesh.positions  # node coordinates

Step 3: Define material properties

import numpy as np

E = 210000.0   # Young's modulus (MPa)
nu = 0.3       # Poisson's ratio
t = 1.0        # thickness (mm)

# G matrix: [E, nu, t, plane_stress_flag]
G = np.array([[E, nu, t, 1]])

Step 4: Define boundary conditions

# C matrix: [node, dof, value]
# Fix left edge (x=0): ux=0, uy=0
left_nodes = np.where(X[:, 0] < 0.01)[0] + 1
C = []
for n in left_nodes:
    C.append([n, 1, 0.0])  # ux = 0
    C.append([n, 2, 0.0])  # uy = 0
C = np.array(C)

Step 5: Define loads

# P matrix: [node, dof, value]
# Apply tension on right edge (x=1)
right_nodes = np.where(X[:, 0] > 0.99)[0] + 1
P = []
for n in right_nodes:
    P.append([n, 1, 100.0])  # Fx = 100 N
P = np.array(P)

Step 6: Solve

result = fp.elastic(T, X, G, C, P, dof=2)
u = result["u"]  # displacements
S = result["S"]  # stresses

Step 7: Visualize

fp.plotu(T, X, u, dof=2)
fp.plotelem(T, X)

Quick Examples

Cantilever beam (packaged example)

import femlabpy as fp

data = fp.canti()
result = fp.elastic(data["T"], data["X"], data["G"], data["C"], data["P"], dof=2)
print("Max displacement:", result["u"].max())
print("Max stress:", result["S"].max())

Potential flow

# Q4 mesh
result_q4 = fp.flowq4(plot=True)
print("Temperature range:", result_q4["u"].min(), "to", result_q4["u"].max())

# T3 mesh
result_t3 = fp.flowt3(plot=True)

Nonlinear truss (large deformation)

from femlabpy.examples import run_bar01_nlbar

result = run_bar01_nlbar(plot=True)
print("Load path:", result["F_path"].ravel())
print("Displacement path:", result["U_path"].ravel())

Plane stress plasticity (von Mises)

from femlabpy.examples import run_square_plastps

result = run_square_plastps(plot=True)
print("Plastic strain:", result["E"].max())

Plane strain plasticity

from femlabpy.examples import run_hole_plastpe

result = run_hole_plastpe(plot=True)

Custom FEM assembly (low-level)

import femlabpy as fp
import numpy as np

# Initialize arrays
nn = 4  # number of nodes
dof = 2  # degrees of freedom per node
K, p, q = fp.init(nn, dof)

# Element connectivity and coordinates
T = np.array([[1, 2, 3, 4, 1]])  # Q4 element
X = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
G = np.array([[210000, 0.3, 1, 1]])

# Assemble element stiffness
K = fp.kq4e(K, T, X, G)

# Apply boundary conditions
C = np.array([[1, 1, 0], [1, 2, 0], [4, 1, 0]])
P = np.array([[2, 1, 1000], [3, 1, 1000]])

p = fp.setload(p, P)
K, p, _ = fp.setbc(K, p, C, dof)

# Solve
u = np.linalg.solve(K, p)
print("Displacements:", u.ravel())

API Reference

Data Loaders

Function	Description
`canti()`	Return cantilever beam benchmark data. Returns dict with T (connectivity), X (coordinates), G (material), C (constraints), P (loads).
`flow()`	Return potential flow benchmark data for Q4 and T3 meshes.
`bar01()`	Return 2-bar truss benchmark for nonlinear analysis.
`bar02()`	Return 3-bar truss benchmark for snap-through analysis.
`bar03()`	Return 12-bar dome truss benchmark.
`square()`	Return square plate benchmark for plasticity analysis.
`hole()`	Return plate with hole benchmark for plasticity analysis.

High-Level Solvers

Function	Description
`elastic(T, X, G, C, P, dof=2, etype=None)`	Solve a linear elastic T3 or Q4 problem. The wrapper auto-detects the element family from `T` unless `etype` is given explicitly.
`flowq4(plot=False)`	Solve potential/thermal problem on Q4 mesh. Returns nodal temperatures and fluxes.
`flowt3(plot=False)`	Solve potential/thermal problem on T3 mesh. Returns nodal temperatures and fluxes.
`nlbar(T, X, G, C, P, *, no_loadsteps, i_max, i_d, plotdof, tol=1e-8, ...)`	Solve a geometrically nonlinear truss with the legacy orthogonal-residual controls.
`plastps(T, X, G, C, P, *, no_loadsteps, i_max, i_d, plotdof, ...)`	Solve the plane-stress elastoplastic legacy driver.
`plastpe(T, X, G, C, P, *, no_loadsteps, i_max, i_d, plotdof, ...)`	Solve the plane-strain elastoplastic legacy driver.

Element Stiffness (ke = element, k = assembly)

Function	Description
`ket3e(Xe, Ge)`	Compute 6x6 stiffness matrix for CST triangle (T3). Xe: 3x2 coordinates, Ge: material [E, nu, t, flag].
`kt3e(K, T, X, G)`	Assemble all T3 element stiffness matrices into global K.
`keq4e(Xe, Ge)`	Compute 8x8 stiffness matrix for bilinear quad (Q4) using 2x2 Gauss integration.
`kq4e(K, T, X, G)`	Assemble all Q4 element stiffness matrices into global K.
`kebar(Xe, A, E, u_e)`	Compute 4x4 tangent stiffness for geometrically nonlinear bar. Includes geometric stiffness.
`kbar(K, T, X, G, u)`	Assemble all bar tangent stiffness matrices into global K.
`keT4e(Xe, Ge)`	Compute 12x12 stiffness matrix for 4-node tetrahedron (T4).
`kT4e(K, T, X, G)`	Assemble all T4 element stiffness matrices into global K.
`keh8e(Xe, Ge)`	Compute 24x24 stiffness matrix for 8-node hexahedron (H8) using 2x2x2 Gauss integration.
`kh8e(K, T, X, G)`	Assemble all H8 element stiffness matrices into global K.
`ket3p(Xe, Ge)`	Compute 3x3 conductivity matrix for T3 potential element.
`kt3p(K, T, X, G)`	Assemble all T3 potential element matrices into global K.
`keq4p(Xe, Ge)`	Compute 4x4 conductivity matrix for Q4 potential element.
`kq4p(K, T, X, G)`	Assemble all Q4 potential element matrices into global K.
`keq4eps(Xe, Ge, u_e, H)`	Compute consistent tangent stiffness for plane stress plastic Q4.
`kq4eps(K, T, X, G, u, H)`	Assemble plane stress plastic Q4 tangent matrices.
`keq4epe(Xe, Ge, u_e, H)`	Compute consistent tangent stiffness for plane strain plastic Q4.
`kq4epe(K, T, X, G, u, H)`	Assemble plane strain plastic Q4 tangent matrices.

Element Response (qe = element, q = assembly)

Function	Description
`qet3e(Xe, Ge, u_e)`	Compute stress and strain for single T3 element. Returns (stress, strain).
`qt3e(q, T, X, G, u)`	Compute T3 stresses for all elements. Returns internal forces and stress/strain arrays.
`qeq4e(Xe, Ge, u_e)`	Compute stress and strain at 4 Gauss points for single Q4. Returns 4x3 arrays.
`qq4e(q, T, X, G, u)`	Compute Q4 stresses for all elements and assemble internal forces.
`qebar(Xe, A, E, u_e)`	Compute internal force for single nonlinear bar. Returns axial force and strain.
`qbar(q, T, X, G, u)`	Compute bar internal forces and assemble into global vector.
`qeT4e(Xe, Ge, u_e)`	Compute stress and strain for single T4 element.
`qT4e(q, T, X, G, u)`	Compute T4 stresses for all elements and assemble internal forces.
`qeh8e(Xe, Ge, u_e)`	Compute stress and strain at 8 Gauss points for single H8.
`qh8e(q, T, X, G, u)`	Compute H8 stresses for all elements and assemble internal forces.
`qeq4eps(Xe, Ge, u_e, H)`	Update plane stress plastic Q4 response. Returns stress, strain, updated history.
`qq4eps(q, T, X, G, u, H)`	Compute plane stress Q4 internal forces with plasticity.
`qeq4epe(Xe, Ge, u_e, H)`	Update plane strain plastic Q4 response. Returns stress, strain, updated history.
`qq4epe(q, T, X, G, u, H)`	Compute plane strain Q4 internal forces with plasticity.

Assembly and Boundary Conditions

Function	Description
`init(nn, dof, *, dynamic=False, use_sparse=None)`	Initialize the global FEM arrays. Returns `(K, p, q)` by default and `(K, M, p, q)` when `dynamic=True`.
`assmk(K, ke, nodes, dof)`	Assemble single element stiffness ke into global K at specified nodes.
`assmq(q, qe, nodes, dof)`	Assemble single element force qe into global internal force vector q.
`setload(p, P)`	Set nodal loads from the FemLab load table. Replaces existing values at the referenced global DOFs.
`addload(p, P)`	Add nodal loads from the FemLab load table. Accumulates with existing values.
`setbc(K, p, C, dof)`	Apply Dirichlet BCs using direct elimination. Zeros rows/columns, sets diagonal.
`solve_lag(K, p, C, dof)`	Solve with Lagrange multipliers for Dirichlet constraints.
`solve_lag_general(K, p, G, Q)`	Solve with general linear constraints Gu = Q.
`reaction(q, C, dof, comp=None)`	Extract support reactions from the assembled internal-force vector at constrained DOFs.
`rnorm(r, C, dof)`	Compute residual norm excluding constrained DOFs.

Material Models

Function	Description
`devstress(S)`	Compute deviatoric stress and mean stress. S: [sxx, syy, sxy] or [sxx, syy, szz, sxy, syz, sxz].
`eqstress(S)`	Compute von Mises equivalent stress from stress vector.
`yieldvm(S, mat, ep, dL)`	Evaluate von Mises yield function f = seq - (sy + H*ep). mat: [E, nu, sy, H].
`dyieldvm(S, mat, ep, dL)`	Derivative of yield function with respect to plastic multiplier dL.
`stressvm(S_trial, mat, ep)`	Perform von Mises radial return mapping. Returns (S_updated, delta_ep).
`stressdp(S_trial, mat, ep)`	Perform Drucker-Prager return mapping with Newton iterations.

Mesh I/O

Function	Description
`load_gmsh(filename)`	Load Gmsh mesh file (.msh). Returns GmshMesh with positions, triangles, quads, etc. Legacy 2.x format.
`load_gmsh2(filename)`	Load Gmsh mesh with flexible format detection. Supports both 2.x and 4.x (requires [mesh] extra).

GmshMesh attributes:

positions: Nx3 node coordinates
triangles: Mx4 triangle connectivity [elem_id, n1, n2, n3]
quads: Mx5 quad connectivity
lines: Mx3 line connectivity
bounds_min, bounds_max: Bounding box

Plotting

Function	Description
`plotelem(T, X, line_style="k-", nonum=False, noelem=False)`	Plot the undeformed mesh and optionally annotate node or element ids.
`plotforces(T, X, P)`	Plot load arrows on the mesh.
`plotbc(T, X, C)`	Plot boundary-condition markers.
`plotu(T, X, u, dof=None, component=0)`	Plot a scalar nodal field or a contour extracted from a flattened displacement vector. `component=0` plots magnitude, `1` plots x, `2` plots y.
`plotq4(T, X, S, scomp)`	Plot a one-based Q4 stress or strain component reconstructed from Gauss-point data.
`plott3(T, X, S, scomp)`	Plot a one-based T3 stress or strain component.

Examples Module

Function	Description
`run_cantilever(plot=False)`	Solve cantilever beam. Returns u, S, E, R.
`run_flow_q4(plot=False)`	Solve Q4 potential flow. Returns u (temperatures).
`run_flow_t3(plot=False)`	Solve T3 potential flow. Returns u (temperatures).
`run_bar01_nlbar(plot=False)`	Solve 2-bar snap-through. Returns load-displacement path.
`run_bar02_nlbar(plot=False)`	Solve 3-bar snap-through. Returns load-displacement path.
`run_bar03_nlbar(plot=False)`	Solve 12-bar dome. Returns load-displacement path.
`run_square_plastps(plot=False)`	Solve square plate plane stress plasticity.
`run_square_plastpe(plot=False)`	Solve square plate plane strain plasticity.
`run_hole_plastps(plot=False)`	Solve plate with hole plane stress plasticity.
`run_hole_plastpe(plot=False)`	Solve plate with hole plane strain plasticity.
`run_ex_lag_mult(plot=False)`	Solve 3-bar truss with displacement constraint.
`run_gmsh_triangle(plot=False)`	Solve imported Gmsh triangle mesh.

Data loaders: cantilever_data(), flow_data(), bar01_data(), bar02_data(), bar03_data(), square_data(plane_strain=False), hole_data(plane_strain=False), ex_lag_mult_data(), gmsh_triangle_data()

Dynamic Analysis

Function	Description
`solve_newmark(M, C, K, p_func, u0, v0, dt, nsteps, ...)`	Solve time history using implicit Newmark-beta method. Returns `TimeHistory` object.
`solve_central_diff(M_lump, C, K, p_func, u0, v0, dt, ...)`	Solve time history using explicit central difference method. Requires lumped mass.
`solve_hht(M, C, K, p_func, u0, v0, dt, nsteps, ...)`	Solve time history using HHT-alpha method for high-frequency numerical damping.
`solve_modal(K, M, n_modes, C_bc, dof)`	Compute natural frequencies and mode shapes. Returns `ModalResult` object.
`seismic_load(M, direction, accel_record, dt_record)`	Create an effective seismic load function $p(t) = -M a_g(t)$ from ground acceleration.
`compute_frf(M, C, K, input_dof, output_dof, freq_range)`	Compute Frequency Response Function (FRF) for harmonic excitation.
`rayleigh_damping(M, K, f1, f2, zeta1, zeta2)`	Compute Rayleigh damping matrix $C = \alpha M + \beta K$ from two modal frequencies.
`plot_time_history(result, dof_index)`	Plot displacement/velocity/acceleration vs time.
`plot_modes(X, result, scale)`	Plot deformed mode shapes.

Periodic Boundaries

Function	Description
`find_periodic_pairs(X, n1_nodes, n2_nodes, tol)`	Find matching node pairs on opposite boundaries for periodic conditions.
`periodic_constraints(pairs, dof)`	Generate linear constraint equations (G matrix) for periodic boundary conditions.
`homogenize(T, X, G, pairs, dof)`	Compute the 3x3 effective homogenized stiffness tensor (C_eff) for a unit cell.

Dynamic Analysis & Benchmarks

Dynamic analysis examples and benchmarks against OpenSeesPy and CalculiX.

2D plane-stress cantilever beam (32x4 Q4 elements). Tested for static tip deflection and natural frequencies.

Metric	Analytical (E-B)	OpenSeesPy	CalculiX (CPS4)	femlabpy (Q4)	femlabpy vs OpenSees
Static Deflection	-9.752e-05 m	-9.544e-05 m	-9.518e-05 m	-9.544e-05 m	0.000% (Exact)
Frequency: Mode 1	26.110 Hz	26.228 Hz	26.283 Hz	26.248 Hz	0.07%
Frequency: Mode 2	163.628 Hz	153.716 Hz	154.668 Hz	154.437 Hz	0.47%
Frequency: Mode 3	458.163 Hz	323.810 Hz	323.892 Hz	323.881 Hz	0.02%

Note: Small differences against 1D Euler-Bernoulli theory are expected for 2D plane-stress elements.

2. Seismic Time History Analysis

2D concrete column (4x32 Q4 elements) under horizontal earthquake (Düzce, BOL090.AT2, PGA = 0.82g). Solved with Newmark-beta average acceleration and 5% Rayleigh damping for 5590 steps (dt = 0.01s).

Metric	OpenSeesPy (UniformExcitation)	CalculiX (*DLOAD GRAV)	femlabpy (seismic_load)
Max Roof Disp. (X)	5.850 mm	4.180 mm	5.848 mm
Solve Time	1.33 seconds	141.37 seconds	1.12 seconds

Conclusion: femlabpy matches OpenSeesPy with a 0.67% RMS difference over 5590 steps. femlabpy solve time is 1.12 seconds.

Development

pytest -q                    # run tests
python -m build              # build package

References

Bathe, K.J. (2014). Finite Element Procedures. Prentice Hall.
Zienkiewicz, Taylor, Zhu (2013). The Finite Element Method. Elsevier.
Hughes, T.J.R. (2000). The Finite Element Method. Dover.
Cook et al. (2002). Concepts and Applications of FEA. Wiley.
de Souza Neto et al. (2008). Computational Methods for Plasticity. Wiley.
Gmsh: https://gmsh.info/

Metadata

Release files for femlabpy 0.6.2

For a detailed explanation of source distributions (sdists) and built distributions (wheels), please see the package formats documentation.

Source distribution (sdist)

Source distribution for femlabpy 0.6.2
File	Size	Uploaded
femlabpy-0.6.2.tar.gz	148.6 kB	Apr 13, 2026	Details

Built distribution (wheel)

Table of built distributions (wheels) for femlabpy 0.6.2
File	Interpreter	ABI	Platform	Reset
femlabpy-0.6.2-py3-none-any.whl 154.0 kB Apr 13, 2026	Python 3	none	any	Details

Total release size: 302.7 kB

Release files / femlabpy-0.6.2.tar.gz

Download URL	femlabpy-0.6.2.tar.gz
Size	148.6 kB
Tags	Source
SHA-256 checksum How to use checksums	`155c9762b850618614f968d79aa03cc07de575bab11d33849db0d45cbe526514`
BLAKE2b-256 checksum How to use checksums	`e013f411656172b07329a779796f0f8e958ac598e2edc622d304dabb7f9f79f8`
Upload date	Apr 13, 2026
Uploaded using Trusted Publishing? What is trusted publishing?	No
Uploaded via	`twine/6.2.0 CPython/3.13.9`

Release files / femlabpy-0.6.2-py3-none-any.whl

Download URL	femlabpy-0.6.2-py3-none-any.whl
Size	154.0 kB
Tags	Python 3
SHA-256 checksum How to use checksums	`c931870537083f171b5ca97466ae1bba32143ea67f8ff7b654596dcb3c5b4a3c`
BLAKE2b-256 checksum How to use checksums	`3b7dc7643394e885070caef3512bfc368a15ccd46e3bfdfb6b7a90eed6d74ccb`
Upload date	Apr 13, 2026
Uploaded using Trusted Publishing? What is trusted publishing?	No
Uploaded via	`twine/6.2.0 CPython/3.13.9`