From 365cb13b07446216c3240fee83a40f47bc692dd8 Mon Sep 17 00:00:00 2001 From: ferrari212 Date: Sat, 18 Jul 2026 10:13:54 +0200 Subject: [PATCH 1/5] Add cubic Hermite basis functions and quadrature rule --- src/mesh/basisFunctions.js | 43 +++++++++++++++++++++++++++-- src/methods/numericalIntegration.js | 14 ++++++++++ 2 files changed, 54 insertions(+), 3 deletions(-) diff --git a/src/mesh/basisFunctions.js b/src/mesh/basisFunctions.js index c5af18f..4fb6c33 100644 --- a/src/mesh/basisFunctions.js +++ b/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 }; } } diff --git a/src/methods/numericalIntegration.js b/src/methods/numericalIntegration.js index 3f05c0a..04def3c 100644 --- a/src/methods/numericalIntegration.js +++ b/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 }; From 066106c51f3a0481545a54f09696ba8420e7d5b0 Mon Sep 17 00:00:00 2001 From: ferrari212 Date: Sat, 18 Jul 2026 10:14:23 +0200 Subject: [PATCH 2/5] Add 1D Euler-Bernoulli beam solver --- src/FEAScript.js | 23 ++++ src/models/beamBoundaryConditions.js | 134 +++++++++++++++++++++ src/models/eulerBernoulliBeam.js | 169 +++++++++++++++++++++++++++ 3 files changed, 326 insertions(+) create mode 100644 src/models/beamBoundaryConditions.js create mode 100644 src/models/eulerBernoulliBeam.js diff --git a/src/FEAScript.js b/src/FEAScript.js index 00ea1e9..966bdee 100644 --- a/src/FEAScript.js +++ b/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"); diff --git a/src/models/beamBoundaryConditions.js b/src/models/beamBoundaryConditions.js new file mode 100644 index 0000000..d879e61 --- /dev/null +++ b/src/models/beamBoundaryConditions.js @@ -0,0 +1,134 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Lightweight Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; + +/** + * Class to handle boundary conditions for the 1D Euler-Bernoulli beam model + * + * Unlike the side-based boundary conditions used by the scalar 1D/2D models + * (which only tag the two ends of a 1D domain as boundary "0" and "1"), + * beam problems routinely carry essential and natural conditions at interior + * nodes as well (e.g. an intermediate roller support, or a point load applied + * partway along the beam). For this reason, `boundaryConditions` here is keyed + * by the 1-based GLOBAL NODE NUMBER rather than by a domain side, and each key + * maps to an ARRAY of condition tuples, since a single node routinely carries + * more than one condition at once (e.g. a spring support plus a point load): + * + * boundaryConditions = { + * "1": [["fixed"]], // clamped: w = 0, theta = 0 + * "2": [["pinned"], ["moment", 1250]], // roller + applied moment + * "3": [["spring", 200], ["force", -2500]], // elastic support + applied force + * }; + * + * Supported condition types: + * - ["fixed"] Essential: w = 0 and theta = 0 (clamped/built-in 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 + * reference deflection uRef (default uRef = 0) + */ +export class BeamBoundaryConditions { + /** + * Constructor to initialize the BeamBoundaryConditions class + * @param {object} boundaryConditions - Object containing boundary conditions, keyed by 1-based node number + * @param {number} totalNodesX - Total number of nodes along the beam + * @param {number} dofsPerNode - Number of degrees of freedom per node (2: deflection, rotation) + */ + constructor(boundaryConditions, totalNodesX, dofsPerNode) { + this.boundaryConditions = boundaryConditions; + this.totalNodesX = totalNodesX; + this.dofsPerNode = dofsPerNode; + } + + /** + * Function to impose natural (point force/moment) and spring (Robin-type) boundary conditions + * This must be called BEFORE imposeEssentialBoundaryConditions(), since essential conditions + * override (zero out) the equation row at a constrained DOF, which must take precedence over + * any natural/spring contribution assembled at the same DOF + * @param {array} residualVector - The residual (load) vector to be modified + * @param {array} jacobianMatrix - The Jacobian (stiffness) matrix to be modified + */ + imposeNaturalAndSpringBoundaryConditions(residualVector, jacobianMatrix) { + Object.keys(this.boundaryConditions).forEach((nodeKey) => { + const globalNodeIndex = Number(nodeKey) - 1; // Convert 1-based node number to 0-based index + const deflectionDOF = this.dofsPerNode * globalNodeIndex; + const rotationDOF = deflectionDOF + 1; + + this.boundaryConditions[nodeKey].forEach((condition) => { + const [conditionType, value1, value2] = condition; + + if (conditionType === "force") { + residualVector[deflectionDOF] += value1; + debugLog(`Node ${nodeKey}: Applied point force ${value1} (natural BC)`); + } else if (conditionType === "moment") { + residualVector[rotationDOF] += value1; + debugLog(`Node ${nodeKey}: Applied point moment ${value1} (natural BC)`); + } else if (conditionType === "spring") { + const springConstant = value1; + const referenceDeflection = value2 ?? 0; + jacobianMatrix[deflectionDOF][deflectionDOF] += springConstant; + residualVector[deflectionDOF] += springConstant * referenceDeflection; + debugLog( + `Node ${nodeKey}: Applied transverse elastic spring, k=${springConstant} (mixed/Robin BC)`, + ); + } + }); + }); + } + + /** + * Function to impose essential (deflection/rotation) boundary conditions (Dirichlet-type) + * This must be called AFTER imposeNaturalAndSpringBoundaryConditions() + * @param {array} residualVector - The residual vector to be modified + * @param {array} jacobianMatrix - The Jacobian matrix to be modified + */ + imposeEssentialBoundaryConditions(residualVector, jacobianMatrix) { + const totalDOFs = residualVector.length; + + const applyDirichlet = (dofIndex, prescribedValue) => { + residualVector[dofIndex] = prescribedValue; + for (let colIndex = 0; colIndex < totalDOFs; colIndex++) { + jacobianMatrix[dofIndex][colIndex] = 0; + } + jacobianMatrix[dofIndex][dofIndex] = 1; + }; + + Object.keys(this.boundaryConditions).forEach((nodeKey) => { + const globalNodeIndex = Number(nodeKey) - 1; // Convert 1-based node number to 0-based index + const deflectionDOF = this.dofsPerNode * globalNodeIndex; + const rotationDOF = deflectionDOF + 1; + + this.boundaryConditions[nodeKey].forEach((condition) => { + const [conditionType, value] = condition; + + if (conditionType === "fixed") { + applyDirichlet(deflectionDOF, 0); + applyDirichlet(rotationDOF, 0); + debugLog(`Node ${nodeKey}: Applied fixed (clamped) support: w=0, theta=0 (essential BC)`); + } else if (conditionType === "pinned" || conditionType === "deflection") { + applyDirichlet(deflectionDOF, value ?? 0); + debugLog(`Node ${nodeKey}: Applied deflection w=${value ?? 0} (essential BC)`); + } else if (conditionType === "rotationFixed" || conditionType === "rotation") { + applyDirichlet(rotationDOF, value ?? 0); + debugLog(`Node ${nodeKey}: Applied rotation theta=${value ?? 0} (essential BC)`); + } else if ( + conditionType !== "force" && + conditionType !== "moment" && + conditionType !== "spring" + ) { + errorLog(`Unknown beam boundary condition type: "${conditionType}"`); + } + }); + }); + } +} diff --git a/src/models/eulerBernoulliBeam.js b/src/models/eulerBernoulliBeam.js new file mode 100644 index 0000000..c061a98 --- /dev/null +++ b/src/models/eulerBernoulliBeam.js @@ -0,0 +1,169 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Lightweight Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { BasisFunctions } from "../mesh/basisFunctions.js"; +import { NumericalIntegration } from "../methods/numericalIntegration.js"; +import { BeamBoundaryConditions } from "./beamBoundaryConditions.js"; +import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; + +/** + * Function to assemble the Jacobian matrix and residual vector for the 1D Euler-Bernoulli beam model + * + * The beam is governed by the fourth-order equation: + * d²/dx²( EI(x) d²w/dx² ) + c0(x) w = q(x) + * with weak form (element level): + * K_ij = ∫ [ EI(x) psi_i'' psi_j'' + c0(x) psi_i psi_j ] dx + * F_i = ∫ q(x) psi_i dx + (point forces/moments/spring contributions from boundary conditions) + * where psi_i are cubic Hermite shape functions (see basisFunctions.js) and c0 is an optional + * elastic-foundation modulus. + * + * Each mesh node carries 2 degrees of freedom, [w, theta] with theta = dw/dx, ordered + * consecutively per node: [w_0, theta_0, w_1, theta_1, ...]. The mesh geometry itself only + * needs the 2 end nodes of each element (i.e. meshConfig.elementOrder must be 'linear'); the + * richer cubic Hermite interpolation is applied internally to the field, independent of the + * geometric element order (a "subparametric" formulation, standard for beam elements). + * + * @param {object} meshData - Object containing prepared mesh data + * @param {object} boundaryConditions - Object containing boundary conditions, keyed by 1-based node + * number (see beamBoundaryConditions.js for the supported condition types) + * @param {object} coefficientFunctions - Functions of x for the beam model: + * - EI(x): bending stiffness (Young's modulus times second moment of area) + * - c0(x): elastic foundation modulus (optional, defaults to 0) + * - q(x): distributed transverse load (optional, defaults to 0) + * @returns {object} An object containing: + * - jacobianMatrix: The assembled Jacobian (stiffness) matrix + * - residualVector: The assembled residual (load) vector + * - dofsPerNode: Number of degrees of freedom per node (2) + * - totalDOFs: Total number of degrees of freedom in the assembled system + * + * For consistency across both linear and nonlinear formulations, + * this project always refers to the assembled right-hand side vector + * as `residualVector` and the assembled system matrix as `jacobianMatrix`. + * + * In linear problems `jacobianMatrix` is equivalent to the + * classic stiffness matrix and `residualVector` corresponds to the traditional load (RHS) vector. + */ +export function assembleEulerBernoulliBeamMat(meshData, boundaryConditions, coefficientFunctions) { + basicLog("Starting Euler-Bernoulli beam matrix assembly..."); + + // Extract mesh data + const { nodesXCoordinates, nop, totalElements, totalNodesX, meshDimension, elementOrder } = meshData; + + if (meshDimension !== "1D") { + errorLog("Euler-Bernoulli beam solver requires a 1D mesh"); + } + if (elementOrder !== "linear") { + errorLog( + "Euler-Bernoulli beam solver requires 'linear' (2-node) elements for the beam geometry; " + + "the cubic Hermite field interpolation for w(x) is applied internally regardless of this setting", + ); + } + + // Extract coefficient functions, with sensible defaults for the optional ones + const { EI, c0 = () => 0, q = () => 0 } = coefficientFunctions; + + const dofsPerNode = 2; // [deflection w, rotation theta = dw/dx] at each node + const totalDOFs = dofsPerNode * totalNodesX; + + // Initialize global Jacobian matrix and residual vector + let residualVector = new Array(totalDOFs).fill(0); + let jacobianMatrix = []; + for (let dofIndex = 0; dofIndex < totalDOFs; dofIndex++) { + jacobianMatrix.push(new Array(totalDOFs).fill(0)); + } + + // Cubic Hermite basis functions for the field, with a 4-point Gauss quadrature rule + const basisFunctions = new BasisFunctions({ meshDimension: "1D", elementOrder: "hermiteCubic" }); + const numericalIntegration = new NumericalIntegration({ meshDimension: "1D", elementOrder: "hermiteCubic" }); + const { gaussPoints, gaussWeights } = numericalIntegration.getGaussPointsAndWeights(); + + // Matrix assembly + for (let elementIndex = 0; elementIndex < totalElements; elementIndex++) { + // Beam elements only use their 2 end (geometric) nodes + const globalNode1 = nop[elementIndex][0] - 1; // Convert to 0-based indexing + const globalNode2 = nop[elementIndex][1] - 1; + const x1 = nodesXCoordinates[globalNode1]; + const x2 = nodesXCoordinates[globalNode2]; + const elementLength = x2 - x1; + + // Map local DOFs [w1, theta1, w2, theta2] to global DOF indices + const dofMap = [ + dofsPerNode * globalNode1, + dofsPerNode * globalNode1 + 1, + dofsPerNode * globalNode2, + dofsPerNode * globalNode2 + 1, + ]; + + // Loop over Gauss points + for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { + const ksi = gaussPoints[gaussPointIndex]; + + // Get Hermite cubic basis functions and their ksi-derivatives at this Gauss point + const { basisFunction, basisFunctionDerivKsi2 } = basisFunctions.getBasisFunctions( + ksi, + null, + elementLength, + ); + + // Straight 2-node element: the Jacobian dx/dksi is the constant element length + const detJacobian = elementLength; + + // Second derivative with respect to x: d2N/dx2 = (1/J^2) d2N/dksi2 (valid since J is + // constant along a straight element) + const basisFunctionDerivXX = basisFunctionDerivKsi2.map( + (secondDerivKsi) => secondDerivKsi / detJacobian ** 2, + ); + + // Physical coordinate of this Gauss point (linear mapping between the 2 end nodes) + const xCoord = x1 + elementLength * ksi; + + // Evaluate coefficient functions at this physical coordinate + const EIVal = EI(xCoord); + const c0Val = c0(xCoord); + const qVal = q(xCoord); + + const weightFactor = gaussWeights[gaussPointIndex] * detJacobian; + + // Computation of the element contributions to the residual vector and Jacobian matrix + for (let localIndex1 = 0; localIndex1 < 4; localIndex1++) { + const globalDOF1 = dofMap[localIndex1]; + + // Distributed load contribution to the residual (load) vector + residualVector[globalDOF1] += weightFactor * qVal * basisFunction[localIndex1]; + + for (let localIndex2 = 0; localIndex2 < 4; localIndex2++) { + const globalDOF2 = dofMap[localIndex2]; + + // Bending stiffness (curvature-curvature) and elastic-foundation terms + jacobianMatrix[globalDOF1][globalDOF2] += + weightFactor * + (EIVal * basisFunctionDerivXX[localIndex1] * basisFunctionDerivXX[localIndex2] + + c0Val * basisFunction[localIndex1] * basisFunction[localIndex2]); + } + } + } + } + + // Apply boundary conditions: natural/spring contributions first, essential conditions last + // (essential conditions override the equation row, taking precedence over any natural/spring + // contribution assembled at the same DOF) + const beamBoundaryConditions = new BeamBoundaryConditions(boundaryConditions, totalNodesX, dofsPerNode); + beamBoundaryConditions.imposeNaturalAndSpringBoundaryConditions(residualVector, jacobianMatrix); + beamBoundaryConditions.imposeEssentialBoundaryConditions(residualVector, jacobianMatrix); + + basicLog("Euler-Bernoulli beam matrix assembly completed"); + + return { + jacobianMatrix, + residualVector, + dofsPerNode, + totalDOFs, + }; +} From 089aa4b08b6a4908c5785388fd214e1ab5d36529 Mon Sep 17 00:00:00 2001 From: ferrari212 Date: Sat, 18 Jul 2026 10:14:57 +0200 Subject: [PATCH 3/5] Add clamped/roller/spring beam example --- examples/Beam1DFEM/Beam1DEuler_Bernoulli.js | 95 ++++++++++++++++++ examples/Beam1DFEM/README.md | 102 ++++++++++++++++++++ 2 files changed, 197 insertions(+) create mode 100644 examples/Beam1DFEM/Beam1DEuler_Bernoulli.js create mode 100644 examples/Beam1DFEM/README.md diff --git a/examples/Beam1DFEM/Beam1DEuler_Bernoulli.js b/examples/Beam1DFEM/Beam1DEuler_Bernoulli.js new file mode 100644 index 0000000..ef22f55 --- /dev/null +++ b/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)}`, + ); +} diff --git a/examples/Beam1DFEM/README.md b/examples/Beam1DFEM/README.md new file mode 100644 index 0000000..923e53c --- /dev/null +++ b/examples/Beam1DFEM/README.md @@ -0,0 +1,102 @@ +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 +``` From 589c1b5102974d950675f79d7fa3c72228841428 Mon Sep 17 00:00:00 2001 From: ferrari212 Date: Sat, 18 Jul 2026 10:15:20 +0200 Subject: [PATCH 4/5] Add tests for the beam solver --- .../EulerBernoulliBeam/REGRESSION.md | 92 ++++++++++ .../EulerBernoulliBeam/regression.test.js | 159 ++++++++++++++++++ tests/unit/eulerBernoulliBeam.test.js | 121 +++++++++++++ 3 files changed, 372 insertions(+) create mode 100644 tests/regression/EulerBernoulliBeam/REGRESSION.md create mode 100644 tests/regression/EulerBernoulliBeam/regression.test.js create mode 100644 tests/unit/eulerBernoulliBeam.test.js diff --git a/tests/regression/EulerBernoulliBeam/REGRESSION.md b/tests/regression/EulerBernoulliBeam/REGRESSION.md new file mode 100644 index 0000000..f78eed9 --- /dev/null +++ b/tests/regression/EulerBernoulliBeam/REGRESSION.md @@ -0,0 +1,92 @@ +# Regression Test — EulerBernoulliBeam + +## Purpose + +This test guards the numerical output of the 1D Euler-Bernoulli beam example against +unintended changes to the beam solver, assembler, or mesh-generation logic. + +It replicates exactly the problem set up in +[`Beam1DEuler_Bernoulli.js`](../../../examples/Beam1DFEM/Beam1DEuler_Bernoulli.js) — 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) — and asserts both a set of +known-good baseline values and, independently of those baseline numbers, that the resulting +finite element solution satisfies global static equilibrium. + +## Problem setup + +| Parameter | Value | +| ----------------------------- | -------------------------------------------------------- | +| Domain | 1D beam, 0 – 10 m | +| Mesh | 2 cubic Hermite beam elements of 5 m each (3 nodes) | +| Bending stiffness EI | 2.0 × 10⁶ N·m² (constant) | +| Distributed load | −1,000 N/m over 0 ≤ x ≤ 5 m only | +| Node 1 (x = 0) | Fixed (clamped): w = 0, theta = 0 | +| Node 2 (x = 5) | Pinned (roller): w = 0, plus an applied moment M = 1,250 N·m | +| Node 3 (x = 10) | Transverse spring k = 200 N/m, plus a point load P = −2,500 N | +| Solver | LU decomposition (`lusolve`) | + +## Expected values + +| Quantity | Value | +| ------------------------ | ------------------------- | +| w₁ (deflection, node 1) | 0 m | +| θ₁ (rotation, node 1) | 0 rad | +| w₂ (deflection, node 2) | 0 m | +| θ₂ (rotation, node 2) | −5.6790761806 × 10⁻³ rad | +| w₃ (deflection, node 3) | −8.0144777663 × 10⁻² m | +| θ₃ (rotation, node 3) | −2.1203895209 × 10⁻² rad | + +Tolerance used in the baseline assertions: `1e-8`. + +These baseline values were derived from the FEAScript implementation itself (the source +material only publishes the `.inp`-style problem setup, not the FEM1D program's numeric +output), and were independently cross-checked three ways before being adopted as the +regression baseline: + +1. The single-element bending stiffness matrix (for constant EI) was verified against the + closed-form Hermite beam element stiffness matrix found in any FEM textbook — see + `tests/unit/eulerBernoulliBeam.test.js`. +2. A simple cantilever case (fixed-free beam under a tip point load) was verified against the + classical closed-form solution `w_tip = P·L³/(3EI)`, `theta_tip = P·L²/(2EI)` — also in + `tests/unit/eulerBernoulliBeam.test.js`. +3. The full statically-indeterminate solution above was checked against global force and + moment equilibrium (see below) — this check is also asserted every run, independently of + the baseline numbers. + +## Equilibrium check + +In addition to the baseline values, this test recomputes the raw (no boundary conditions +applied) element assembly, uses it to recover the support reactions and spring force from the +solved DOF values, and asserts: + +- **Force equilibrium**: sum of the clamped reaction, roller reaction, spring force, distributed + load resultant, and applied point load is ~0. +- **Moment equilibrium** about x = 0: sum of the clamped reaction moment, roller reaction moment + arm, spring force moment arm, distributed load resultant moment arm, point load moment arm, + and applied moment is ~0. + +Tolerance used for the equilibrium assertions: `1e-6`. + +## How to run + +From the repository root: + +```bash +node tests/regression/EulerBernoulliBeam/regression.test.js +``` + +The `test` script in `package.json` also runs this file, so `npm test` works too. + +## After modifying the code + +| Situation | Action | +| ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------ | +| Bug fix that should not change results | Run the test — it must still pass. | +| Intentional algorithm change (new integration rule, new element type, etc.) | Re-derive the expected values, update `EXPECTED` in `regression.test.js`, and document the reason here. | +| New boundary condition type | Update both the test and `Beam1DEuler_Bernoulli.js` together. | + +## Change log + +| Date | Change | New expected values | +| ---------- | ---------------------------- | -------------------- | +| 2026-07-17 | Initial regression baseline | See table above | diff --git a/tests/regression/EulerBernoulliBeam/regression.test.js b/tests/regression/EulerBernoulliBeam/regression.test.js new file mode 100644 index 0000000..a228fdf --- /dev/null +++ b/tests/regression/EulerBernoulliBeam/regression.test.js @@ -0,0 +1,159 @@ +/** + * Regression test for the 1D Euler-Bernoulli beam model (EulerBernoulliBeam) + * + * Replicates the exact setup from Beam1DEuler_Bernoulli.js — 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) — and asserts: + * 1) the deflection/rotation solution vector against known-good baseline values + * 2) global force and moment equilibrium of the resulting FE solution, which + * holds regardless of the specific numeric baseline and independently + * confirms the assembled system is physically consistent + * + * Run: node tests/regression/EulerBernoulliBeam/regression.test.js + */ + +import * as mathjs from "mathjs"; +import { FEAScriptModel } from "../../../src/FEAScript.js"; +import { assembleEulerBernoulliBeamMat } from "../../../src/models/eulerBernoulliBeam.js"; +import { basicLog, errorLog } from "../../../src/utilities/logging.js"; + +// FEAScript.js references `math` as a global (loaded via CDN in browser). +// Set it here before any solve() call. +globalThis.math = mathjs; + +// Baseline values (see REGRESSION.md for derivation and verification) +const EXPECTED = { + w1: 0, + theta1: 0, + w2: 0, + theta2: -0.0056790761806273645, + w3: -0.08014477766287427, + theta3: -0.02120389520854876, +}; +const TOLERANCE = 1e-8; + +const EI = 2.0e6; +const springConstant = 200; +const distributedLoad = -1000; +const distributedLoadSpan = [0, 5]; +const appliedMoment = 1250; +const appliedForce = -2500; + +function runSimulation() { + const model = new FEAScriptModel(); + + model.setModelConfig("eulerBernoulliBeamScript", { + coefficientFunctions: { + EI: (x) => EI, + q: (x) => (x <= distributedLoadSpan[1] ? distributedLoad : 0), + }, + }); + model.setMeshConfig({ + meshDimension: "1D", + elementOrder: "linear", + numElementsX: 2, + maxX: 10, + }); + + model.addBoundaryCondition("1", [["fixed"]]); + model.addBoundaryCondition("2", [["pinned"], ["moment", appliedMoment]]); + model.addBoundaryCondition("3", [["spring", springConstant], ["force", appliedForce]]); + model.setSolverMethod("lusolve"); + + return model.solve(); +} + +function assert(condition, message) { + if (!condition) { + errorLog(`FAIL: ${message}`); + process.exit(1); + } + basicLog(`PASS: ${message}`); +} + +basicLog(""); +basicLog("================================"); +basicLog("Starting regression test for the 1D Euler-Bernoulli beam model..."); + +const { solutionVector } = runSimulation(); +// mathjs' lusolve returns a nested array: [[v0], [v1], ...] +const flatSolution = solutionVector.map((entry) => (Array.isArray(entry) ? entry[0] : entry)); +const [w1, theta1, w2, theta2, w3, theta3] = flatSolution; + +// --------------------------------------------------------------------------- +// 1) Baseline values +// --------------------------------------------------------------------------- +assert(Math.abs(w1 - EXPECTED.w1) < TOLERANCE, `w1: expected ~${EXPECTED.w1}, got ${w1}`); +assert(Math.abs(theta1 - EXPECTED.theta1) < TOLERANCE, `theta1: expected ~${EXPECTED.theta1}, got ${theta1}`); +assert(Math.abs(w2 - EXPECTED.w2) < TOLERANCE, `w2: expected ~${EXPECTED.w2}, got ${w2}`); +assert( + Math.abs(theta2 - EXPECTED.theta2) < TOLERANCE, + `theta2: expected ${EXPECTED.theta2}, got ${theta2}`, +); +assert(Math.abs(w3 - EXPECTED.w3) < TOLERANCE, `w3: expected ${EXPECTED.w3}, got ${w3}`); +assert( + Math.abs(theta3 - EXPECTED.theta3) < TOLERANCE, + `theta3: expected ${EXPECTED.theta3}, got ${theta3}`, +); + +// --------------------------------------------------------------------------- +// 2) Global equilibrium check, independent of the specific baseline values above. +// Recompute the raw (no boundary conditions applied) element assembly and use +// it to recover support reactions: at any DOF, K_raw.u - F_raw equals whatever +// external generalized force (reaction, spring force, or applied load) is +// required to balance that row. +// --------------------------------------------------------------------------- +const meshDataForCheck = { + nodesXCoordinates: [0, 5, 10], + nop: [ + [1, 2], + [2, 3], + ], + totalElements: 2, + totalNodesX: 3, + meshDimension: "1D", + elementOrder: "linear", +}; +const { jacobianMatrix: Kraw, residualVector: Fraw } = assembleEulerBernoulliBeamMat( + meshDataForCheck, + {}, + { + EI: (x) => EI, + q: (x) => (x <= distributedLoadSpan[1] ? distributedLoad : 0), + }, +); + +function matVec(matrix, vector) { + return matrix.map((row) => row.reduce((sum, value, colIndex) => sum + value * vector[colIndex], 0)); +} +const residualAtEachDOF = matVec(Kraw, flatSolution).map((value, i) => value - Fraw[i]); + +const reactionForceNode1 = residualAtEachDOF[0]; +const reactionMomentNode1 = residualAtEachDOF[1]; +const reactionForceNode2 = residualAtEachDOF[2]; +const springForceOnBeam = -springConstant * w3; + +const totalDistributedLoad = distributedLoad * (distributedLoadSpan[1] - distributedLoadSpan[0]); +const distributedLoadCentroidX = (distributedLoadSpan[0] + distributedLoadSpan[1]) / 2; + +const sumVerticalForces = + reactionForceNode1 + reactionForceNode2 + springForceOnBeam + totalDistributedLoad + appliedForce; +const sumMomentsAboutOrigin = + reactionMomentNode1 + + reactionForceNode2 * 5 + + springForceOnBeam * 10 + + totalDistributedLoad * distributedLoadCentroidX + + appliedForce * 10 + + appliedMoment; + +const equilibriumTolerance = 1e-6; +assert( + Math.abs(sumVerticalForces) < equilibriumTolerance, + `Global force equilibrium holds (sum=${sumVerticalForces})`, +); +assert( + Math.abs(sumMomentsAboutOrigin) < equilibriumTolerance, + `Global moment equilibrium about x=0 holds (sum=${sumMomentsAboutOrigin})`, +); + +basicLog("================================"); diff --git a/tests/unit/eulerBernoulliBeam.test.js b/tests/unit/eulerBernoulliBeam.test.js new file mode 100644 index 0000000..8b71882 --- /dev/null +++ b/tests/unit/eulerBernoulliBeam.test.js @@ -0,0 +1,121 @@ +/** + * Unit tests for the 1D Euler-Bernoulli beam model (assembleEulerBernoulliBeamMat) + * + * Covers: + * - Single-element bending stiffness matrix vs. the closed-form textbook + * Hermite beam element stiffness matrix (2 independent EI/L combinations) + * - Cantilever beam under a tip point load vs. the classical closed-form + * solution (w_tip = P*L^3/(3EI), theta_tip = P*L^2/(2EI)) + * + * Run: node tests/unit/eulerBernoulliBeam.test.js + */ + +import * as mathjs from "mathjs"; +globalThis.math = mathjs; + +import { assembleEulerBernoulliBeamMat } from "../../src/models/eulerBernoulliBeam.js"; +import { prepareMesh } from "../../src/mesh/meshUtils.js"; +import { solveLinearSystem } from "../../src/methods/linearSystemSolver.js"; +import { basicLog, errorLog } from "../../src/utilities/logging.js"; + +basicLog(""); +basicLog("================================"); +basicLog("Unit tests: assembleEulerBernoulliBeamMat"); + +let passed = 0; +let failed = 0; + +function assert(condition, message) { + if (!condition) { + errorLog(`FAIL: ${message}`); + failed++; + } else { + basicLog(`PASS: ${message}`); + passed++; + } +} + +/** + * Closed-form Euler-Bernoulli beam element stiffness matrix for constant EI, + * DOF order [w1, theta1, w2, theta2] (see e.g. any FEM textbook) + */ +function closedFormBeamStiffness(EI, L) { + const c = EI / L ** 3; + return [ + [12 * c, 6 * L * c, -12 * c, 6 * L * c], + [6 * L * c, 4 * L ** 2 * c, -6 * L * c, 2 * L ** 2 * c], + [-12 * c, -6 * L * c, 12 * c, -6 * L * c], + [6 * L * c, 2 * L ** 2 * c, -6 * L * c, 4 * L ** 2 * c], + ]; +} + +function maxAbsDiff(A, B) { + let diff = 0; + for (let i = 0; i < A.length; i++) { + for (let j = 0; j < A[i].length; j++) { + diff = Math.max(diff, Math.abs(A[i][j] - B[i][j])); + } + } + return diff; +} + +// --------------------------------------------------------------------------- +// 1. Single-element stiffness matrix vs. closed-form matrix +// --------------------------------------------------------------------------- +basicLog(""); +basicLog("[1] Single-element stiffness matrix vs. closed-form Hermite beam matrix"); + +for (const [EI, L] of [ + [2.0e6, 5], + [7.5e4, 2.5], +]) { + const meshData = prepareMesh({ meshDimension: "1D", elementOrder: "linear", numElementsX: 1, maxX: L }); + const { jacobianMatrix } = assembleEulerBernoulliBeamMat(meshData, {}, { EI: () => EI }); + const diff = maxAbsDiff(jacobianMatrix, closedFormBeamStiffness(EI, L)); + assert(diff < 1e-6, `Element stiffness matches closed-form matrix for EI=${EI}, L=${L} (diff=${diff})`); +} + +// --------------------------------------------------------------------------- +// 2. Cantilever beam under a tip point load vs. classical closed-form solution +// --------------------------------------------------------------------------- +basicLog(""); +basicLog("[2] Cantilever beam under a tip point load"); + +const L = 4; +const EI = 1.0e5; +const P = -1000; // Downward tip load + +const meshData = prepareMesh({ meshDimension: "1D", elementOrder: "linear", numElementsX: 1, maxX: L }); +const { jacobianMatrix, residualVector } = assembleEulerBernoulliBeamMat( + meshData, + { 1: [["fixed"]], 2: [["force", P]] }, + { EI: () => EI }, +); +const { solutionVector } = solveLinearSystem("lusolve", jacobianMatrix, residualVector); +const flatSolution = solutionVector.map((entry) => (Array.isArray(entry) ? entry[0] : entry)); + +const wExact = (P * L ** 3) / (3 * EI); +const thetaExact = (P * L ** 2) / (2 * EI); +const tolerance = 1e-9; + +assert( + Math.abs(flatSolution[0]) < tolerance && Math.abs(flatSolution[1]) < tolerance, + "Clamped end has zero deflection and rotation", +); +assert( + Math.abs(flatSolution[2] - wExact) < tolerance, + `Tip deflection matches closed form (got ${flatSolution[2]}, expected ${wExact})`, +); +assert( + Math.abs(flatSolution[3] - thetaExact) < tolerance, + `Tip rotation matches closed form (got ${flatSolution[3]}, expected ${thetaExact})`, +); + +// --------------------------------------------------------------------------- +basicLog(""); +if (failed > 0) { + errorLog(`${passed} passed, ${failed} failed.`); + process.exit(1); +} else { + basicLog(`${passed} passed, ${failed} failed.`); +} From f8158a136185606ccec0cfade88130c9119bc584 Mon Sep 17 00:00:00 2001 From: ferrari212 Date: Sat, 18 Jul 2026 10:23:09 +0200 Subject: [PATCH 5/5] Improving the general form description --- src/models/generalFormPDE.js | 63 ++++++++++++++++++++++++++++++------ 1 file changed, 53 insertions(+), 10 deletions(-) diff --git a/src/models/generalFormPDE.js b/src/models/generalFormPDE.js index dce8538..fc25ec3 100644 --- a/src/models/generalFormPDE.js +++ b/src/models/generalFormPDE.js @@ -13,13 +13,40 @@ import { GenericBoundaryConditions } from "./genericBoundaryConditions.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; /** - * Function to assemble the Jacobian matrix and residuals vector for the general form PDE model - * @param {object} meshData - Object containing prepared mesh data - * @param {object} boundaryConditions - Object containing boundary conditions - * @param {object} coefficientFunctions - Functions A(x), B(x), C(x), D(x) for the PDE - * @returns {object} An object containing: - * - jacobianMatrix: The assembled Jacobian matrix - * - residualVector: The assembled residual vector + * Assembles the discrete weak form of a general 1D scalar linear PDE of the form + * + * -d/dx (A(x) du/dx) + B(x) du/dx + C(x) u = D(x) + * + * where: + * - u(x) is the scalar unknown field being solved for, + * - A(x) is the diffusion or stiffness coefficient, + * - B(x) is the advection or transport coefficient, + * - C(x) is the reaction or damping coefficient, + * - D(x) is the source or forcing term. + * + * This is a general-purpose finite-element model rather than a fluid-only + * formulation. In practice, the same structure can represent diffusion, + * transport, reaction, damping, or stiffness-dominated scalar field problems. + * + * In the finite-element implementation, the residual and Jacobian are assembled + * from the Galerkin weak form at each Gauss point. The code evaluates the + * coefficient functions A(x), B(x), C(x), and D(x) at the mapped physical + * coordinate x and accumulates the corresponding diffusion, advection, + * reaction, and source terms into the global system. + * + * The resulting algebraic system is: + * + * K(u) = F + * + * where K is the assembled Jacobian matrix and F is the residual forcing term. + * + * @param {object} meshData - Prepared mesh data containing nodal coordinates, + * element connectivity, element order, and dimension. + * @param {object} boundaryConditions - Boundary condition definition object. + * @param {object} coefficientFunctions - Functions A(x), B(x), C(x), D(x) + * describing the PDE coefficients. + * @returns {{ jacobianMatrix: number[][], residualVector: number[] }} + * The assembled global Jacobian matrix and residual vector. */ export function assembleGeneralFormPDEMat(meshData, boundaryConditions, coefficientFunctions) { basicLog("Starting general form PDE matrix assembly..."); @@ -153,9 +180,25 @@ export function assembleGeneralFormPDEMat(meshData, boundaryConditions, coeffici } /** - * Function to assemble the frontal solver matrix for the general form PDE model - * @param {object} data - Object containing element data for the frontal solver - * @returns {object} An object containing local Jacobian matrix and residual vector + * Assembles the local element contribution for the same general 1D PDE model + * used by the frontal solver. The element residual and local Jacobian are + * formed by evaluating the same diffusion, advection, and reaction terms at + * each Gauss point, but only for the current element rather than the full mesh. + * + * The returned local operators are later assembled into the global system by + * the frontal solver. This allows the same PDE structure to be used in both + * the full global assembly path and the element-by-element frontal assembly path. + * + * @param {object} data - Element assembly input for the frontal solver, containing: + * - elementIndex: index of the current element in the mesh, + * - nop: element connectivity array, + * - meshData: mesh coordinates and dimensionality metadata, + * - basisFunctions: interpolation and derivative providers for the element, + * - FEAData: numerical integration data such as Gauss points, weights, + * and nodes-per-element, + * - coefficientFunctions: callback functions A(x), B(x), C(x), and D(x). + * @returns {{ localJacobianMatrix: number[][], localResidualVector: number[], ngl: number[] }} + * Local Jacobian matrix, local residual vector, and global node list. */ export function assembleGeneralFormPDEFront({ elementIndex,