beamfeapy 0.2.0

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

beamfeapy

A Python finite-element solver for the static 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 I/O, and Plotly visualization.

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)
• Plotly plots: loads (per load case), internal-force diagrams (European convention), deformed shape, reactions
• Sparse solver (COO→CSR assembly + scipy SuperLU) for large models

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]

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.