Re-organize folder stucture

CONTRIBUTORS.md

@@ -0,0 +1,10 @@
+## Contributors
+
+(alphabetical order)
+
+- [Akhilesh Tiwari](https://github.com/silky-x0)
+- [Edoardo R](https://github.com/Ed6003)
+- [Felipe Ferrari](https://github.com/ferrari212)
+- [Nikolaos Chamakos](https://www.npmjs.com/~nikoscham)
+- [Priyank Bansal](https://github.com/Stu-ops)
+- [sridhar-mani](https://github.com/sridhar-mani)

examples/creepingFlowScript/lidDrivenCavity2DCreepingFlow/lidDrivenCavity2DCreepingFlow.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/eulerBernoulliBeamScript/clampedSpringSupportedBeam1D/clampedSpringSupportedBeam1D.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 
@@ -50,9 +46,7 @@ model.setSolverMethod("lusolve");
 const { solutionVector } = model.solve();
 
 // Print results
-const flatSolution = solutionVector.map((entry) =>
-  Array.isArray(entry) ? entry[0] : entry
-);
+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)");
@@ -61,10 +55,8 @@ 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
+    `  ${nodeIndex + 1}  | ${nodeXCoordinates[nodeIndex].toFixed(2).padStart(8)} | ${w
       .toExponential(4)
-      .padStart(20)}`
+      .padStart(17)} | ${theta.toExponential(4).padStart(20)}`,
   );
 }

examples/frontPropagationScript/solidificationFront2D/solidificationFront2D.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/generalFormPDEScript/advectionDiffusion1D/advectionDiffusion1D.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/generalFormPDEScript/nonlinearReactionDiffusion1D/nonlinearReactionDiffusion1D.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction1DWall/heatConduction1DVaryingK.js

@@ -19,10 +19,6 @@
  *   - Right (x = 0.15): constant temperature, T = 20 °C
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction1DWall/heatConduction1DWall.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction2DFin/heatConduction2DFin.js

@@ -7,10 +7,6 @@
  * ════════════════════════════════════════════════════════════════
  */
 
-// Import Math.js
-import * as math from "mathjs";
-globalThis.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction2DFin/heatConduction2DFinGmsh.js

@@ -12,10 +12,6 @@ import fs from "fs";
 import path from "path";
 import { fileURLToPath } from "url";
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, importGmshMesh, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction2DFin/heatConduction2DRhomFinGmsh.js

@@ -12,10 +12,6 @@ import fs from "fs";
 import path from "path";
 import { fileURLToPath } from "url";
 
-// Import Math.js
-import * as math from "mathjs";
-global.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, importGmshMesh, printVersion } from "feascript";
 

examples/heatConductionScript/heatConduction2DFin/heatConduction2DVaryingKQ.js

@@ -27,10 +27,6 @@
  *   - Right  (x = 4):  constant temperature, T = 200 °C
  */
 
-// Import Math.js
-import * as math from "mathjs";
-globalThis.math = math;
-
 // Import FEAScript library
 import { FEAScriptModel, printVersion } from "feascript";
 

package.json

@@ -39,24 +39,24 @@
   "author": "Nikolaos Chamakos <nikoscham@gmail.com> (https://www.npmjs.com/~nikoscham)",
   "contributors": [
     {
-      "name": "sridhar-mani",
-      "url": "https://www.npmjs.com/~sridhar-mani"
+      "name": "Akhilesh Tiwari",
+      "url": "https://github.com/silky-x0"
     },
     {
-      "name": "Felipe Ferrari",
-      "url": "https://github.com/ferrari212"
+      "name": "Edoardo R",
+      "url": "https://github.com/Ed6003"
     },
     {
-      "name": "Akhilesh Tiwari",
-      "url": "https://github.com/silky-x0"
+      "name": "Felipe Ferrari",
+      "url": "https://github.com/ferrari212"
     },
     {
       "name": "Priyank Bansal",
       "url": "https://github.com/Stu-ops"
     },
     {
-      "name": "Edoardo R",
-      "url": "https://github.com/Ed6003"
+      "name": "sridhar-mani",
+      "url": "https://github.com/sridhar-mani"
     }
   ],
   "license": "MIT",

src/base/linalg/index.js

@@ -0,0 +1,10 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+export { euclideanNorm, dotProduct, copyVector, axpy } from "./blasUtilities.js";

src/base/solvers/index.js

@@ -0,0 +1,11 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+export { jacobiSolver } from "./jacobiSolver.js";
+export { mathjsLUSolver } from "./mathjsLUSolver.js";

src/methods/jacobiSolver.js → src/base/solvers/jacobiSolver.js

@@ -8,48 +8,48 @@
  */
 
 // Internal imports
-import { dotProduct, copyVector, euclideanNorm } from "./blasUtilities.js";
+import { dotProduct, copyVector, euclideanNorm } from "../linalg/blasUtilities.js";
 
 /**
  * Function to solve a system of linear equations using the Jacobi iterative method (CPU synchronous version)
- * @param {array} A - The system matrix
- * @param {array} b - The right-hand side vector
- * @param {array} x0 - Initial guess for solution vector
+ * @param {array} systemMatrix - The system matrix
+ * @param {array} rightHandSideVector - The right-hand side vector
+ * @param {array} initialGuess - Initial guess for solution vector
  * @param {object} [options] - Optional parameters for the solver, such as `maxIterations` and `tolerance`
  * @returns {object} An object containing:
  *  - solutionVector: The solution vector
  *  - iterations: The number of iterations performed
  *  - converged: Boolean indicating whether the method converged
  */
-export function jacobiSolver(A, b, x0, options = {}) {
+export function jacobiSolver(systemMatrix, rightHandSideVector, initialGuess, options = {}) {
   // Extract options
   const { maxIterations, tolerance } = options;
 
-  const n = A.length;
+  const n = systemMatrix.length;
 
   // Convert inputs to Float64Arrays for BLAS operations
-  const Arows = A.map((row) => new Float64Array(row));
-  const bVec = new Float64Array(b);
-  let x = new Float64Array(x0);
-  let xNew = new Float64Array(n);
+  const rows = systemMatrix.map((row) => new Float64Array(row));
+  const rhs = new Float64Array(rightHandSideVector);
+  let solutionVector = new Float64Array(initialGuess);
+  let updatedSolutionVector = new Float64Array(n);
   const diff = new Float64Array(n);
 
   // Jacobi update: xNew[i] = (b[i] - (A[i] · x) + A[i][i] * x[i]) / A[i][i]
   for (let iter = 0; iter < maxIterations; iter++) {
     for (let i = 0; i < n; i++) {
-      const rowDot = dotProduct(Arows[i], x);
-      xNew[i] = (bVec[i] - rowDot + Arows[i][i] * x[i]) / Arows[i][i];
+      const rowDot = dotProduct(rows[i], solutionVector);
+      updatedSolutionVector[i] = (rhs[i] - rowDot + rows[i][i] * solutionVector[i]) / rows[i][i];
     }
 
-    // Compute diff and copy xNew into x
-    for (let i = 0; i < n; i++) diff[i] = xNew[i] - x[i];
+    // Compute diff and copy updatedSolutionVector into solutionVector
+    for (let i = 0; i < n; i++) diff[i] = updatedSolutionVector[i] - solutionVector[i];
     const residual = euclideanNorm(diff);
-    copyVector(xNew, x);
+    copyVector(updatedSolutionVector, solutionVector);
 
     if (residual < tolerance) {
-      return { solutionVector: x, iterations: iter + 1, converged: true };
+      return { solutionVector, iterations: iter + 1, converged: true };
     }
   }
 
-  return { solutionVector: x, iterations: maxIterations, converged: false };
+  return { solutionVector, iterations: maxIterations, converged: false };
 }

src/base/solvers/mathjsLUSolver.js

@@ -0,0 +1,39 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+// External imports
+import { create, all } from "mathjs";
+const math = create(all);
+
+/**
+ * Function to solve a system of linear equations using LU decomposition (mathjs)
+ * @param {array} systemMatrix - The system matrix
+ * @param {array} rightHandSideVector - The right-hand side vector
+ * @param {object} [options] - Optional parameters for the solver, such as `matrixType` ("sparse" or "dense")
+ * @returns {object} An object containing:
+ *  - solutionVector: The solution vector
+ */
+export function mathjsLUSolver(systemMatrix, rightHandSideVector, options = {}) {
+  // Extract options
+  const { matrixType } = options;
+
+  let solutionVector;
+  if (matrixType === "sparse") {
+    const sparseMatrix = math.sparse(systemMatrix);
+    const luFactorization = math.slu(sparseMatrix, 1, 1); // order=1, threshold=1 for pivoting
+    const solutionMatrix = math.lusolve(luFactorization, rightHandSideVector);
+    solutionVector = math.squeeze(solutionMatrix).valueOf();
+  } else {
+    // Dense matrix
+    const solutionMatrix = math.lusolve(systemMatrix, rightHandSideVector);
+    solutionVector = math.squeeze(solutionMatrix).valueOf();
+  }
+
+  return { solutionVector };
+}

src/methods/linearSystemSolver.js

@@ -8,7 +8,7 @@
  */
 
 // Internal imports
-import { jacobiSolver } from "./jacobiSolver.js";
+import { jacobiSolver, mathjsLUSolver } from "../base/solvers/index.js";
 import { basicLog, debugLog, errorLog } from "../utilities/logging.js";
 import * as Comlink from "../vendor/comlink.mjs";
 
@@ -37,11 +37,8 @@ export function solveLinearSystem(solverMethod, jacobianMatrix, residualVector,
 
   if (solverMethod === "lusolve") {
     // Use LU decomposition method
-    const jacobianMatrixSparse = math.sparse(jacobianMatrix);
-    const luFactorization = math.slu(jacobianMatrixSparse, 1, 1); // order=1, threshold=1 for pivoting
-    let solutionMatrix = math.lusolve(luFactorization, residualVector);
-    solutionVector = math.squeeze(solutionMatrix).valueOf();
-    //solutionVector = math.lusolve(jacobianMatrix, residualVector); // In the case of a dense matrix
+    const luSolverResult = mathjsLUSolver(jacobianMatrix, residualVector, { matrixType: "sparse" });
+    solutionVector = luSolverResult.solutionVector;
   } else if (solverMethod === "jacobi") {
     // Use Jacobi method
     const initialGuess = new Array(residualVector.length).fill(0);
@@ -102,8 +99,12 @@ export async function solveLinearSystemAsync(solverMethod, jacobianMatrix, resid
   console.time("systemSolving");
 
   // Normalize inputs
-  const A = Array.isArray(jacobianMatrix) ? jacobianMatrix : (jacobianMatrix?.toArray?.() ?? jacobianMatrix);
-  const b = Array.isArray(residualVector) ? residualVector : (residualVector?.toArray?.() ?? residualVector);
+  const systemMatrix = Array.isArray(jacobianMatrix)
+    ? jacobianMatrix
+    : jacobianMatrix?.toArray?.() ?? jacobianMatrix;
+  const rightHandSideVector = Array.isArray(residualVector)
+    ? residualVector
+    : residualVector?.toArray?.() ?? residualVector;
 
   let created = null;
   let computeEngine = null;
@@ -117,10 +118,13 @@ export async function solveLinearSystemAsync(solverMethod, jacobianMatrix, resid
     created = await createDefaultComputeEngine();
     computeEngine = created.computeEngine;
 
-    const x0 = new Array(b.length).fill(0);
+    const initialGuess = new Array(rightHandSideVector.length).fill(0);
     let result;
 
-    result = await computeEngine.webgpuJacobiSolver(A, b, x0, { maxIterations, tolerance });
+    result = await computeEngine.webgpuJacobiSolver(systemMatrix, rightHandSideVector, initialGuess, {
+      maxIterations,
+      tolerance,
+    });
     solutionVector = result.solutionVector;
     converged = result.converged;
     iterations = result.iterations;

src/methods/newtonRaphson.js

@@ -8,7 +8,7 @@
  */
 
 // Internal imports
-import { euclideanNorm, axpy } from "../methods/blasUtilities.js";
+import { euclideanNorm, axpy } from "../base/linalg/index.js";
 import { solveLinearSystem } from "./linearSystemSolver.js";
 import { basicLog, debugLog, errorLog } from "../utilities/logging.js";
 import { runFrontalSolver } from "./frontalSolver.js";