Metadata-Version: 2.4
Name: beamfeapy
Version: 0.2.0
Summary: 3D FEM beam solver: Euler-Bernoulli & Timoshenko beams, tapered sections, releases, thermal/prestress loads, modal & buckling analysis, section groups, Excel I/O & Plotly plots
Author: Domenico Gaudioso
License-Expression: MIT
Project-URL: Homepage, https://github.com/DomenicoGaudioso/beamfeapy
Project-URL: Repository, https://github.com/DomenicoGaudioso/beamfeapy
Project-URL: Documentation, https://github.com/DomenicoGaudioso/beamfeapy#readme
Project-URL: Issues, https://github.com/DomenicoGaudioso/beamfeapy/issues
Keywords: FEM,beam,frame,structural,Euler-Bernoulli,Timoshenko,prestress,tapered,finite-element
Classifier: Development Status :: 4 - Beta
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Typing :: Typed
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Provides-Extra: plot
Requires-Dist: plotly>=6.1.1; extra == "plot"
Requires-Dist: kaleido>=1.0; extra == "plot"
Provides-Extra: excel
Requires-Dist: pandas>=1.5; extra == "excel"
Requires-Dist: openpyxl>=3.1; extra == "excel"
Provides-Extra: all
Requires-Dist: plotly>=6.1.1; extra == "all"
Requires-Dist: kaleido>=1.0; extra == "all"
Requires-Dist: pandas>=1.5; extra == "all"
Requires-Dist: openpyxl>=3.1; extra == "all"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: plotly>=6.1.1; extra == "dev"
Requires-Dist: kaleido>=1.0; extra == "dev"
Requires-Dist: pandas>=1.5; extra == "dev"
Requires-Dist: openpyxl>=3.1; extra == "dev"
Dynamic: license-file

# beamfeapy

<div align="center">
  <img src="img/Design%206.png" alt="beamfeapy Logo" width="216">
</div>

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

```bash
# 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

```python
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

```python
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

```python
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

```python
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

```python
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]`)

```python
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]`)

```python
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

```python
# 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

```python
# 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

```python
# 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

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

### Prestress

```python
# 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)

```python
# 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

```python
# 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)

```python
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

```python
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`:

```python
# 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:

```python
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:

```bash
cd FEM
python usage_examples/01_cantilever_nodal_load.py
```

## Testing

```bash
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`.
