Add 1D Euler-Bernoulli beam solver (cubic Hermite elements)

examples/Beam1DFEM/Beam1DEuler_Bernoulli.js

@@ -0,0 +1,95 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+/**
+ * Clamped and spring-supported Euler-Bernoulli beam
+ *
+ * Reproduces the "Bending of a Beam" example from J.N. Reddy, "An Introduction to
+ * the Finite Element Method", 3rd ed., McGraw-Hill, 2006 (FEM1D example problems,
+ * Chapter 7), solved there with the reference FEM1D Fortran program (MODEL=3,
+ * NTYPE=0, IELEM=0 i.e. Hermite cubic elements).
+ *
+ * Geometry: a 10 m beam, clamped at x=0, resting on a roller at midspan (x=5 m),
+ * and connected to a linear transverse spring at the free end (x=10 m).
+ * Mesh: 2 Hermite cubic beam elements of 5 m each -> 3 nodes: [1, 2, 3]
+ *
+ *        1,000 N/m (on 0 <= x <= 5)      2,500 N (down, at node 3)
+ *        v v v v v v v v v v            |
+ *   /////|--------------------|-------------------|       ~~~~ spring, k = 1e-4*EI
+ *   ///// 1  (clamped)         2 (roller)          3 (free end, spring)
+ *         |<-------- 5 m ----->|<------- 5 m ------>|
+ *                    ^ moment 1,250 N-m applied at node 2
+ *
+ * EI = 2e6 N-m^2 (constant), k_spring = 1e-4 * EI = 200 N/m
+ *
+ * Boundary conditions (see beamBoundaryConditions.js for the condition syntax):
+ *  - Node 1 (x=0):  fixed            -> w=0, theta=0     (clamped support)
+ *  - Node 2 (x=5):  pinned + moment  -> w=0, applied moment M=1250 N-m
+ *  - Node 3 (x=10): spring + force   -> k=200 N/m, applied point load P=-2500 N
+ */
+
+// Import Math.js
+import * as math from "mathjs";
+global.math = math;
+
+// Import FEAScript library
+import { FEAScriptModel, printVersion } from "feascript";
+
+console.log("FEAScript Version:", printVersion);
+
+// Create a new FEAScript model
+const model = new FEAScriptModel();
+
+// Select physics/PDE
+model.setModelConfig("eulerBernoulliBeamScript", {
+  coefficientFunctions: {
+    EI: (x) => 2.0e6, // Bending stiffness E*I (N-m^2), constant along the beam
+    // Distributed transverse load: -1,000 N/m over the first (clamped) span only
+    q: (x) => (x <= 5 ? -1000 : 0),
+    // c0 defaults to 0 (no elastic foundation) when omitted
+  },
+});
+
+// Define mesh configuration
+// elementOrder is 'linear' because that only describes the 2-node beam geometry;
+// the field itself is always interpolated with cubic Hermite shape functions internally
+model.setMeshConfig({
+  meshDimension: "1D",
+  elementOrder: "linear",
+  numElementsX: 2,
+  maxX: 10,
+});
+
+// Define boundary conditions (keyed by 1-based global node number)
+model.addBoundaryCondition("1", [["fixed"]]); // Clamped support
+model.addBoundaryCondition("2", [["pinned"], ["moment", 1250]]); // Roller + applied moment
+model.addBoundaryCondition("3", [["spring", 200], ["force", -2500]]); // Spring support + point load
+
+// Set solver method
+model.setSolverMethod("lusolve");
+
+// Solve the problem
+const { solutionVector } = model.solve();
+
+// The solution vector is ordered [w_0, theta_0, w_1, theta_1, w_2, theta_2, ...]
+// (mathjs' lusolve returns a nested array, so flatten defensively before reading it)
+const flatSolution = solutionVector.map((entry) => (Array.isArray(entry) ? entry[0] : entry));
+
+const nodeXCoordinates = [0, 5, 10];
+console.log("\nNode |    x (m) | Deflection w (m) | Rotation theta (rad)");
+console.log("-----|----------|-------------------|----------------------");
+for (let nodeIndex = 0; nodeIndex < nodeXCoordinates.length; nodeIndex++) {
+  const w = flatSolution[2 * nodeIndex];
+  const theta = flatSolution[2 * nodeIndex + 1];
+  console.log(
+    `  ${nodeIndex + 1}  | ${nodeXCoordinates[nodeIndex].toFixed(2).padStart(8)} | ${w
+      .toExponential(4)
+      .padStart(17)} | ${theta.toExponential(4).padStart(20)}`,
+  );
+}

examples/Beam1DFEM/README.md

@@ -0,0 +1,102 @@
+<img src="https://feascript.github.io/FEAScript-website/assets/feascript-structural-mechanics.png" width="80" alt="FEAScript Beam1DFEM Logo">
+
+# 1D Euler-Bernoulli Beam Examples
+
+This directory contains Node.js examples demonstrating how to use the FEAScript library to solve
+1D beam bending problems with the Euler-Bernoulli beam theory, using cubic Hermite finite elements.
+
+## Examples
+
+#### 1. Clamped and Spring-Supported Beam (`Beam1DEuler_Bernoulli.js`)
+
+Reproduces the "Bending of a Beam" example from J.N. Reddy, _An Introduction to the Finite Element
+Method_, 3rd ed., McGraw-Hill, 2006 (FEM1D example problems, Chapter 7). A 10 m beam is clamped at
+one end, supported by a roller at midspan, and connected to a linear elastic spring at the free
+end. It carries a uniformly distributed load over the clamped span, a concentrated moment at the
+roller, and a point load at the free end.
+
+## Model overview
+
+The 1D Euler-Bernoulli beam solver (`eulerBernoulliBeamScript`) assembles the fourth-order beam
+bending equation
+
+```
+d²/dx²( EI(x) d²w/dx² ) + c0(x) w = q(x)
+```
+
+using cubic Hermite shape functions, so that both the deflection `w` and the rotation
+`theta = dw/dx` are continuous across element boundaries (C¹ continuity), as required for a
+fourth-order (curvature-based) formulation. Each node therefore carries **2 degrees of freedom**,
+ordered as `[w_0, theta_0, w_1, theta_1, ...]` in the solution vector.
+
+### Mesh configuration
+
+The beam geometry only needs each element's 2 end nodes, so `elementOrder` should be set to
+`"linear"` when configuring the mesh — this describes the geometry only. The field itself is
+always interpolated with cubic Hermite shape functions internally, independent of this setting
+(a "subparametric" element formulation, standard for beam elements).
+
+```javascript
+model.setMeshConfig({
+  meshDimension: "1D",
+  elementOrder: "linear",
+  numElementsX: 2,
+  maxX: 10,
+});
+```
+
+### Coefficient functions
+
+```javascript
+model.setModelConfig("eulerBernoulliBeamScript", {
+  coefficientFunctions: {
+    EI: (x) => 2.0e6, // Bending stiffness E*I(x) [required]
+    c0: (x) => 0, // Elastic foundation modulus [optional, defaults to 0]
+    q: (x) => -1000, // Distributed transverse load [optional, defaults to 0]
+  },
+});
+```
+
+### Boundary conditions
+
+Unlike the scalar 1D models (which tag only the two ends of the domain as boundary `"0"`/`"1"`),
+beam problems commonly need essential and natural conditions at interior nodes as well (e.g. a
+midspan support, or a point load applied partway along the beam). Boundary conditions for the beam
+solver are therefore keyed by the **1-based global node number**, and each key maps to an **array**
+of condition tuples, since a node can carry more than one condition at once (e.g. a spring support
+plus a point load):
+
+```javascript
+model.addBoundaryCondition("1", [["fixed"]]); // w=0, theta=0 (clamped)
+model.addBoundaryCondition("2", [["pinned"], ["moment", 1250]]); // w=0, plus an applied moment
+model.addBoundaryCondition("3", [["spring", 200], ["force", -2500]]); // elastic support, plus a point load
+```
+
+| Condition type                      | Kind             | Effect                                                    |
+| ------------------------------------ | ---------------- | ----------------------------------------------------------- |
+| `["fixed"]`                          | Essential         | `w = 0` and `theta = 0` (clamped support)                   |
+| `["pinned"]` / `["deflection", v]`   | Essential         | `w = v` (default `v = 0`; roller/pin support)                |
+| `["rotationFixed"]` / `["rotation", v]` | Essential      | `theta = v` (default `v = 0`)                                |
+| `["force", v]`                       | Natural           | Applies a concentrated transverse force `v` at the node       |
+| `["moment", v]`                      | Natural           | Applies a concentrated moment `v` at the node                 |
+| `["spring", k, uRef]`                | Mixed (Robin)     | Transverse elastic support of stiffness `k` about `uRef` (default `uRef = 0`) |
+
+## Running the Node.js Examples
+
+#### 1. Create package.json with ES module support:
+
+```bash
+echo '{"type":"module"}' > package.json
+```
+
+#### 2. Install dependencies:
+
+```bash
+npm install feascript
+```
+
+#### 3. Run the example:
+
+```bash
+node Beam1DEuler_Bernoulli.js
+```

src/FEAScript.js

@@ -16,6 +16,7 @@ import { assembleFrontPropagationMat } from "./models/frontPropagation.js";
 import { assembleGeneralFormPDEMat, assembleGeneralFormPDEFront } from "./models/generalFormPDE.js";
 import { assembleHeatConductionMat, assembleHeatConductionFront } from "./models/heatConduction.js";
 import { assembleCreepingFlowMatrix } from "./models/creepingFlow.js";
+import { assembleEulerBernoulliBeamMat } from "./models/eulerBernoulliBeam.js";
 import { runFrontalSolver } from "./methods/frontalSolver.js";
 import { basicLog, debugLog, warnLog, errorLog } from "./utilities/logging.js";
 
@@ -212,7 +213,29 @@ export class FEAScriptModel {
         totalNodesPressure: creepingFlowResult.totalNodesPressure,
         pressureNodeIndices: creepingFlowResult.pressureNodeIndices,
       };
+    } else if (this.solverConfig === "eulerBernoulliBeamScript") {
+      // Use regular linear solver methods for the 1D Euler-Bernoulli beam model
+      const beamResult = assembleEulerBernoulliBeamMat(
+        meshData,
+        this.boundaryConditions,
+        this.coefficientFunctions,
+      );
+      jacobianMatrix = beamResult.jacobianMatrix;
+      residualVector = beamResult.residualVector;
+
+      const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
+        maxIterations: options.maxIterations ?? this.maxIterations,
+        tolerance: options.tolerance ?? this.tolerance,
+      });
+      solutionVector = linearSystemResult.solutionVector;
+
+      // Store beam-specific metadata for solution extraction (2 DOFs per node: deflection, rotation)
+      this._eulerBernoulliBeamMetadata = {
+        dofsPerNode: beamResult.dofsPerNode,
+        totalNodesX: meshData.totalNodesX,
+      };
     }
+
     console.timeEnd("totalSolvingTime");
     basicLog("Solving process completed");
 

src/mesh/basisFunctions.js

@@ -26,17 +26,26 @@ export class BasisFunctions {
 
   /**
    * Function to calculate basis functions and their derivatives based on the dimension and order
-   * @param {number} ksi - Natural coordinate (for both 1D and 2D)
+   * @param {number} ksi - Natural coordinate (for both 1D and 2D), in [0, 1]
    * @param {number} [eta] - Second natural coordinate (only for 2D elements)
+   * @param {number} [elementLength] - Physical element length, only required for 1D 'hermiteCubic' elements
    * @returns {object} An object containing:
    *  - basisFunction: Array of evaluated basis functions
    *  - basisFunctionDerivKsi: Array of derivatives of basis functions with respect to ksi
    *  - basisFunctionDerivEta: Array of derivatives of basis functions with respect to eta (only for 2D elements)
+   *  - basisFunctionDerivKsi2: Array of second derivatives of basis functions with respect to ksi
+   *    (only for 1D 'hermiteCubic' elements)
+   *
+   * 'hermiteCubic' is a general-purpose interpolation type, not specific to beams: it makes both
+   * the value and the slope continuous across elements (unlike the C0 Lagrange types above), which
+   * any fourth-order problem needs (e.g. beam or plate bending). eulerBernoulliBeam.js is currently
+   * the only place in this codebase that uses it.
    */
-  getBasisFunctions(ksi, eta = null) {
+  getBasisFunctions(ksi, eta = null, elementLength = null) {
     let basisFunction = [];
     let basisFunctionDerivKsi = [];
     let basisFunctionDerivEta = [];
+    let basisFunctionDerivKsi2 = [];
 
     if (this.meshDimension === "1D") {
       if (this.elementOrder === "linear") {
@@ -57,6 +66,34 @@ export class BasisFunctions {
         basisFunctionDerivKsi[0] = -3 + 4 * ksi;
         basisFunctionDerivKsi[1] = 4 - 8 * ksi;
         basisFunctionDerivKsi[2] = -1 + 4 * ksi;
+      } else if (this.elementOrder === "hermiteCubic") {
+        // Cubic Hermite basis functions for 1D Euler-Bernoulli beam elements
+        // DOF order: [w1, theta1, w2, theta2], with w = transverse deflection and
+        // theta = dw/dx the (physical) slope. The h-scaling on the theta-associated
+        // functions converts the interpolated dw/dksi at the nodes into the actual
+        // rotation DOF, so no extra scaling is required when mapping derivatives to x
+        if (elementLength === null) {
+          errorLog("elementLength is required to evaluate 'hermiteCubic' basis functions");
+          return;
+        }
+        const h = elementLength;
+
+        basisFunction[0] = 2 * ksi ** 3 - 3 * ksi ** 2 + 1;
+        basisFunction[1] = h * (ksi ** 3 - 2 * ksi ** 2 + ksi);
+        basisFunction[2] = -2 * ksi ** 3 + 3 * ksi ** 2;
+        basisFunction[3] = h * (ksi ** 3 - ksi ** 2);
+
+        // First derivatives of basis functions with respect to ksi
+        basisFunctionDerivKsi[0] = 6 * ksi ** 2 - 6 * ksi;
+        basisFunctionDerivKsi[1] = h * (3 * ksi ** 2 - 4 * ksi + 1);
+        basisFunctionDerivKsi[2] = -6 * ksi ** 2 + 6 * ksi;
+        basisFunctionDerivKsi[3] = h * (3 * ksi ** 2 - 2 * ksi);
+
+        // Second derivatives of basis functions with respect to ksi
+        basisFunctionDerivKsi2[0] = 12 * ksi - 6;
+        basisFunctionDerivKsi2[1] = h * (6 * ksi - 4);
+        basisFunctionDerivKsi2[2] = -12 * ksi + 6;
+        basisFunctionDerivKsi2[3] = h * (6 * ksi - 2);
       }
     } else if (this.meshDimension === "2D") {
       if (eta === null) {
@@ -152,6 +189,6 @@ export class BasisFunctions {
       }
     }
 
-    return { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta };
+    return { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta, basisFunctionDerivKsi2 };
   }
 }

src/methods/numericalIntegration.js

@@ -43,6 +43,20 @@ export class NumericalIntegration {
       gaussWeights[0] = 5 / 18;
       gaussWeights[1] = 8 / 18;
       gaussWeights[2] = 5 / 18;
+    } else if (this.elementOrder === "hermiteCubic") {
+      // For cubic Hermite elements, use 6-point Gauss quadrature (exact up to degree 11)
+      gaussPoints[0] = 0.03376524289842399;
+      gaussPoints[1] = 0.16939530676686775;
+      gaussPoints[2] = 0.38069040695840155;
+      gaussPoints[3] = 0.61930959304159845;
+      gaussPoints[4] = 0.8306046932331322;
+      gaussPoints[5] = 0.966234757101576;
+      gaussWeights[0] = 0.08566224618958517;
+      gaussWeights[1] = 0.1803807865240693;
+      gaussWeights[2] = 0.23395696728634552;
+      gaussWeights[3] = 0.23395696728634552;
+      gaussWeights[4] = 0.1803807865240693;
+      gaussWeights[5] = 0.08566224618958517;
     }
 
     return { gaussPoints, gaussWeights };