beamfeapy 0.4.1

3D FEM beam solver: Euler-Bernoulli & Timoshenko beams, tapered sections, releases, thermal/prestress loads, modal & buckling analysis, section groups, Excel & HDF5 I/O, Plotly plots

beamfeapy

A Python finite-element solver for the static, modal and buckling analysis of 3D frame structures composed of Euler-Bernoulli and Timoshenko beams — including tapered (non-prismatic) elements, end releases, thermal loads, prestress, nodal settlements, load cases, Excel/HDF5 I/O, export to external solvers, Plotly visualization and a Streamlit web UI.

📖 Documentation

Full documentation is published as a bilingual website (same topics in both languages):

• 🌐 Site: https://domenicogaudioso.github.io/beamfeapy/
• 🇬🇧 English: https://domenicogaudioso.github.io/beamfeapy/en.html
• 🇮🇹 Italiano: https://domenicogaudioso.github.io/beamfeapy/it.html

The source pages live in docs/ (Jekyll + just-the-docs).

Features

• 3D Euler-Bernoulli beam element (12 DOFs per element: axial + bi-axial bending + torsion)
• Timoshenko beam element (shear deformability), shear=True
• Tapered beam element (non-prismatic / variable section) — exact force-based stiffness, one element with no mesh required
• Nodal loads (forces and moments)
• Concentrated loads in span (force/moment at internal point ξ∈[0,1])
• Distributed loads: uniform, partial, trapezoidal (forces and distributed moments)
• Thermal loads: uniform, linear gradient, and generic nonlinear profiles along the section height (eigenstress, EN 1991-1-5)
• Nodal settlements (imposed displacements/rotations)
• Prestress: equivalent-load method for parabolic/straight/eccentric cables, and 3D cable geometry (anchor + deviation forces from polyline)
• End releases (hinges / internal releases via static condensation)
• Support reactions and internal forces (N, Vy, Vz, T, My, Mz)
• Load cases: assign loads to cases; solve single cases or combinations with multiplicative coefficients (solve(cases={"G": 1.35, "Q": 1.5}))
• Modal analysis: masses derived from chosen load cases (distributed + concentrated) with coefficients; natural frequencies, periods, mode shapes, mass participation (validated vs OpenSees)
• Excel I/O (Node/Material/Section/Element/Support/Load sheets, WOBridge-style)
• Results I/O: save & read results in Excel (.xlsx, static) or HDF5 (.h5, big-data format for static + modal + buckling in one file)
• Model export to external solvers: OpenSees (TCL & OpenSeesPy), SAP2000 (.s2k), MIDAS (.mct), Robot (.str), Straus7 (text) — with exact local-axis mapping
• Springs & cables: axial springs (bilateral / tension-only / compression-only), ground elastic supports, and cable elements (tension-only co-rotational bar with Ernst modulus, and exact elastic catenary) solved with a non-linear Newton scheme
• Modal & buckling analysis: natural frequencies/periods/mode shapes and critical load multipliers (λ), with optional per-case section groups — also with springs and cables present (cables linearized about their non-linear equilibrium, so the analysis always runs). Cable self-weight is included in the mass; the suspension-bridge example matches the classical beam-under-tension theory to ≈2%
• Plotly plots: loads (per load case), internal-force diagrams (European convention), deformed shape, reactions
• Sparse solver (COO→CSR assembly + scipy SuperLU, optional MKL Pardiso via pip install beamfeapy[fast]) for large models, with vectorized batched assembly, closed-form distributed loads and sparse ARPACK eigensolvers for modal & buckling — within ~2× of OpenSees on large models; see Performance & Benchmarks
• solve_many: solve many load combinations reusing a single factorization
• Web UI (Streamlit) (app.py): import a model from Excel, edit it in tables, run analyses and visualize 3D results in the browser — see the Web UI guide
• Command-line interface (beamfeapy / bfp): run solve, modal, buckling, export and report directly from an Excel model, with a colored ASCII logo — the python snake as a simply-supported beam on two tree trunks, body shaded like its bending-moment diagram (beamfeapy logo) — see Command-Line Interface

Installation

# From source (development), with all extras (plot + Excel):
pip install -e ".[all]"

# Base only (numpy, scipy):
pip install -e .

Optional extras: plot (Plotly + kaleido), excel (pandas + openpyxl), all, dev.

Once published on PyPI: pip install beamfeapy[all]

Requirements: Python ≥ 3.9, numpy ≥ 1.24, scipy ≥ 1.10

Quick Start

from beamfeapy import Model, Material, Section

m = Model()
m.add_node(1, 0, 0, 0)
m.add_node(2, 4, 0, 0)

mat = Material(E=210e9, nu=0.3, alpha=1.2e-5)
sec = Section(A=1e-2, Iy=2e-5, Iz=3e-5, J=1e-5)
m.add_beam(1, 1, 2, mat, sec)

m.fix(1)                        # fixed support at node 1
m.add_nodal_load(2, Fy=-10000)  # vertical force at tip

res = m.solve()
print(res.displacements(2))     # [ux, uy, uz, rx, ry, rz]
print(res.reactions(1))         # [Fx, Fy, Fz, Mx, My, Mz]

Command-Line Interface (CLI)

Installing the package adds a stylized beamfeapy command (alias bfp) so you can run analyses straight from an Excel model — no Python code required:

beamfeapy logo                                 # show the colored app logo
beamfeapy new model.xlsx                        # build a model with a guided wizard
beamfeapy template model.xlsx                 # generate an input template
beamfeapy info model.xlsx                      # model summary (nodes, loads, cases)
beamfeapy solve model.xlsx -o out.xlsx --cases G=1.35 Q=1.5
beamfeapy modal model.xlsx -n 12 --mass-cases G=1.0 Q=0.3 -o modal.h5
beamfeapy buckling model.xlsx --cases P
beamfeapy plot model.xlsx --what deformed --cases G=1.35 Q=1.5 --open   # interactive HTML, opens in browser
beamfeapy plot model.xlsx --what forces --component Mz --png            # static PNG instead
beamfeapy export model.xlsx out.tcl            # → OpenSees / SAP2000 / MIDAS / Robot / Straus7
beamfeapy report model.xlsx -o report.docx     # Word calc report

Build a model from scratch without leaving the terminal with beamfeapy new model.xlsx: a guided wizard asks for nodes, material, section, elements, supports and nodal loads, then writes the Excel model — ready to solve or plot.

Visualize the model and results with beamfeapy plot. --what selects model (geometry + supports), loads, deformed (deflected shape, auto-scaled), forces (internal-force diagrams N/V/T/M, pick one with --component Mz), reactions, or all (writes every figure to a folder). By default it writes an interactive 3D Plotly HTML and prints a clickable file:// link to open in the browser (add --open to launch it automatically); pass --png for a static Matplotlib image instead:

Run beamfeapy logo to play the animated logo: a little rainbow python that slithers across the terminal from right to left, with the beamfeapy wordmark shimmering below; it animates for a few seconds in an interactive terminal — use --loop to keep it going, --still for the static logo. The static logo is the python snake drawn as a simply-supported beam resting on two cut-trunk supports. Its body carries the FEM stress-contour gradient — blue at the supports, red at midspan — so it doubles as the bending-moment diagram of a simply-supported beam, while the beamfeapy wordmark uses the same blue→red gradient. A compact "mini snake" header precedes the operational commands. beamfeapy logo --image shows instead the photo-realistic logo rendered in truecolor half-blocks (from img/2.png, embedded compressed via scripts/_gen_cli_logo.py). beamfeapy <command> --help shows per-command options.

Load combinations are passed as CASE=COEF pairs (G=1.35 Q=1.5 → dict with multiplicative coefficients) or plain names (G Q). Output is colorized automatically (24-bit ANSI, with Virtual-Terminal support on Windows); use --color never or set NO_COLOR to disable (falls back to a stylized box banner), and --no-banner to hide the header. The CLI needs no image library at runtime.

Web UI (Streamlit)

Prefer a graphical interface? A Streamlit app (app.py) lets you import a model from Excel, edit it in editable tables, run analyses (static / modal / buckling) and visualize 3D results in the browser:

pip install beamfeapy streamlit plotly openpyxl
streamlit run app.py

See the Web UI guide (italiano) for the full walkthrough with screenshots.

API Reference

Model Construction

Method	Description
Model()	Create an empty model
m.add_node(id, x, y, z)	Add a node with ID and coordinates
m.add_beam(id, ni, nj, mat, sec, ...)	Add a 3D beam element
m.add_tapered_beam(id, ni, nj, mat, ...)	Add a tapered (variable-section) beam element
m.add_section(id, A=..., Iy=..., ...)	Register a section by ID for reuse

Materials & Sections

mat = Material(E=210e9, nu=0.3, alpha=1.2e-5)   # steel, SI
mat = Material(E=30e9, nu=0.2, alpha=1.0e-5)     # concrete

sec = Section(A=1e-2, Iy=2e-5, Iz=3e-5, J=1e-5)                # basic (EB)
sec = Section(A=0.18, Iy=5.4e-3, Iz=1.35e-3, J=2e-3,
              Asy=5/6*0.18, Asz=5/6*0.18)                       # Timoshenko
sec = Section(A=..., Iy=..., Iz=..., J=..., h_y=0.6, h_z=0.3)  # with heights for thermal

Section parameters:

Parameter	Description
A	Cross-sectional area
Iy	Moment of inertia about local y-axis (bending in x-z plane)
Iz	Moment of inertia about local z-axis (bending in x-y plane)
J	Torsional constant
Asy, Asz	Effective shear areas (Timoshenko only)
h_y, h_z	Section heights along local y/z (for thermal gradients)

Supports

m.fix(1)                                    # fixed (all 6 DOFs restrained)
m.pin(2)                                    # pin (3 translations restrained)
m.support(3, ux=True, uy=True, rz=True)    # custom: restrain specific DOFs

Loads

Method	Description
m.add_nodal_load(node, Fx=..., case=...)	Force/moment at a node (global axes)
m.add_distributed_load(elem, comp, qi, [qj], [a], [b], [frame], [case])	Distributed load (uniform/partial/trapezoidal)
m.add_concentrated_load(elem, xi, Fx=..., frame=..., case=...)	Concentrated load at internal point ξ∈[0,1]
m.add_thermal_load(elem, dT_axial=..., dT_grad_y=..., dT_grad_z=..., case=...)	Thermal load (uniform + linear gradient)
m.add_thermal_profile(elem, profile, axis=..., width=..., case=...)	Nonlinear thermal profile along section height
m.add_prestress(elem, P, e_i=..., e_j=..., plane=..., sag=..., case=...)	Prestress via equivalent loads
m.add_cable_prestress(P, points, [elements], case=...)	Prestress from 3D cable geometry
m.add_settlement(node, dof, value)	Imposed displacement at a node

Solution & Results

res = m.solve()                    # dense solver (default)
res = m.solve(sparse=True)         # sparse solver (large models)
res = m.solve(cases=["G", "Q"])    # load-case combination
res = m.solve(cases="G")           # single load case

res.displacements(node)            # array [ux, uy, uz, rx, ry, rz]
res.displacement(node, "uy")       # single DOF value
res.reactions(node)                 # array [Fx, Fy, Fz, Mx, My, Mz]
res.element_forces[elem_id]        # end forces in local coords (12,)

Post-processing

from beamfeapy import postprocess

di = postprocess.internal_forces(res, elem_id, n=101)
# Returns dict with keys: x, N, Vy, Vz, T, My, Mz

dd = postprocess.element_displacements(res, elem_id, n=51)
# Returns dict with keys: x, u_local (n×6 array)

pts = postprocess.deformed_shape_global(res, elem_id, n=51, scale=100)
# Returns n×3 array of global coordinates (deformed, scaled)

Plotting (requires pip install beamfeapy[plot])

from beamfeapy.plotting import (
    plot_model, plot_loads, plot_diagram,
    plot_deformed, plot_reactions, plot_internal_forces,
)

plot_loads(m, case="G").show()        # structure + loads for load case "G"
plot_diagram(res, "Mz").show()        # Mz diagram (European convention)
plot_deformed(res, scale=200).show()   # deformed shape
plot_reactions(res).show()             # support reactions
plot_internal_forces(res, elem_id=1)   # all 6 diagrams for one element

• plot_diagram(result, component) — component is one of: N, Vy, Vz, T, My, Mz
• European convention: negative moment drawn at the extrados (top)

Excel I/O (requires pip install beamfeapy[excel])

from beamfeapy import Model, read_excel
from beamfeapy.io_excel import write_template

write_template("input.xlsx")             # generate a fillable template
m = read_excel("input.xlsx")             # = Model.from_excel("input.xlsx")
res = m.solve()
res.to_excel("results.xlsx", n_diagram=21)  # displacements, reactions, forces, diagrams

Detailed Feature Guide

Distributed Loads

# Uniform load over the entire element (local y, -3 kN/m)
m.add_distributed_load(1, "fy", -3000)

# Partial load from x=1m to x=3m (local y, -5 kN/m)
m.add_distributed_load(1, "fy", -5000, a=1.0, b=3.0)

# Trapezoidal load: from 0 to -8 kN/m over the full element
m.add_distributed_load(1, "fy", 0, -8000)

# Trapezoidal partial load from x=2m to x=5m
m.add_distributed_load(1, "fy", -2000, -8000, a=2.0, b=5.0)

# Distributed moments (torsion, bending)
m.add_distributed_load(1, "mx", 500)                    # torsional moment
m.add_distributed_load(1, "mz", 0, -1200, a=2, b=5)     # partial trapezoidal bending moment

# Global-frame distributed load
m.add_distributed_load(1, "fy", -3000, frame="global")

Component ∈ {fx, fy, fz, mx, my, mz}.
Frame: "local" (default) or "global".

Concentrated Loads in Span

# Force at 35% of element length
m.add_concentrated_load(1, 0.35, Fy=-50000)

# Moment at 70% of element length
m.add_concentrated_load(1, 0.70, Mz=80000)

# Force in global coordinates
m.add_concentrated_load(1, 0.5, Fz=-20000, frame="global")

Thermal Loads

# Uniform temperature increase (+20°C)
m.add_thermal_load(1, dT_axial=20.0)

# Temperature gradient along z (requires section.h_z)
m.add_thermal_load(1, dT_axial=20.0, dT_grad_z=15.0)

# Nonlinear thermal profile (EN 1991-1-5 style)
m.add_thermal_profile(1, lambda s: 15 * (0.5 + s / 0.3), axis="z")
# Or from discrete points:
m.add_thermal_profile(1, [(-0.15, 0), (0.05, 2.5), (0.15, 15)], axis="z", width=0.30)

Settlements

# Vertical settlement of 5 mm at node 2
m.add_settlement(2, "uz", -0.005)

Prestress

# Parabolic cable: P=2 MN, sag=0.35 m, eccentricity zero at ends
m.add_prestress(1, P=2.0e6, sag=0.35)

# Straight eccentric cable: e=0.20 m (constant)
m.add_prestress(2, P=1.5e6, e_i=0.20, e_j=0.20, plane="z")

# Generic eccentricity profile e(ξ):
m.add_prestress(3, P=1.0e6, profile=lambda xi: 0.3 * (1 - (2*xi - 1)**2))

# 3D cable geometry
pts = [(0, 0.3, 0), (5, -0.05, 0), (10, 0.3, 0)]
m.add_cable_prestress(P=3.0e6, points=pts)

End Releases (Hinges)

# Moment-free (hinge) at end j (Mz = 0)
m.add_beam(1, 1, 2, mat, sec, releases_j=["rz"])

# Releases at both ends
m.add_beam(2, 2, 3, mat, sec, releases_i=["rz"], releases_j=["ry", "rz"])

Allowed release names: ux, uy, uz, rx, ry, rz (local DOFs).

Timoshenko Beam

# Shear-deformable beam: provide effective shear areas
sec = Section(A=0.18, Iy=5e-3, Iz=2e-3, J=3e-3,
              Asy=5/6*0.18, Asz=5/6*0.18)
m.add_beam(1, 1, 2, mat, sec, shear=True)

Tapered Beam (Variable Section)

from beamfeapy import VariableSection

# Method 1: continuous function via VariableSection
vs = VariableSection.rectangular(b=0.30, h=lambda xi: 0.70 * (1 - 0.6 * xi))
m.add_tapered_beam(1, 1, 2, mat, vs)

# Method 2: sections at i and j ends (linear interpolation)
m.add_section(1, A=1.2e-2, Iy=4e-5, Iz=6e-5, J=3e-5)
m.add_section(2, A=0.6e-2, Iy=1e-5, Iz=1.5e-5, J=0.8e-5)
m.add_tapered_beam(1, 1, 2, mat, section_i=1, section_j=2)

# Method 3: sections at multiple stations
m.add_section(3, A=1.0e-2, Iy=3e-5, Iz=4e-5, J=2e-5)
m.add_tapered_beam(2, 2, 3, mat, stations={0.0: 1, 0.4: 3, 1.0: 2})

Load Cases & Combinations

m.add_distributed_load(2, "fy", -20e3, case="G")   # permanent loads
m.add_nodal_load(2, Fx=30e3, case="Q")              # variable loads
m.add_thermal_load(1, dT_axial=15, case="T")         # thermal

m.load_cases()                    # → ['G', 'Q', 'T']
res = m.solve(cases=["G", "Q"])   # combine G + Q
res = m.solve(cases="G")          # single case
res = m.solve()                    # all cases combined

Section Orientation (ref_vector and roll)

By default, the local y-axis is determined automatically. For 3D frames (e.g. portal in the x-y plane), use ref_vector or roll:

# Portal frame in x-y plane: local z = global Z
m.add_beam(2, 2, 3, mat, sec, ref_vector=(0, 0, 1))

# Roll angle (radians) rotates the section about x-local
m.add_beam(1, 1, 2, mat, sec, roll=0.15)

Convention: e_y = ref_vector × e_x (right-hand rule). The roll angle rotates the section about the local x-axis after the default orientation.

Sparse Solver

For models with many DOFs, use the sparse solver for significant speed and memory savings:

res = m.solve(sparse=True)

# many combinations, a single factorization
results = m.solve_many({"ULS": {"G": 1.35, "Q": 1.5}, "SLS": {"G": 1, "Q": 1}})

# sparse eigensolvers (default for large models)
modal = m.modal(n_modes=6, mass_source={"G": 1.0}, sparse=True)
buck  = m.buckling(n_modes=4, cases="G", sparse=True)

Install the optional MKL Pardiso backend for an extra speed-up on large systems:

pip install beamfeapy[fast]   # uses pypardiso automatically when available

See Performance & Benchmarks for the full comparison against OpenSees, PyNite and anastruct.

Conventions

• Nodal DOFs: [ux, uy, uz, rx, ry, rz] (global system)
• Local axes: x from node i to node j; y, z are principal axes of the section.
  Iz → bending in x-y plane; Iy → bending in x-z plane.
• Normal force: positive in tension.
• Units: user's choice (consistent, e.g. SI: N, m, Pa).

Project Structure

FEM/
├── beamfeapy/                  # core library
│   ├── material.py              # Material (E, nu, alpha, G)
│   ├── section.py               # Section (A, Iy, Iz, J, Asy, Asz, h_y, h_z)
│   ├── node.py                  # Node (id, x, y, z)
│   ├── element.py               # BeamElement3D (stiffness, transformation, releases)
│   ├── tapered.py               # VariableSection, TaperedBeamElement3D
│   ├── loads.py                 # NodalLoad, DistributedLoad, ConcentratedLoad, Thermal, Prestress, Settlement
│   ├── integration.py           # Gauss-Legendre quadrature
│   ├── model.py                 # Model: assembly, constraints, load cases, solution, Result
│   ├── postprocess.py           # internal forces, deformed shape
│   ├── io_excel.py              # Excel I/O (read_model, write_results, write_template)
│   └── plotting/                # Plotly visualizations
├── examples/                   # basic examples (save HTML to output/)
├── usage_examples/              # comprehensive usage examples (see below)
├── tests/                      # pytest tests vs analytical solutions
├── validation/                  # validation scripts vs OpenSees / analytical
├── benchmark/                   # benchmarks vs OpenSees / PyNite / CALFEM / pystran / anastruct (static, modal, buckling)
└── pyproject.toml               # packaging (pip install -e ".[all]")

Usage Examples

The usage_examples/ directory contains self-contained scripts covering every feature:

#	File	Description
01	01_cantilever_nodal_load.py	Cantilever beam with tip force — the "hello world"
02	02_simply_supported_beam.py	Simply supported beam with uniform distributed load
03	03_distributed_loads.py	Partial and trapezoidal distributed loads
04	04_concentrated_loads_in_span.py	Concentrated forces and moments at internal points
05	05_thermal_loads.py	Uniform and gradient thermal loads
06	06_thermal_profile.py	Nonlinear thermal profile (EN 1991-1-5, eigenstress)
07	07_settlements.py	Imposed displacements (settlements)
08	08_prestress_parabolic.py	Prestress with parabolic cable (equivalent loads)
09	09_prestress_cable_geometry.py	Prestress from 3D cable polyline geometry
10	10_timoshenko.py	Timoshenko vs Euler-Bernoulli comparison
11	11_end_releases.py	End releases (hinges, pins, partial restraints)
12	12_tapered_beam.py	Tapered (variable-section) beam, single element
13	13_tapered_beam_stations.py	Tapered beam with section IDs and multiple stations
14	14_3d_portal_frame.py	3D portal frame: columns + beam, multiple load types
15	15_load_cases.py	Load cases and combinations
16	16_ref_vector_and_roll.py	Section orientation with ref_vector and roll angle
17	17_support_types.py	Fix, pin, roller, and custom support conditions
18	18_internal_forces.py	Post-processing internal forces (N, V, M, T diagrams)
19	19_plotting.py	All Plotly visualization functions
20	20_excel_io.py	Excel input/output workflow
21	21_sparse_solver.py	Sparse solver for large models
22	22_continuous_beam.py	Multi-span continuous beam with various load patterns
23	23_frame_with_hinges.py	Frame with internal hinges and distributed loads
24	24_prestress_secondary_moments.py	Prestress in a hyperstatic structure (secondary moments)

Run any example:

cd FEM
python usage_examples/01_cantilever_nodal_load.py

Testing

pip install -e ".[dev]"
python -m pytest tests -q

Tests verify results against known analytical solutions: cantilever deflection (PL³/3EI), simply supported beam with uniform load (5qL⁴/384EI), thermal expansion, Timoshenko shear deformability, prestress, and mode shapes (validated vs OpenSees).

Validation

See validation/ for cross-validation against OpenSees and analytical benchmarks (errors at machine precision, ≈1e-10%).

License

PolyForm Noncommercial License 1.0.0 — see LICENSE.

beamfeapy is source-available: you may use, modify and share it freely for any noncommercial purpose (personal projects, research, teaching, public institutions, etc.). Commercial use requires a separate license — contact the author. See https://polyformproject.org/licenses/noncommercial/1.0.0 for the full terms.