Assemble and verify a beam finite-element model
Assemble Hermite beam elements, enforce the clamp and compare nodal deflection and reactions with a closed reference.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from IPython.display import display as _display
def display(value):
if isinstance(value,pd.DataFrame):
table=value.to_html(max_rows=12,escape=True)
style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
return _display(value)
RECIPE_ID='en04-beam-fem'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['models']
Read the method functions
def cantilever_fem(elements,length_m,young_Pa,second_moment_m4,tip_force_N):
"""Euler–Bernoulli Hermite beam elements: transverse displacement and rotation DOFs."""
if not isinstance(elements,(int,np.integer)) or elements<1 or min(length_m,young_Pa,second_moment_m4)<=0:
raise ValueError('Need positive beam geometry/material and an integer element count')
nodes=elements+1;ell=length_m/elements;stiffness=np.zeros((2*nodes,2*nodes));load=np.zeros(2*nodes)
local=young_Pa*second_moment_m4/ell**3*np.array([[12,6*ell,-12,6*ell],[6*ell,4*ell**2,-6*ell,2*ell**2],[-12,-6*ell,12,-6*ell],[6*ell,2*ell**2,-6*ell,4*ell**2]])
for element in range(elements):
indices=np.arange(2*element,2*element+4);stiffness[np.ix_(indices,indices)]+=local
load[-2]=tip_force_N;displacement=np.zeros(2*nodes)
displacement[2:]=np.linalg.solve(stiffness[2:,2:],load[2:])
reaction=stiffness@displacement-load
return dict(x=np.linspace(0,length_m,nodes),displacement_m=displacement[::2],rotation_rad=displacement[1::2],stiffness=stiffness,load=load,all_DOFs=displacement,reaction=reaction)
print("Method ready: cantilever_fem")
Method ready: cantilever_fem
Data and model boundary
The building record contains real ten-minute energy and environmental observations from one source building. Airport weather fields include upstream interpolation and are not colocated calibration references. Beam, circuit and PID examples are declared numerical design models; computed responses do not establish hardware performance or commissioning acceptance.
Steps
1. Assemble material, geometry and boundary conditions
Each Hermite beam node has transverse displacement (m) and rotation (rad). The clamp sets both root DOFs to zero. Nominal E and I are declared teaching parameters; the model assumes a slender, linear elastic beam and small displacement.
p=snapshot("models")["beam"]
beam=cantilever_fem(10,p["length_m"],p["young_modulus_Pa"],p["second_moment_m4"],p["tip_force_N"])
dense=np.linspace(0,p["length_m"],501)
analytic=p["tip_force_N"]*dense**2*(3*p["length_m"]-dense)/(6*p["young_modulus_Pa"]*p["second_moment_m4"])
fig,ax=plt.subplots(figsize=(8,4));ax.plot(dense,1000*analytic,label="Closed tip-load solution")
ax.plot(beam["x"],1000*beam["displacement_m"],"o",label="Finite-element nodes")
ax.set(xlabel="Beam coordinate (m)",ylabel="Transverse deflection (mm)",title="Declared cantilever tip-load model");ax.legend();plt.tight_layout();plt.show()
2. Inspect the assembled connectivity
The stiffness sparsity pattern reflects adjacent-element coupling before the clamp is removed. Because displacement and rotation DOFs have different units, a raw condition number depends on scaling and is not a material property.
fig,ax=plt.subplots(figsize=(6,5));ax.spy(beam["stiffness"],markersize=3)
ax.set(xlabel="DOF column",ylabel="DOF row",title="Assembled beam stiffness connectivity");plt.tight_layout();plt.show()
3. Check meshes for this exactly represented load case
The closed deflection is cubic, matching the Hermite basis, so nodal tip-load results are exact up to rounding even on a coarse mesh. This special case is a correctness benchmark, not proof that arbitrary loads, geometries or dynamic analyses need only one element.
records=[]
reference_tip=p["tip_force_N"]*p["length_m"]**3/(3*p["young_modulus_Pa"]*p["second_moment_m4"])
for count in [1,2,4,8,16]:
trial=cantilever_fem(count,p["length_m"],p["young_modulus_Pa"],p["second_moment_m4"],p["tip_force_N"])
records.append({"elements":count,"tip_deflection_m":trial["displacement_m"][-1],"reference_m":reference_tip,"relative_error":abs(trial["displacement_m"][-1]/reference_tip-1)})
mesh=pd.DataFrame(records);display(mesh)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(mesh.elements,1000*mesh.tip_deflection_m,"o-")
ax.axhline(1000*reference_tip,color="gray",linestyle="--")
ax.set(xlabel="Elements",ylabel="Tip deflection (mm)",title="Special cubic benchmark across meshes");plt.tight_layout();plt.show()
| elements | tip_deflection_m | reference_m | relative_error | |
|---|---|---|---|---|
| 0 | 1 | 0.00162 | 0.00162 | 2.220446e-16 |
| 1 | 2 | 0.00162 | 0.00162 | 2.220446e-15 |
| 2 | 4 | 0.00162 | 0.00162 | 2.375877e-14 |
| 3 | 8 | 0.00162 | 0.00162 | 1.620926e-14 |
| 4 | 16 | 0.00162 | 0.00162 | 2.161826e-12 |
4. Reconcile equilibrium reactions and export
The root reaction force balances the tip force, and the root reaction moment balances F·L under the declared sign convention. Interior residuals should vanish. These are model-equilibrium checks, not load-test acceptance.
assert np.allclose(beam["stiffness"],beam["stiffness"].T)
assert np.max(abs(beam["reaction"][2:]))<1e-6
assert np.isclose(beam["reaction"][0],-p["tip_force_N"],rtol=1e-8)
assert np.isclose(beam["reaction"][1],-p["tip_force_N"]*p["length_m"],rtol=1e-8)
assert mesh.relative_error.max()<1e-8
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"x_m":beam["x"],"deflection_m":beam["displacement_m"],"rotation_rad":beam["rotation_rad"]}).to_csv(OUTPUT_DIR/"beam-nodes.csv",index=False)
mesh.to_csv(OUTPUT_DIR/"cubic-benchmark-meshes.csv",index=False)
print("Assembly, reaction and closed-solution checks passed; exports:",OUTPUT_DIR)
Assembly, reaction and closed-solution checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/en04-beam-fem
Exercises and reference explanation
- Why does the cubic tip-load benchmark look mesh-independent?
- Why are displacement and rotation DOFs not interchangeable?
Reference explanation. Hermite elements represent the cubic closed displacement for this special load case. Displacement is in m and rotation in rad; the stiffness couples different physical quantities.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Using NumPy/Matplotlib, add the analytical rotation curve and compare with all nodal rotations.
Prompt 2
Using the same packages, vary E and I separately and check inverse-stiffness scaling while retaining force/moment equilibrium.
Human acceptance. Clamped DOFs remain zero, reactions balance F and F·L and nodal references match within numerical tolerance.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/linalg.html',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256
a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256
4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.