Abstract

1. Introduction

Structural simulation is a cornerstone of modern civil, mechanical, and aerospace engineering. The ability to predict deformation, stress, and failure in complex structures before physical prototyping reduces cost, accelerates design cycles, and improves safety. Historically, commercial software such as ANSYS, ABAQUS, and NASTRAN dominated this space. However, the rapid maturation of the Python ecosystem has democratized access to powerful simulation tools, offering transparency, extensibility, and integration with broader scientific workflows.

Python's appeal for structural simulation stems from several factors:

1. Rich numerical ecosystem: NumPy, SciPy, and sparse matrix libraries provide high-performance linear algebra.
2. Visualization: Matplotlib, PyVista, and Paraview Python bindings enable rich post-processing.
3. Machine learning integration: Seamless interoperability with TensorFlow, PyTorch, and scikit-learn opens pathways to physics-informed neural networks (PINNs).
4. Open-source transparency: Researchers can inspect, modify, and extend solver internals.
5. Reproducibility: Python scripts serve as executable documentation for simulation workflows.

This survey targets researchers and engineers seeking to select the appropriate Python library for structural simulation tasks ranging from simple beam analysis to nonlinear transient dynamics of large-scale structures.

1.1 Scope

This review covers:

• Linear and nonlinear static analysis
• Modal and transient dynamic analysis
• Geometric and material nonlinearity
• Multiphysics coupling relevant to structural problems
• Mesh generation and post-processing

2. Background: Fundamentals of Structural Simulation

2.1 Governing Equations

Structural simulation is rooted in the equilibrium of internal and external forces. For a continuum body occupying domain Ω, the governing partial differential equation (strong form) is:

∇·σ + b = ρü in Ω

where:

• σ = Cauchy stress tensor
• b = body force vector
• ρ = material density
• ü = acceleration vector

Boundary conditions:

• Essential (Dirichlet): u = ū on Γ_u
• Natural (Neumann): σ·n = t̄ on Γ_t

2.2 Finite Element Discretization

The Finite Element Method (FEM) approximates the displacement field using shape functions N_i:

u(x) ≈ Σ N_i(x) · u_i

This yields the discrete system:

[K]{u} = {f} (static)
[M]{ü} + [C]{u̇} + [K]{u} = {f(t)} (dynamic)

where K, M, C are stiffness, mass, and damping matrices respectively.

3. Library Survey

3.1 OpenSeesPy

Repository: https://github.com/zhuminjie/OpenSeesPy
Documentation: https://openseespydoc.readthedocs.io
License: BSD
Primary Domain: Earthquake engineering, geomechanics

OpenSeesPy is the Python interpreter interface to OpenSees (Open System for Earthquake Engineering Simulation), developed at UC Berkeley. It provides access to a vast library of element formulations, material models, and nonlinear solvers specifically calibrated for structural and geotechnical earthquake engineering.

Key Features:

• Nonlinear fiber-section beam-column elements
• Nonlinear material models (Steel02, Concrete01, Hysteretic)
• Integration algorithms: Newmark, HHT, GeneralizedAlpha
• Parallel execution via OpenMP
• Direct access to element state variables

Strengths:

• Industry-validated for seismic analysis
• Rich material library for concrete and steel
• Direct equivalence to peer-reviewed OpenSees models

Limitations:

• Primarily 1D/2D element focus (beams, columns, frames)
• Less suited for detailed 3D continuum problems
• Steeper learning curve for non-earthquake engineers

Sample Code 3.1 — Nonlinear Static Pushover of a Steel Frame

import openseespy.opensees as opsimport numpy as npimport matplotlib.pyplot as plt# ─────────────────────────────────────────# MODEL SETUP: 2D Steel Moment Frame# Units: kN, m# ─────────────────────────────────────────ops.wipe()ops.model('basic', '-ndm', 2, '-ndf', 3)# --- Nodes ---# Ground levelops.node(1, 0.0, 0.0)ops.node(2, 6.0, 0.0)# First floorops.node(3, 0.0, 4.0)ops.node(4, 6.0, 4.0)# Second floorops.node(5, 0.0, 8.0)ops.node(6, 6.0, 8.0)# --- Boundary conditions (fixed base) ---ops.fix(1, 1, 1, 1)ops.fix(2, 1, 1, 1)# --- Material: Steel02 (Giuffre-Menegotto-Pinto) ---matTag = 1Fy = 250e3    # kPa (yield strength)E  = 200e6    # kPa (elastic modulus)b  = 0.01     # strain hardening ratioops.uniaxialMaterial('Steel02', matTag, Fy, E, b, 20.0, 0.925, 0.15)# --- Cross-section: W-section via fiber discretization ---secTag = 1H  = 0.300   # m (section depth)Bf = 0.150   # m (flange width)tf = 0.012   # m (flange thickness)tw = 0.008   # m (web thickness)ops.section('Fiber', secTag)# Top flange (5 fibers)for i in range(5):    y = H/2 - tf/2    z = -Bf/2 + (i + 0.5) * Bf/5    A = Bf/5 * tf    ops.fiber(y, z, A, matTag)# Bottom flange (5 fibers)for i in range(5):    y = -H/2 + tf/2    z = -Bf/2 + (i + 0.5) * Bf/5    A = Bf/5 * tf    ops.fiber(y, z, A, matTag)# Web (10 fibers)for i in range(10):    y = -H/2 + tf + (i + 0.5) * (H - 2*tf) / 10    z = 0.0    A = tw * (H - 2*tf) / 10    ops.fiber(y, z, A, matTag)# --- Geometric transformation ---ops.geomTransf('PDelta', 1)  # P-Delta formulation# --- Elements: Nonlinear Beam-Column ---numIntPts = 5beamTag = 1# Columns (floor 1)ops.element('nonlinearBeamColumn', 1, 1, 3, numIntPts, secTag, beamTag)ops.element('nonlinearBeamColumn', 2, 2, 4, numIntPts, secTag, beamTag)# Beams (floor 1)ops.element('nonlinearBeamColumn', 3, 3, 4, numIntPts, secTag, beamTag)# Columns (floor 2)ops.element('nonlinearBeamColumn', 4, 3, 5, numIntPts, secTag, beamTag)ops.element('nonlinearBeamColumn', 5, 4, 6, numIntPts, secTag, beamTag)# Beams (floor 2)ops.element('nonlinearBeamColumn', 6, 5, 6, numIntPts, secTag, beamTag)# ─────────────────────────────────────────# LOADING: Gravity first, then lateral pushover# ─────────────────────────────────────────# Gravity loadsops.timeSeries('Constant', 1)ops.pattern('Plain', 1, 1)ops.eleLoad('-ele', 3, '-type', '-beamUniform', -30.0)  # kN/mops.eleLoad('-ele', 6, '-type', '-beamUniform', -30.0)ops.constraints('Plain')ops.numberer('RCM')ops.system('BandGeneral')ops.test('NormUnbalance', 1e-8, 25)ops.algorithm('Newton')ops.integrator('LoadControl', 0.1)ops.analysis('Static')ops.analyze(10)ops.loadConst('-time', 0.0)# Lateral pushover pattern (inverted triangular)ops.timeSeries('Linear', 2)ops.pattern('Plain', 2, 2)ops.load(3, 1.0, 0.0, 0.0)   # kN at floor 1ops.load(5, 2.0, 0.0, 0.0)   # kN at floor 2 (2× for inverted triangle)# ─────────────────────────────────────────# PUSHOVER ANALYSIS# ─────────────────────────────────────────dU      = 0.001    # displacement increment (m)maxU    = 0.15     # maximum roof displacement (m)nSteps  = int(maxU / dU)ops.integrator('DisplacementControl', 5, 1, dU)ops.analysis('Static')disp_history = [0.0]base_shear   = [0.0]for i in range(nSteps):    ok = ops.analyze(1)    if ok != 0:        print(f"Analysis failed at step {i}")        break    d = ops.nodeDisp(5, 1)    V = -(ops.eleResponse(1, 'localForce')[0] + ops.eleResponse(2, 'localForce')[0])    disp_history.append(d)    base_shear.append(V)# ─────────────────────────────────────────# POST-PROCESSING# ─────────────────────────────────────────plt.figure(figsize=(8, 5))plt.plot(disp_history, base_shear, 'b-', linewidth=2)plt.xlabel('Roof Displacement (m)')plt.ylabel('Base Shear (kN)')plt.title('Nonlinear Pushover Curve — 2-Story Steel Frame')plt.grid(True, linestyle='--', alpha=0.6)plt.savefig('pushover_curve.png', dpi=150, bbox_inches='tight')plt.show()print("Peak base shear:", max(base_shear), "kN")

3.2 FEniCSx (DOLFINx)

Repository: https://github.com/FEniCS/dolfinx
Documentation: https://docs.fenicsproject.org
License: LGPL-3.0
Primary Domain: Continuum mechanics, multiphysics, PDE solvers

FEniCSx (previously FEniCS/DOLFIN) is a high-performance computing platform for solving PDEs via the finite element method. It employs the Unified Form Language (UFL) to express variational formulations symbolically, enabling automatic generation of low-level assembly code.

Key Features:

• Symbolic variational form specification via UFL
• MPI-parallel distributed computing
• Arbitrary-order polynomial basis functions
• Mixed-element formulations
• Automatic differentiation for sensitivity analysis
• Built-in mesh generation and refinement

Strengths:

• Research-grade flexibility for custom constitutive models
• Excellent for 3D continuum problems
• Publication-quality scientific outputs
• Active development with growing Python API

Limitations:

• Installation can be complex (typically Docker-recommended)
• Limited built-in structural elements (no beam/shell by default)
• Steeper learning curve for non-FEM experts

Sample Code 3.2 — Linear Elasticity: Cantilever Beam Under Tip Load

import numpy as npfrom mpi4py import MPIfrom dolfinx import mesh, fem, default_scalar_typefrom dolfinx.fem.petsc import LinearProblemimport ufl# ─────────────────────────────────────────# GEOMETRY: Cantilever beam# Dimensions: L=1m, H=0.1m, W=0.1m# ─────────────────────────────────────────L, H, W = 1.0, 0.1, 0.1nx, ny, nz = 20, 4, 4domain = mesh.create_box(    MPI.COMM_WORLD,    [[0.0, 0.0, 0.0], [L, H, W]],    [nx, ny, nz],    mesh.CellType.hexahedron)# ─────────────────────────────────────────# FUNCTION SPACE: Lagrange P1 vector elements# ─────────────────────────────────────────V = fem.functionspace(domain, ("Lagrange", 1, (3,)))# ─────────────────────────────────────────# MATERIAL: Isotropic linear elastic# E = 200 GPa, nu = 0.3 (structural steel)# ─────────────────────────────────────────E_mod = 200e9   # Panu    = 0.3mu    = fem.Constant(domain, E_mod / (2 * (1 + nu)))lam   = fem.Constant(domain, E_mod * nu / ((1 + nu) * (1 - 2 * nu)))def epsilon(u):    """Small strain tensor."""    return ufl.sym(ufl.grad(u))def sigma(u):    """Cauchy stress tensor — isotropic linear elasticity."""    return lam * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)# ─────────────────────────────────────────# BOUNDARY CONDITIONS# ─────────────────────────────────────────def clamped_boundary(x):    """Left face: x ≈ 0"""    return np.isclose(x[0], 0)# Fix all DOFs on clamped faceclamped_dofs = fem.locate_dofs_geometrical(V, clamped_boundary)u_bc = np.array([0, 0, 0], dtype=default_scalar_type)bc   = fem.dirichletbc(u_bc, clamped_dofs, V)# ─────────────────────────────────────────# VARIATIONAL FORMULATION# ─────────────────────────────────────────u = ufl.TrialFunction(V)v = ufl.TestFunction(V)# Body force: gravity in -y directionf = fem.Constant(domain, default_scalar_type((0.0, -7850 * 9.81, 0.0)))  # steel density# Traction: tip load on right faceT = fem.Constant(domain, default_scalar_type((0.0, -1e5, 0.0)))  # 100 kN/m² downward# Locate right face for Neumann BCdef right_face(x):    return np.isclose(x[0], L)facet_indices, facet_markers = [], []fdim = domain.topology.dim - 1right_facets = mesh.locate_entities_boundary(domain, fdim, right_face)facet_indices.extend(right_facets)facet_markers.extend([1] * len(right_facets))facet_tag = mesh.meshtags(domain, fdim,                            np.array(facet_indices, dtype=np.int32),                            np.array(facet_markers, dtype=np.int32))ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)# Bilinear and linear formsa = ufl.inner(sigma(u), epsilon(v)) * ufl.dxL_form = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(1)# ─────────────────────────────────────────# SOLVE# ─────────────────────────────────────────problem = LinearProblem(a, L_form, bcs=[bc],                        petsc_options={"ksp_type": "preonly",                                       "pc_type":  "lu"})uh = problem.solve()uh.name = "Displacement"# ─────────────────────────────────────────# POST-PROCESSING# ─────────────────────────────────────────# Compute von Mises stresss     = sigma(uh) - (1.0/3) * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))von_mises = ufl.sqrt((3.0/2) * ufl.inner(s, s))W_stress = fem.functionspace(domain, ("DG", 0))stress_expr = fem.Expression(von_mises, W_stress.element.interpolation_points())vm_stress = fem.Function(W_stress)vm_stress.interpolate(stress_expr)vm_stress.name = "VonMises"# Export to XDMF for Paraview visualizationfrom dolfinx.io import XDMFFilewith XDMFFile(MPI.COMM_WORLD, "cantilever_results.xdmf", "w") as xdmf:    xdmf.write_mesh(domain)    xdmf.write_function(uh)    xdmf.write_function(vm_stress)# Tip displacementtip_node = np.array([L, H/2, W/2])from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cellstree    = bb_tree(domain, domain.topology.dim)cells   = compute_collisions_points(tree, tip_node.reshape(1,3))collide = compute_colliding_cells(domain, cells, tip_node.reshape(1,3))tip_disp = uh.eval(tip_node.reshape(1,3), collide.links(0))[0]print(f"Tip displacement (y): {tip_disp[1]*1000:.4f} mm")print(f"Theoretical (EulerBernoulli): {-1e5*H*W*L**3/(3*200e9*(H**3*W/12))*1000:.4f} mm")

3.3 PyNite (Structural Frame Analysis)

Repository: https://github.com/JWock82/PyNite
License: MIT
Primary Domain: 3D frame analysis, beams, trusses, plates

PyNite is a Python library purpose-built for structural engineers requiring 3D frame and finite element analysis without the complexity of research-grade solvers. It implements classical beam theory elements, spring elements, and plate bending elements in an accessible, Pythonic API.

Key Features:

• 3D beam-column elements (Euler-Bernoulli / Timoshenko)
• Plate/shell elements (MITC4 formulation)
• Spring supports (linear and nonlinear)
• Load combinations and enveloping
• Built-in rendering via Matplotlib and PyVista

Strengths:

• Extremely accessible for structural engineers
• Clear and readable API
• Suitable for everyday structural design checks
• No complex installation requirements

Limitations:

• Linear analysis only (no material nonlinearity)
• Not suitable for large-scale continuum problems
• Limited element library

Sample Code 3.3 — 3D Portal Frame Analysis

from PyNite import FEModel3Dimport matplotlib.pyplot as plt# ─────────────────────────────────────────# BUILD MODEL: 2-bay, 1-story portal frame# Units: kips, inches# ─────────────────────────────────────────frame = FEModel3D()# --- Material: Structural Steel A992 ---frame.add_material('Steel', E=29000, G=11200, nu=0.3, rho=0.000284)# --- Section properties: W12x96 ---# Ix, Iy, J for torsion, A cross-sectionframe.add_section('W12x96',                  A=28.2,    # in²                  Iy=270.0,  # in⁴ (weak axis)                  Iz=833.0,  # in⁴ (strong axis)                  J=2.82,    # in⁴ (torsional constant)                  material='Steel')# --- Nodes ---# Column bases (z=0)frame.add_node('N1', 0,   0, 0)frame.add_node('N2', 240, 0, 0)  # 20 ft bayframe.add_node('N3', 480, 0, 0)  # 20 ft bay# Column tops (z=144 = 12 ft)frame.add_node('N4', 0,   0, 144)frame.add_node('N5', 240, 0, 144)frame.add_node('N6', 480, 0, 144)# --- Supports: Fixed bases ---for n in ['N1', 'N2', 'N3']:    frame.def_support(n, True, True, True, True, True, True)# --- Members ---# Columnsframe.add_member('Col1', 'N1', 'N4', 'Steel', 'W12x96')frame.add_member('Col2', 'N2', 'N5', 'Steel', 'W12x96')frame.add_member('Col3', 'N3', 'N6', 'Steel', 'W12x96')# Beamsframe.add_member('Beam1', 'N4', 'N5', 'Steel', 'W12x96')frame.add_member('Beam2', 'N5', 'N6', 'Steel', 'W12x96')# ─────────────────────────────────────────# LOADS# ─────────────────────────────────────────# Load case 'D' — Dead loadframe.add_member_dist_load('Beam1', 'Fy', -0.8, -0.8, 'D')  # kip/in (gravity)frame.add_member_dist_load('Beam2', 'Fy', -0.8, -0.8, 'D')# Load case 'W' — Wind (lateral point load at roof)frame.add_node_load('N4', 'FX', 15.0, 'W')  # kipsframe.add_node_load('N6', 'FX', 15.0, 'W')# Load combinations (LRFD)frame.add_load_combo('1.2D+1.6W',  {'D': 1.2, 'W': 1.6})frame.add_load_combo('1.2D+1.0W',  {'D': 1.2, 'W': 1.0})# ─────────────────────────────────────────# ANALYZE# ─────────────────────────────────────────frame.analyze()# ─────────────────────────────────────────# RESULTS# ─────────────────────────────────────────print("n=== Portal Frame Analysis Results ===")print(f"{'Node':<8} {'DX (in)':<12} {'DY (in)':<12} {'DZ (in)':<12}")for node_id in ['N4', 'N5', 'N6']:    n = frame.nodes[node_id]    dx = n.DX['1.2D+1.6W']    dy = n.DY['1.2D+1.6W']    dz = n.DZ['1.2D+1.6W']    print(f"{node_id:<8} {dx:<12.4f} {dy:<12.4f} {dz:<12.4f}")# Member end forcesbeam1 = frame.members['Beam1']print("n=== Beam1 End Forces (1.2D+1.6W) ===")print(f"  Start Shear Fy : {beam1.shear('Fy', 0, '1.2D+1.6W'):.2f} kips")print(f"  Mid   Moment Mz: {beam1.moment('Mz', 120, '1.2D+1.6W'):.2f} kip·in")print(f"  End   Shear Fy : {beam1.shear('Fy', 240, '1.2D+1.6W'):.2f} kips")# Render modelframe.render_model(render_loads=True, combo_name='1.2D+1.6W')

3.4 SfePy (Simple Finite Elements in Python)

Repository: https://github.com/sfepy/sfepy
Documentation: https://sfepy.org
License: BSD
Primary Domain: Continuum mechanics, acoustics, multiphysics

SfePy is a general-purpose FEM solver written entirely in Python with performance-critical routines in C/Cython. It supports a wide range of physics and is particularly strong for custom PDE formulations.

Key Features:

• Supports 1D/2D/3D problems
• Rich element library (Lagrange, Hermite, Raviart-Thomas)
• Homogenization framework for multiscale problems
• Piezoelectric, acoustic, fluid-structure coupling
• Problem definition via Python dictionaries or input files

Sample Code 3.4 — Plane Stress Analysis of a Perforated Plate

# sfepy_plate_with_hole.py# Plane stress analysis: square plate with central circular hole# Loaded in tension — stress concentration factor verificationimport numpy as npfrom sfepy.discrete import FieldVariable, Material, Integral, Equation, Equations, Problemfrom sfepy.discrete.fem import FEDomain, Fieldfrom sfepy.mesh.mesh_generators import gen_block_meshfrom sfepy.terms import Termfrom sfepy.solvers.ls import ScipyDirectfrom sfepy.solvers.nls import Newton# ─────────────────────────────────────────# MESH: rectangular plate (will need hole separately or via external mesh)# Here we demonstrate the SfePy API structure# ─────────────────────────────────────────# In practice, load a mesh with hole cut (e.g., from Gmsh):# from sfepy.mesh.mesh_io import MeshIO# mesh = MeshIO.any_from_filename('plate_with_hole.vtk').read()# For demonstration: simple block meshmesh = gen_block_mesh(    dims=[1.0, 1.0],    shape=[20, 20],    centre=[0.5, 0.5])domain = FEDomain('domain', mesh)omega  = domain.create_region('Omega', 'all')left   = domain.create_region('Left',  'vertices in (x < 0.001)', 'facet')right  = domain.create_region('Right', 'vertices in (x > 0.999)', 'facet')# ─────────────────────────────────────────# FIELD AND VARIABLES# ─────────────────────────────────────────field = Field.from_args('vector', np.float64, 'vector', omega, approx_order=1)u = FieldVariable('u', 'unknown', field)v = FieldVariable('v', 'test',    field, primary_var_name='u')# ─────────────────────────────────────────# MATERIAL: Steel, plane stress# ─────────────────────────────────────────m = Material('m', D=None)  # Will be set via constitutive hookE_val, nu_val = 200e9, 0.3lam_val = E_val * nu_val / ((1 + nu_val) * (1 - 2*nu_val))mu_val  = E_val / (2 * (1 + nu_val))def get_D(ts, coors, mode=None, **kwargs):    """Plane stress elasticity tensor (Voigt notation)."""    if mode == 'qp':        D = E_val / (1 - nu_val**2) * np.array([            [1,     nu_val, 0              ],            [nu_val, 1,    0              ],            [0,     0,     (1-nu_val)/2   ]        ])        nqp = coors.shape[0]        return {'D': np.tile(D, (nqp, 1, 1))}m = Material('m', function=get_D)# Traction on right facetraction = Material('traction', val=[[100e6], [0.0]])  # 100 MPa in x# ─────────────────────────────────────────# EQUATIONS# ─────────────────────────────────────────integral = Integral('i', order=2)t1 = Term.new('dw_lin_elastic(m.D, v, u)',              integral, omega, m=m, v=v, u=u)t2 = Term.new('dw_surface_ltr(traction.val, v)',              integral, right, traction=traction, v=v)eq = Equation('balance', t1 - t2)eqs = Equations([eq])# ─────────────────────────────────────────# BOUNDARY CONDITIONS# ─────────────────────────────────────────from sfepy.discrete.conditions import Conditions, EssentialBCbc_left = EssentialBC('bc_left', left, {'u.0': 0.0})   # no x-displacement on leftbc_sym  = EssentialBC('bc_sym',  left, {'u.1': 0.0})   # symmetry# ─────────────────────────────────────────# SOLVER SETUP# ─────────────────────────────────────────ls  = ScipyDirect({})nls = Newton({}, lin_solver=ls)pb = Problem('elastic', equations=eqs)pb.set_bcs(ebcs=Conditions([bc_left, bc_sym]))pb.set_solver(nls)status = pb.solve()# Extract displacementsu_val = pb.get_variables()['u']print("Max x-displacement:", u_val().reshape(-1, 2)[:, 0].max(), "m")print("Max y-displacement:", u_val().reshape(-1, 2)[:, 1].max(), "m")

3.5 Anastruct (2D Frame & Truss)

Repository: https://github.com/ritchie46/anaStruct
License: MIT
Primary Domain: 2D structural analysis, teaching, rapid prototyping

Anastruct provides a clean, Pythonic interface for 2D frame and truss analysis with an emphasis on visualization and ease of use for structural engineering education.

Sample Code 3.5 — Continuous Beam Analysis with Influence Lines

from anastruct import SystemElementsimport matplotlib.pyplot as plt# ─────────────────────────────────────────# 3-SPAN CONTINUOUS BEAM# Span: 5m + 7m + 5m = 17m total# ─────────────────────────────────────────ss = SystemElements()# E = 200 GPa, I = 3.33e-4 m⁴ (IPE400)EI = 200e9 * 3.33e-4  # N·m²# Add beam elements (0.5m segments for accuracy)nodes = [0, 5, 12, 17]x = 0for i in range(len(nodes)-1):    span   = nodes[i+1] - nodes[i]    nsegs  = int(span / 0.25)    dx     = span / nsegs    for _ in range(nsegs):        ss.add_element(            location=[[x, 0], [x + dx, 0]],            EI=EI,            EA=200e9 * 0.01  # A = 100 cm²        )        x += dx# --- Supports ---ss.add_support_hinged(node_id=ss.find_node_id(vertex=[0,   0]))   # leftss.add_support_roll(  node_id=ss.find_node_id(vertex=[5,   0]))   # mid-leftss.add_support_roll(  node_id=ss.find_node_id(vertex=[12,  0]))   # mid-rightss.add_support_hinged(node_id=ss.find_node_id(vertex=[17,  0]))   # right# --- Loads: distributed on all spans ---ss.q_load(q=-15000, element_id='all', direction='element')  # 15 kN/m downward# --- Solve ---ss.solve()# --- Post-processing ---print("n=== Support Reactions ===")for node_id, reaction in ss.get_support_force_sum().items():    print(f"  Node {node_id}: Fy = {reaction:.2f} N")# --- Plot results ---fig, axes = plt.subplots(3, 1, figsize=(12, 9))ss.show_bending_moment(figure=fig, ax=axes[0], factor=1e-6, title="Bending Moment (kN·m)")ss.show_shear_force(   figure=fig, ax=axes[1], title="Shear Force (kN)")ss.show_displacement(  figure=fig, ax=axes[2], factor=100, title="Displacement (×100)")plt.tight_layout()plt.savefig('continuous_beam.png', dpi=150)plt.show()

3.6 SciPy + NumPy (Direct FEM Implementation)

For researchers requiring complete control over the simulation pipeline, SciPy and NumPy can be combined to implement custom finite element solvers. This approach is ideal for pedagogical purposes, prototyping novel element formulations, and integrating with machine learning pipelines.

Sample Code 3.6 — Direct FEM: 2D Plane Stress with Custom Assembler

import numpy as npfrom scipy.sparse import lil_matrixfrom scipy.sparse.linalg import spsolveimport matplotlib.pyplot as pltimport matplotlib.tri as mtri# ─────────────────────────────────────────# CST ELEMENT: Constant Strain Triangle# Plane stress formulation# ─────────────────────────────────────────def get_B_matrix(coords):    """    Strain-displacement matrix B for CST element.    coords: (3,2) array of node coordinates    Returns: B (3,6), area    """    x = coords[:, 0]    y = coords[:, 1]        # Shape function derivatives    b = np.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]])    c = np.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]])        area = 0.5 * abs((x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]))        B = np.zeros((3, 6))    for i in range(3):        B[0, 2*i  ] = b[i]        B[1, 2*i+1] = c[i]        B[2, 2*i  ] = c[i]        B[2, 2*i+1] = b[i]        return B / (2 * area), areadef get_D_matrix(E, nu):    """Plane stress constitutive matrix."""    factor = E / (1 - nu**2)    return factor * np.array([        [1,  nu, 0           ],        [nu, 1,  0           ],        [0,  0,  (1-nu)/2   ]    ])def stiffness_matrix(coords, D, thickness=1.0):    """Local stiffness matrix for CST element."""    B, area = get_B_matrix(coords)    return thickness * area * (B.T @ D @ B)# ─────────────────────────────────────────# MESH GENERATION: Square plate 1m × 1m# Structured triangular mesh# ─────────────────────────────────────────def generate_mesh(Lx, Ly, nx, ny):    """Structured triangular mesh for rectangle."""    dx = Lx / nx    dy = Ly / ny        # Nodes    nodes = []    for j in range(ny + 1):        for i in range(nx + 1):            nodes.append([i * dx, j * dy])    nodes = np.array(nodes)        # Elements (2 triangles per rectangle)    elements = []    for j in range(ny):        for i in range(nx):            n0 = j * (nx + 1) + i            n1 = n0 + 1            n2 = n0 + (nx + 1)            n3 = n2 + 1            elements.append([n0, n1, n3])  # lower triangle            elements.append([n0, n3, n2])  # upper triangle        return nodes, np.array(elements)# ─────────────────────────────────────────# PROBLEM SETUP# ─────────────────────────────────────────E  = 200e9    # Panu = 0.3t  = 0.01     # m (thickness)Lx, Ly = 2.0, 1.0   # plate dimensionsnx, ny = 20, 10      # mesh divisionsnodes, elements = generate_mesh(Lx, Ly, nx, ny)n_nodes    = len(nodes)n_elements = len(elements)n_dof      = 2 * n_nodesD = get_D_matrix(E, nu)print(f"Mesh: {n_nodes} nodes, {n_elements} elements, {n_dof} DOFs")# ─────────────────────────────────────────# GLOBAL ASSEMBLY# ─────────────────────────────────────────K = lil_matrix((n_dof, n_dof))f = np.zeros(n_dof)for elem in elements:    coords = nodes[elem]  # (3,2)    Ke = stiffness_matrix(coords, D, t)        # Degree-of-freedom mapping    dofs = np.array([2*elem[0], 2*elem[0]+1,                     2*elem[1], 2*elem[1]+1,                     2*elem[2], 2*elem[2]+1])        for i, di in enumerate(dofs):        for j, dj in enumerate(dofs):            K[di, dj] += Ke[i, j]# ─────────────────────────────────────────# BOUNDARY CONDITIONS# ─────────────────────────────────────────tol = 1e-10fixed_dofs = []for i, (x, y) in enumerate(nodes):    if x < tol:  # Left edge: fix x-displacement (roller)        fixed_dofs.append(2*i)    if x < tol and abs(y - Ly/2) < tol:  # Center of left: fix y        fixed_dofs.append(2*i + 1)# Apply traction on right edge: 50 MPa tensionP = 50e6 * t  # N/m distributed as nodal forcesright_nodes = [i for i, (x, y) in enumerate(nodes) if abs(x - Lx) < tol]right_nodes.sort(key=lambda i: nodes[i, 1])for k, ni in enumerate(right_nodes):    if k == 0 or k == len(right_nodes)-1:        f[2*ni] += P * (Ly/ny) / 2   # half segment at ends    else:        f[2*ni] += P * (Ly/ny)       # full segment at interior# ─────────────────────────────────────────# APPLY BCS AND SOLVE# ─────────────────────────────────────────K_csr = K.tocsr()all_dofs  = np.arange(n_dof)free_dofs = np.setdiff1d(all_dofs, fixed_dofs)K_red = K_csr[np.ix_(free_dofs, free_dofs)]f_red = f[free_dofs]u_red = spsolve(K_red, f_red)u = np.zeros(n_dof)u[free_dofs] = u_red# ─────────────────────────────────────────# POST-PROCESSING: Stress calculation# ─────────────────────────────────────────sigma_x   = np.zeros(n_nodes)node_count = np.zeros(n_nodes)for elem in elements:    coords = nodes[elem]    B, _ = get_B_matrix(coords)        dofs = np.array([2*elem[0], 2*elem[0]+1,                     2*elem[1], 2*elem[1]+1,                     2*elem[2], 2*elem[2]+1])    u_elem = u[dofs]    eps    = B @ u_elem    sig    = D @ eps        for ni in elem:        sigma_x[ni]    += sig[0]        node_count[ni] += 1sigma_x /= node_count  # Average at nodes# ─────────────────────────────────────────# VISUALIZATION# ─────────────────────────────────────────triang = mtri.Triangulation(nodes[:, 0], nodes[:, 1], elements)fig, axes = plt.subplots(1, 2, figsize=(14, 5))# Displacement contouru_x = u[::2]ax1 = axes[0]tcf1 = ax1.tricontourf(triang, u_x * 1000, levels=20, cmap='RdBu_r')plt.colorbar(tcf1, ax=ax1, label='u_x (mm)')ax1.set_title('X-Displacement Field')ax1.set_aspect('equal')ax1.set_xlabel('x (m)'); ax1.set_ylabel('y (m)')# Stress contourax2 = axes[1]tcf2 = ax2.tricontourf(triang, sigma_x / 1e6, levels=20, cmap='jet')plt.colorbar(tcf2, ax=ax2, label='σ_x (MPa)')ax2.set_title('X-Normal Stress Field')ax2.set_aspect('equal')ax2.set_xlabel('x (m)'); ax2.set_ylabel('y (m)')plt.tight_layout()plt.savefig('plate_stress_field.png', dpi=150, bbox_inches='tight')plt.show()