beamfeapy 0.3.0

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
• Modal & buckling analysis: natural frequencies/periods/mode shapes and critical load multipliers (λ), with optional per-case section groups
• Plotly plots: loads (per load case), internal-force diagrams (European convention), deformed shape, reactions
• Sparse solver (COO→CSR assembly + scipy SuperLU) for large models
• 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

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]

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)

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/                   # performance benchmark vs PyNite
└── 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

MIT — see LICENSE.