Add support 1D nonlinear problems.

README.md

@@ -17,7 +17,7 @@
 
 ## Features
 
-- <b>Physics models:</b> creeping (Stokes) flow, front propagation, heat conduction
+- <b>Physics models:</b> creeping (Stokes) flow, Euler-Bernoulli beam bending, front propagation, heat conduction, general form PDE (linear and nonlinear)
 - <b>Meshing:</b> simple 1D/2D mesh generation, unstructured mesh import from Gmsh (`.msh`)
 - <b>Solvers:</b> frontal, Jacobi (CPU/WebGPU) and LU, Newton–Raphson for nonlinear systems
 - <b>Performance:</b> web worker support for multi-threaded computation

examples/generalFormPDEScript/README.md

@@ -10,6 +10,10 @@ This directory contains Node.js examples demonstrating how to use the FEAScript
 
 This example demonstrates solving a one-dimensional advection-diffusion problem with a Gaussian source term. The problem models the transport of a substance under the effects of both diffusion and advection. For detailed information on the model setup, refer to the corresponding [tutorial](https://feascript.com/tutorials/advection-diffusion-1d.html) in the FEAScript website.
 
+#### 2. Nonlinear Reaction-Diffusion (`nonlinearReactionDiffusion1D.js`)
+
+This example demonstrates solving a one-dimensional nonlinear reaction-diffusion problem with a quadratic reaction term, using the Newton-Raphson method (`nonlinear: true`).
+
 ## Running the Node.js Examples
 
 #### 1. Create package.json with ES module support:

examples/generalFormPDEScript/nonlinearReactionDiffusion1D/nonlinearReactionDiffusion1D.js

@@ -0,0 +1,62 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+// 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();
+
+// Reaction rate coefficient for the nonlinear source term
+const Da = 1;
+
+// Select physics/PDE
+model.setModelConfig("generalFormPDEScript", {
+  nonlinear: true, // Solve with the Newton-Raphson method
+  coefficientFunctions: {
+    // Equation d²u/dx² - Da * u² = 0
+    A: (x) => 1, // Diffusion coefficient
+    B: (x) => 0, // Advection coefficient
+    C: (x) => 0, // Linear reaction coefficient
+    D: (x, u) => Da * u ** 2, // Nonlinear reaction/source term
+    dDdu: (x, u) => 2 * Da * u, // Derivative of D with respect to u, required for the Jacobian
+  },
+});
+
+// Define mesh configuration
+model.setMeshConfig({
+  meshDimension: "1D",
+  elementOrder: "linear",
+  numElementsX: 20,
+  maxX: 10.0,
+});
+
+// Define boundary conditions
+model.addBoundaryCondition("0", ["constantValue", 1]); // Left boundary
+model.addBoundaryCondition("1", "zeroGradient"); // Right boundary
+
+// Set solver method
+model.setSolverMethod("lusolve");
+
+// Solve the problem
+const { solutionVector, nodesCoordinates } = model.solve({
+  maxIterations: 100,
+  tolerance: 1e-5,
+});
+
+// Print results
+console.log(`Number of nodes in mesh: ${nodesCoordinates.nodesXCoordinates.length}`);
+console.log("Node coordinates:", nodesCoordinates);
+console.log("Solution vector:", solutionVector);

src/FEAScript.js

@@ -13,7 +13,11 @@ import { solveLinearSystem } from "./methods/linearSystemSolver.js";
 import { solveLinearSystemAsync } from "./methods/linearSystemSolver.js";
 import { prepareMesh } from "./mesh/meshUtils.js";
 import { assembleFrontPropagationMat } from "./models/frontPropagation.js";
-import { assembleGeneralFormPDEMat, assembleGeneralFormPDEFront } from "./models/generalFormPDE.js";
+import {
+  assembleGeneralFormPDEMat,
+  assembleGeneralFormPDEFront,
+  assembleGeneralFormPDENonlinearMat,
+} from "./models/generalFormPDE.js";
 import { assembleHeatConductionMat, assembleHeatConductionFront } from "./models/heatConduction.js";
 import { assembleCreepingFlowMatrix } from "./models/creepingFlow.js";
 import { assembleEulerBernoulliBeamMat } from "./models/eulerBernoulliBeam.js";
@@ -53,6 +57,11 @@ export class FEAScriptModel {
       this.coefficientFunctions = options.coefficientFunctions;
       debugLog("coefficientFunctions set");
     }
+    // Flag to solve the PDE with the Newton-Raphson method instead of a direct linear solve
+    if (options?.nonlinear !== undefined) {
+      this.nonlinear = options.nonlinear;
+      debugLog(`nonlinear set to ${this.nonlinear}`);
+    }
     // Only update if a value is provided
     if (options?.maxIterations !== undefined) {
       this.maxIterations = options.maxIterations;
@@ -181,6 +190,30 @@ export class FEAScriptModel {
         errorLog(
           "Frontal solver is not yet supported for generalFormPDEScript. Please use 'lusolve' or 'jacobi'.",
         );
+      } else if (this.nonlinear) {
+        // Solve the nonlinear PDE with the Newton-Raphson method
+        const context = {
+          meshData,
+          boundaryConditions: this.boundaryConditions,
+          solverMethod: this.solverMethod,
+          maxIterations: options.maxIterations ?? this.maxIterations,
+          tolerance: options.tolerance ?? this.tolerance,
+        };
+
+        const newtonRaphsonResult = newtonRaphson(
+          (meshDataArg, boundaryConditionsArg, solutionVectorArg) =>
+            assembleGeneralFormPDENonlinearMat(
+              meshDataArg,
+              boundaryConditionsArg,
+              this.coefficientFunctions,
+              solutionVectorArg,
+            ),
+          context,
+        );
+
+        jacobianMatrix = newtonRaphsonResult.jacobianMatrix;
+        residualVector = newtonRaphsonResult.residualVector;
+        solutionVector = newtonRaphsonResult.solutionVector;
       } else {
         // Use regular linear solver methods
         ({ jacobianMatrix, residualVector } = assembleGeneralFormPDEMat(

src/models/generalFormPDE.js

@@ -259,3 +259,146 @@ export function assembleGeneralFormPDEFront({
     ngl,
   };
 }
+
+/**
+ * Function to assemble the Jacobian matrix and residual vector for the nonlinear general form PDE model,
+ * to be used with the Newton-Raphson solver
+ * @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) for the linear terms, plus D(x, u) and
+ *  dDdu(x, u) for the nonlinear reaction/source term and its derivative with respect to u
+ * @param {array} solutionVector - The current solution vector (Newton-Raphson iterate)
+ * @returns {object} An object containing:
+ *  - jacobianMatrix: The assembled Jacobian matrix (negative of dResidual/du)
+ *  - residualVector: The assembled residual vector
+ */
+export function assembleGeneralFormPDENonlinearMat(
+  meshData,
+  boundaryConditions,
+  coefficientFunctions,
+  solutionVector,
+) {
+  basicLog("Starting nonlinear general form PDE matrix assembly...");
+
+  // Extract mesh data
+  const {
+    nodesXCoordinates,
+    nodesYCoordinates,
+    nop,
+    boundaryElements,
+    totalElements,
+    meshDimension,
+    elementOrder,
+  } = meshData;
+
+  // Extract coefficient functions
+  const { A, B, C, D, dDdu } = coefficientFunctions;
+
+  // Initialize FEA components
+  const FEAData = initializeFEA(meshData);
+  const {
+    residualVector,
+    jacobianMatrix,
+    localToGlobalMap,
+    basisFunctions,
+    gaussPoints,
+    gaussWeights,
+    nodesPerElement,
+  } = FEAData;
+
+  if (meshDimension === "1D") {
+    // 1D nonlinear general form PDE
+
+    // Matrix assembly
+    for (let elementIndex = 0; elementIndex < totalElements; elementIndex++) {
+      // Map local element nodes to global mesh nodes
+      for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) {
+        // Convert to 0-based indexing
+        localToGlobalMap[localNodeIndex] = Math.abs(nop[elementIndex][localNodeIndex]) - 1;
+      }
+
+      // Loop over Gauss points
+      for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) {
+        // Get basis functions for the current Gauss point
+        const { basisFunction, basisFunctionDerivKsi } = basisFunctions.getBasisFunctions(
+          gaussPoints[gaussPointIndex],
+        );
+
+        // Perform isoparametric mapping
+        const { detJacobian, basisFunctionDerivX } = performIsoparametricMapping1D({
+          basisFunction,
+          basisFunctionDerivKsi,
+          nodesXCoordinates,
+          localToGlobalMap,
+          nodesPerElement,
+        });
+
+        // Calculate the physical coordinate, solution value and solution derivative at this Gauss point
+        let xCoord = 0;
+        let uValue = 0;
+        let uDerivX = 0;
+        for (let i = 0; i < nodesPerElement; i++) {
+          xCoord += nodesXCoordinates[localToGlobalMap[i]] * basisFunction[i];
+          uValue += solutionVector[localToGlobalMap[i]] * basisFunction[i];
+          uDerivX += solutionVector[localToGlobalMap[i]] * basisFunctionDerivX[i];
+        }
+
+        // Evaluate coefficient functions at this physical coordinate and solution state
+        const a = A(xCoord);
+        const b = B(xCoord);
+        const c = C(xCoord);
+        const d = D(xCoord, uValue);
+        const dDduVal = dDdu(xCoord, uValue);
+
+        // Computation of the residual vector and the Newton-Raphson Jacobian matrix
+        for (let localNodeIndex1 = 0; localNodeIndex1 < nodesPerElement; localNodeIndex1++) {
+          const globalNodeIndex1 = localToGlobalMap[localNodeIndex1];
+
+          // Residual contribution (diffusion, advection, reaction and nonlinear source terms)
+          residualVector[globalNodeIndex1] +=
+            gaussWeights[gaussPointIndex] *
+            detJacobian *
+            (a * uDerivX * basisFunctionDerivX[localNodeIndex1] -
+              b * uDerivX * basisFunction[localNodeIndex1] -
+              c * uValue * basisFunction[localNodeIndex1] +
+              d * basisFunction[localNodeIndex1]);
+
+          for (let localNodeIndex2 = 0; localNodeIndex2 < nodesPerElement; localNodeIndex2++) {
+            const globalNodeIndex2 = localToGlobalMap[localNodeIndex2];
+
+            // Jacobian is the negative of the residual derivative, matching the Newton-Raphson solver convention
+            jacobianMatrix[globalNodeIndex1][globalNodeIndex2] +=
+              -gaussWeights[gaussPointIndex] *
+              detJacobian *
+              (a * basisFunctionDerivX[localNodeIndex1] * basisFunctionDerivX[localNodeIndex2] -
+                b * basisFunctionDerivX[localNodeIndex2] * basisFunction[localNodeIndex1] -
+                c * basisFunction[localNodeIndex1] * basisFunction[localNodeIndex2] +
+                dDduVal * basisFunction[localNodeIndex1] * basisFunction[localNodeIndex2]);
+          }
+        }
+      }
+    }
+  } else if (meshDimension === "2D") {
+    errorLog("2D nonlinear general form PDE is not yet supported in assembleGeneralFormPDENonlinearMat.");
+    // 2D nonlinear general form PDE - empty for now
+  }
+
+  // Apply boundary conditions
+  const genericBoundaryConditions = new GenericBoundaryConditions(
+    boundaryConditions,
+    boundaryElements,
+    nop,
+    meshDimension,
+    elementOrder,
+  );
+
+  // Apply Dirichlet boundary conditions only (as a Newton-Raphson increment, since solutionVector is passed)
+  genericBoundaryConditions.imposeDirichletBoundaryConditions(residualVector, jacobianMatrix, solutionVector);
+
+  basicLog("Nonlinear general form PDE matrix assembly completed");
+
+  return {
+    jacobianMatrix,
+    residualVector,
+  };
+}

src/models/genericBoundaryConditions.js

@@ -34,6 +34,8 @@ export class GenericBoundaryConditions {
    * Function to impose Dirichlet boundary conditions
    * @param {array} residualVector - The residual vector to be modified
    * @param {array} jacobianMatrix - The Jacobian matrix to be modified
+   * @param {array} [solutionVector] - Current solution (Newton-Raphson iterate); when provided, the
+   *  residual is set to the increment needed to reach the prescribed value instead of the value itself
    *
    * For consistency across both linear and nonlinear formulations,
    * this project always refers to the assembled right-hand side vector
@@ -43,7 +45,7 @@ export class GenericBoundaryConditions {
    * classic stiffness/conductivity matrix and `residualVector`
    * corresponds to the traditional load (RHS) vector.
    */
-  imposeDirichletBoundaryConditions(residualVector, jacobianMatrix) {
+  imposeDirichletBoundaryConditions(residualVector, jacobianMatrix, solutionVector) {
     if (this.meshDimension === "1D") {
       Object.keys(this.boundaryConditions).forEach((boundaryKey) => {
         if (this.boundaryConditions[boundaryKey][0] === "constantValue") {
@@ -62,8 +64,10 @@ export class GenericBoundaryConditions {
                     elementIndex + 1
                   }, local node ${nodeIndex + 1})`,
                 );
-                // Set the residual vector to the value
-                residualVector[globalNodeIndex] = value;
+                // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson
+                residualVector[globalNodeIndex] = solutionVector
+                  ? value - solutionVector[globalNodeIndex]
+                  : value;
                 // Set the Jacobian matrix row to zero
                 for (let colIndex = 0; colIndex < residualVector.length; colIndex++) {
                   jacobianMatrix[globalNodeIndex][colIndex] = 0;
@@ -83,8 +87,10 @@ export class GenericBoundaryConditions {
                     elementIndex + 1
                   }, local node ${nodeIndex + 1})`,
                 );
-                // Set the residual vector to the value
-                residualVector[globalNodeIndex] = value;
+                // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson
+                residualVector[globalNodeIndex] = solutionVector
+                  ? value - solutionVector[globalNodeIndex]
+                  : value;
                 // Set the Jacobian matrix row to zero
                 for (let colIndex = 0; colIndex < residualVector.length; colIndex++) {
                   jacobianMatrix[globalNodeIndex][colIndex] = 0;
@@ -116,8 +122,10 @@ export class GenericBoundaryConditions {
                     elementIndex + 1
                   }, local node ${nodeIndex + 1})`,
                 );
-                // Set the residual vector to the value
-                residualVector[globalNodeIndex] = value;
+                // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson
+                residualVector[globalNodeIndex] = solutionVector
+                  ? value - solutionVector[globalNodeIndex]
+                  : value;
                 // Set the Jacobian matrix row to zero
                 for (let colIndex = 0; colIndex < residualVector.length; colIndex++) {
                   jacobianMatrix[globalNodeIndex][colIndex] = 0;
@@ -139,8 +147,10 @@ export class GenericBoundaryConditions {
                     elementIndex + 1
                   }, local node ${nodeIndex + 1})`,
                 );
-                // Set the residual vector to the value
-                residualVector[globalNodeIndex] = value;
+                // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson
+                residualVector[globalNodeIndex] = solutionVector
+                  ? value - solutionVector[globalNodeIndex]
+                  : value;
                 // Set the Jacobian matrix row to zero
                 for (let colIndex = 0; colIndex < residualVector.length; colIndex++) {
                   jacobianMatrix[globalNodeIndex][colIndex] = 0;

src/visualization/plotSolution.js

@@ -8,8 +8,7 @@
  */
 
 /**
- * Plotly-based plotting functions are provided by plotlyPlot.js
- * VTK.js-based plotting and data transformation functions are provided by vtkPlot.js
+ * Function to re-export plotting functions in order to visualize solution fields
  */
 
 export { plotSolution, plotInterpolatedSolution } from "./plotlyPlot.js";