Feature/universal frontal solver

.prettierignore

@@ -0,0 +1,5 @@
+# Build output
+dist/
+
+# Third-party bundles
+src/vendor/

.prettierrc

@@ -0,0 +1,4 @@
+{
+  "printWidth": 110,
+  "trailingComma": "all"
+}

src/FEAScript.js

@@ -12,11 +12,13 @@ import { newtonRaphson } from "./solvers/newtonRaphson.js";
 import { solveLinearSystem } from "./solvers/linearSystemSolver.js";
 import { solveLinearSystemAsync } from "./solvers/linearSystemSolver.js";
 import { prepareMesh } from "./mesh/meshUtils.js";
-import { assembleFrontPropagationMat } from "./models/frontPropagation.js";
-import { assembleGeneralFormPDEMat, assembleGeneralFormPDENonlinearMat } from "./models/generalFormPDE.js";
-import { assembleHeatConductionMat, assembleHeatConductionFront } from "./models/heatConduction.js";
-import { assembleCreepingFlowMatrix } from "./models/creepingFlow.js";
-import { assembleEulerBernoulliBeamMat } from "./models/eulerBernoulliBeam.js";
+import { assembleFrontPropagationMat, assembleFrontPropagationFront } from "./physics/frontPropagation.js";
+import { assembleGeneralFormPDEMat, assembleGeneralFormPDENonlinearMat } from "./physics/generalFormPDE.js";
+import { assembleHeatConductionMat, assembleHeatConductionFront } from "./physics/heatConduction.js";
+import { ThermalBoundaryConditions } from "./physics/thermalBoundaryConditions.js";
+import { GenericBoundaryConditions } from "./physics/genericBoundaryConditions.js";
+import { assembleCreepingFlowMatrix } from "./physics/creepingFlow.js";
+import { assembleEulerBernoulliBeamMat } from "./physics/eulerBernoulliBeam.js";
 import { runFrontalSolver } from "./solvers/frontalSolver.js";
 import { basicLog, debugLog, warnLog, errorLog } from "./utilities/logging.js";
 
@@ -141,11 +143,22 @@ export class FEAScriptModel {
     if (this.solverConfig === "heatConductionScript") {
       // Check if using frontal solver
       if (this.solverMethod === "frontal") {
+        const thermalBoundaryConditions = new ThermalBoundaryConditions(
+          this.boundaryConditions,
+          meshData.boundaryElements,
+          meshData.nop,
+          meshData.meshDimension,
+          meshData.elementOrder,
+          meshData.elementShape,
+        );
         const frontalResult = runFrontalSolver(
-          assembleHeatConductionFront,
+          (elementArgs) =>
+            assembleHeatConductionFront({
+              ...elementArgs,
+              coefficientFunctions: this.coefficientFunctions,
+            }),
           meshData,
-          this.boundaryConditions,
-          { coefficientFunctions: this.coefficientFunctions },
+          thermalBoundaryConditions,
         );
         solutionVector = frontalResult.solutionVector;
       } else {
@@ -176,6 +189,20 @@ export class FEAScriptModel {
         // TODO: Consider using different maxIterations/tolerance for Newton-Raphson and linear solver
         maxIterations: options.maxIterations ?? this.maxIterations,
         tolerance: options.tolerance ?? this.tolerance,
+        // Used only by the frontal solver
+        assembleFront: (elementArgs) =>
+          assembleFrontPropagationFront({
+            ...elementArgs,
+            eikonalActivationFlag: context.eikonalActivationFlag,
+          }),
+        boundaryConditionsHandler: new GenericBoundaryConditions(
+          this.boundaryConditions,
+          meshData.boundaryElements,
+          meshData.nop,
+          meshData.meshDimension,
+          meshData.elementOrder,
+          meshData.elementShape,
+        ),
       };
 
       // Count the incremental steps with an integer, so that the last step activates the eikonal equation

src/mesh/basisFunctions.js → src/elements/basisFunctions.js

@@ -28,6 +28,53 @@ export const triangleBoundarySides = {
   },
 };
 
+/**
+ * Local node indices on each side of the 1D reference element
+ */
+const lineBoundarySides = {
+  linear: {
+    0: [0], // Node at the left side of the reference element
+    1: [1], // Node at the right side of the reference element
+  },
+  quadratic: {
+    0: [0], // Node at the left side of the reference element
+    1: [2], // Node at the right side of the reference element
+  },
+};
+
+/**
+ * Local node indices on each side of the reference quadrilateral
+ * (side 0: bottom, eta = 0; side 1: left, ksi = 0; side 2: top, eta = 1; side 3: right, ksi = 1)
+ * Nodes are listed in order along each side
+ */
+const quadBoundarySides = {
+  linear: {
+    0: [0, 2], // Nodes at the bottom side of the reference element
+    1: [0, 1], // Nodes at the left side of the reference element
+    2: [1, 3], // Nodes at the top side of the reference element
+    3: [2, 3], // Nodes at the right side of the reference element
+  },
+  quadratic: {
+    0: [0, 3, 6], // Nodes at the bottom side of the reference element
+    1: [0, 1, 2], // Nodes at the left side of the reference element
+    2: [2, 5, 8], // Nodes at the top side of the reference element
+    3: [6, 7, 8], // Nodes at the right side of the reference element
+  },
+};
+
+/**
+ * Function to get the local node indices on each side of the reference element
+ * @param {string} meshDimension - The dimension of the mesh, either '1D' or '2D'
+ * @param {string} elementOrder - The order of elements, either 'linear' or 'quadratic'
+ * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle'
+ * @returns {object} Map from side index to the local node indices on that side
+ */
+export function getBoundarySides(meshDimension, elementOrder, elementShape = "quad") {
+  if (meshDimension === "1D") return lineBoundarySides[elementOrder];
+  if (elementShape === "triangle") return triangleBoundarySides[elementOrder];
+  return quadBoundarySides[elementOrder];
+}
+
 /**
  * Class to handle basis functions and their derivatives based on element configuration
  */

src/elements/elementGeometry.js

@@ -0,0 +1,129 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+/**
+ * Function to test if a point is inside a triangle using barycentric coordinates,
+ * also returning the natural coordinates (ksi, eta)
+ * @param {number} x - X-coordinate of the point
+ * @param {number} y - Y-coordinate of the point
+ * @param {array} vertices - Triangle vertices [[x0,y0],[x1,y1],[x2,y2]]
+ * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta}
+ */
+export function pointInsideTriangle(x, y, vertices) {
+  const tolerance = 1e-12;
+  const [v0, v1, v2] = vertices;
+
+  const denom = (v1[1] - v2[1]) * (v0[0] - v2[0]) + (v2[0] - v1[0]) * (v0[1] - v2[1]);
+
+  const ksi = ((v1[1] - v2[1]) * (x - v2[0]) + (v2[0] - v1[0]) * (y - v2[1])) / denom;
+  const eta = ((v2[1] - v0[1]) * (x - v2[0]) + (v0[0] - v2[0]) * (y - v2[1])) / denom;
+  const gamma = 1 - ksi - eta;
+
+  const inside = ksi >= -tolerance && eta >= -tolerance && gamma >= -tolerance;
+  return { inside, ksi, eta };
+}
+
+/**
+ * Function to test if a point is inside a triangular element, also returning the natural coordinates
+ * (ksi, eta) of the FEAScript reference triangle
+ * @param {number} x - X-coordinate of the point
+ * @param {number} y - Y-coordinate of the point
+ * @param {array} vertices - Triangle vertices in FEAScript order [[x0,y0],[x1,y1],[x2,y2]]
+ * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta}
+ */
+export function pointInsideTriangularElement(x, y, vertices) {
+  // Vertices order (reference coordinates):
+  //   1 (0,1)
+  //   | \
+  //   |   \
+  //   0 --- 2 (1,0)
+  // The barycentric weights of vertices 2 and 1 are the natural coordinates ksi and eta
+  const { inside, ksi: weight0, eta: weight1 } = pointInsideTriangle(x, y, vertices);
+  return { inside, ksi: 1 - weight0 - weight1, eta: weight1 };
+}
+
+/**
+ * Function to test if a point is inside a quadrilateral by spliting it into triangles and using barycentric coordinates
+ * @param {number} x - X-coordinate of the point
+ * @param {number} y - Y-coordinate of the point
+ * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]]
+ * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta}
+ */
+export function pointInsideQuadrilateral(x, y, vertices) {
+  const [firstTriangleVertices, secondTriangleVertices] = splitQuadrilateral(vertices);
+  const pointInsideFirstTriangle = pointInsideTriangle(x, y, firstTriangleVertices);
+  const pointInsideSecondTriangle = pointInsideTriangle(x, y, secondTriangleVertices);
+
+  const inside = pointInsideFirstTriangle.inside || pointInsideSecondTriangle.inside;
+  let ksi = 0;
+  let eta = 0;
+
+  if (inside) {
+    const [v0, v1, v2, v3] = vertices;
+
+    // Function to calculate distance from point to line segment
+    const getDistanceFromLine = (p1, p2) => {
+      const num = Math.abs((p2[0] - p1[0]) * (p1[1] - y) - (p1[0] - x) * (p2[1] - p1[1]));
+      const den = Math.sqrt((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2);
+      return num / den;
+    };
+
+    // Calculate distances to edges based on vertex order:
+    //   1 (v1) --- 3 (v3)
+    //   |          |
+    //   0 (v0) --- 2 (v2)
+
+    const distLeft = getDistanceFromLine(v0, v1);
+    const distRight = getDistanceFromLine(v2, v3);
+    const distBottom = getDistanceFromLine(v0, v2);
+    const distTop = getDistanceFromLine(v1, v3);
+
+    ksi = distLeft / (distLeft + distRight);
+    eta = distBottom / (distBottom + distTop);
+
+    // The distance ratios are exact only for parallelograms, so refine them with Newton iterations
+    // on the bilinear mapping x(ksi, eta) = v0 (1-ksi)(1-eta) + v1 (1-ksi) eta + v2 ksi (1-eta) + v3 ksi eta
+    for (let iteration = 0; iteration < 20; iteration++) {
+      const mappedX =
+        v0[0] * (1 - ksi) * (1 - eta) + v1[0] * (1 - ksi) * eta + v2[0] * ksi * (1 - eta) + v3[0] * ksi * eta;
+      const mappedY =
+        v0[1] * (1 - ksi) * (1 - eta) + v1[1] * (1 - ksi) * eta + v2[1] * ksi * (1 - eta) + v3[1] * ksi * eta;
+      const ksiDerivX = (v2[0] - v0[0]) * (1 - eta) + (v3[0] - v1[0]) * eta;
+      const ksiDerivY = (v2[1] - v0[1]) * (1 - eta) + (v3[1] - v1[1]) * eta;
+      const etaDerivX = (v1[0] - v0[0]) * (1 - ksi) + (v3[0] - v2[0]) * ksi;
+      const etaDerivY = (v1[1] - v0[1]) * (1 - ksi) + (v3[1] - v2[1]) * ksi;
+      const detJacobian = ksiDerivX * etaDerivY - etaDerivX * ksiDerivY;
+      const residualX = x - mappedX;
+      const residualY = y - mappedY;
+      const deltaKsi = (etaDerivY * residualX - etaDerivX * residualY) / detJacobian;
+      const deltaEta = (ksiDerivX * residualY - ksiDerivY * residualX) / detJacobian;
+      ksi += deltaKsi;
+      eta += deltaEta;
+      if (Math.abs(deltaKsi) + Math.abs(deltaEta) < 1e-14) break;
+    }
+  }
+
+  return { inside, ksi, eta };
+}
+
+/**
+ * Function to split the quadrilateral elements into two triangles
+ * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]]
+ * @returns {array} Array of triangle vertices: [[v0,v1,v3], [v0,v2,v3]]
+ */
+export function splitQuadrilateral(vertices) {
+  const [v0, v1, v2, v3] = vertices;
+  // Vertices order:
+  //   1 --- 3
+  //   |     |
+  //   0 --- 2
+  return [
+    [v0, v1, v3],
+    [v0, v2, v3],
+  ];
+}

src/elements/feaInitialization.js

@@ -0,0 +1,71 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+// Internal imports
+import { BasisFunctions } from "./basisFunctions.js";
+import { NumericalIntegration } from "./numericalIntegration.js";
+import { getIntegrationPoints2D } from "./isoparametricMapping.js";
+import { SparseMatrix } from "../linalg/sparseMatrix.js";
+/**
+ * Function to initialize the FEA matrices and numerical tools
+ * @param {object} meshData - Object containing mesh data from prepareMesh()
+ * @returns {object} An object containing initialized matrices and numerical tools
+ */
+export function initializeFEA(meshData) {
+  const { totalNodes, nop, meshDimension, elementOrder, elementShape } = meshData;
+
+  // Initialize variables for matrix assembly
+  let residualVector = new Array(totalNodes).fill(0);
+  let jacobianMatrix = new SparseMatrix(totalNodes);
+  let localToGlobalMap = [];
+
+  // Initialize the BasisFunctions class
+  const basisFunctions = new BasisFunctions({
+    meshDimension,
+    elementOrder,
+    elementShape,
+  });
+
+  // Initialize the NumericalIntegration class
+  const numericalIntegration = new NumericalIntegration({
+    meshDimension,
+    elementOrder,
+    elementShape,
+  });
+
+  // Calculate Gauss points and weights
+  let gaussPointsAndWeights = numericalIntegration.getGaussPointsAndWeights();
+  let gaussPoints = gaussPointsAndWeights.gaussPoints;
+  let gaussWeights = gaussPointsAndWeights.gaussWeights;
+
+  // Integration points for 2D elements (tensor product for quads, collapsed rule for triangles)
+  const integrationPoints2D =
+    meshDimension === "2D" ? getIntegrationPoints2D(gaussPoints, gaussWeights, elementShape) : [];
+
+  // Basis functions on the reference element at each integration point (the same for every element)
+  const referenceBasisFunctions =
+    meshDimension === "1D"
+      ? gaussPoints.map((gaussPoint) => basisFunctions.getBasisFunctions(gaussPoint))
+      : integrationPoints2D.map(([ksi, eta]) => basisFunctions.getBasisFunctions(ksi, eta));
+
+  // Determine the number of nodes in the reference element based on the first element in the nop array
+  const nodesPerElement = nop[0].length;
+
+  return {
+    residualVector,
+    jacobianMatrix,
+    localToGlobalMap,
+    basisFunctions,
+    gaussPoints,
+    gaussWeights,
+    integrationPoints2D,
+    referenceBasisFunctions,
+    nodesPerElement,
+  };
+}