FEAScript.js

/**
 * ════════════════════════════════════════════════════════════════
 *  FEAScript Core Library
 *  Finite Element Simulation in JavaScript
 *  Version: 0.3.0 (RC) | https://feascript.com
 *  MIT License © 2023–2026 FEAScript
 * ════════════════════════════════════════════════════════════════
 */

// Internal imports
import { newtonRaphson } from "./solvers/newtonRaphson.js";
import { solveLinearSystem } from "./solvers/linearSystemSolver.js";
import { solveLinearSystemAsync } from "./solvers/linearSystemSolver.js";
import { prepareMesh, resolveMeshConfig } from "./mesh/meshUtils.js";
import { assembleFrontPropagationMat, assembleFrontPropagationFront } from "./physics/frontPropagation.js";
import {
  assembleGeneralFormPDEMat,
  assembleGeneralFormPDEFront,
  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, extractCreepingFlowFields } from "./physics/creepingFlow.js";
import {
  assembleEulerBernoulliBeamMat,
  extractEulerBernoulliBeamFields,
} from "./physics/eulerBernoulliBeam.js";
import { runFrontalSolver } from "./solvers/frontalSolver.js";
import { basicLog, debugLog, warnLog, errorLog } from "./utilities/logging.js";

// Boundary condition types of each model (the Euler-Bernoulli beam model validates its own conditions)
const supportedBoundaryConditionTypes = {
  heatConductionScript: ["constantTemperature", "convection", "symmetry"],
  generalFormPDEScript: ["constantValue", "zeroGradient"],
  frontPropagationScript: ["constantValue", "zeroGradient"],
  creepingFlowScript: ["constantVelocity", "stressFree"],
};

/**
 * Class to implement finite element analysis in JavaScript
 * @param {string} solverConfig - Parameter specifying the type of solver
 * @param {object} meshConfig - Object containing computational mesh details
 * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis
 * @returns {object} An object containifng the solution vector and mesh information
 */
export class FEAScriptModel {
  constructor() {
    this.solverConfig = null;
    this.meshConfig = {};
    this.boundaryConditions = {};
    this.solverMethod = "lusolve"; // Default solver method
    this.coefficientFunctions = null; // Add storage for coefficient functions
    warnLog(
      "FEAScript is provided “as is” without any warranty. The authors are not responsible for any damages or losses that may result from using the software. See the license for more details: https://github.com/FEAScript/FEAScript-core/blob/main/LICENSE",
    );
    basicLog("FEAScriptModel instance created");
  }

  /**
   * Method to set the model configuration
   * @param {string} modelConfig - Parameter specifying the type of model
   * @param {object} [options] - Optional additional configuration
   */
  setModelConfig(modelConfig, options = {}) {
    this.solverConfig = modelConfig;

    // Store coefficient functions if provided
    if (options?.coefficientFunctions !== undefined) {
      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;
      debugLog(`maxIterations set to ${this.maxIterations}`);
    }
    if (options?.tolerance !== undefined) {
      this.tolerance = options.tolerance;
      debugLog(`tolerance set to ${this.tolerance}`);
    }

    debugLog(`solverConfig set to ${modelConfig}`);
  }

  setMeshConfig(meshConfig) {
    // Settings that a parsed mesh defines (dimension, element order) are filled in from the mesh
    this.meshConfig = resolveMeshConfig(meshConfig);
    debugLog(`meshConfig set with dimensions: ${this.meshConfig.meshDimension}`);
  }

  addBoundaryCondition(boundaryKey, condition) {
    // A condition without parameters may be given as a plain string, e.g. "zeroGradient" instead of ["zeroGradient"]
    if (typeof condition === "string") {
      warnLog(
        `Boundary condition for boundary ${boundaryKey} should be an array: use ["${condition}"] instead of "${condition}".`,
      );
      condition = [condition];
    }
    if (!Array.isArray(condition)) {
      errorLog(
        `Boundary condition for boundary ${boundaryKey} must be an array with the condition type first, e.g. ["constantValue", 0].`,
      );
      return;
    }

    // Normalize deprecated boundary condition type strings and emit deprecation warnings
    const deprecatedBoundaryConditionTypes = {
      constantTemp: "constantTemperature",
    };
    const originalType = condition[0];
    if (Object.prototype.hasOwnProperty.call(deprecatedBoundaryConditionTypes, originalType)) {
      const normalizedType = deprecatedBoundaryConditionTypes[originalType];
      warnLog(
        `Boundary condition type "${originalType}" is deprecated and will be removed in a future version. ` +
          `Use "${normalizedType}" instead.`,
      );
      condition = [normalizedType, ...condition.slice(1)];
    }
    this.boundaryConditions[boundaryKey] = condition;
    debugLog(`boundaryConditions added for boundary: ${boundaryKey}, type: ${condition[0]}`);
  }

  /**
   * Method to report boundary conditions that the selected model does not recognize, which would otherwise
   * be ignored (a boundary without a recognized condition gets the natural condition of the model)
   */
  validateBoundaryConditions() {
    if (this.solverConfig === "eulerBernoulliBeamScript") {
      Object.keys(this.boundaryConditions).forEach((nodeKey) => {
        const conditions = this.boundaryConditions[nodeKey];
        if (!conditions.every(Array.isArray)) {
          errorLog(
            `Boundary conditions of node ${nodeKey} must be an array of conditions, e.g. [["fixed"]] or [["pinned"], ["moment", 1250]].`,
          );
        }
      });
      return;
    }

    const supportedTypes = supportedBoundaryConditionTypes[this.solverConfig];
    if (!supportedTypes) return;
    Object.keys(this.boundaryConditions).forEach((boundaryKey) => {
      const boundaryConditionType = this.boundaryConditions[boundaryKey][0];
      if (!supportedTypes.includes(boundaryConditionType)) {
        errorLog(
          `Unknown boundary condition type "${boundaryConditionType}" on boundary ${boundaryKey} for ${
            this.solverConfig
          }; it is ignored. Supported types: ${supportedTypes.join(", ")}.`,
        );
      }
    });
  }

  setSolverMethod(solverMethod) {
    this.solverMethod = solverMethod;
    debugLog(`solverMethod set to: ${solverMethod}`);
  }

  /**
   * Method to solve the finite element problem synchronously
   * @param {object} [options] - Additional parameters for the solver, such as `maxIterations` and `tolerance`
   * @returns {object} An object containing the solution vector (Float64Array, the same for every solver method) and mesh information.
   *  Models with more than one unknown per node (creeping flow, Euler-Bernoulli beam) also return `fields`, an object of
   *  named nodal fields (one value per mesh node) that can be plotted directly, e.g. `plotSolution(model, result, "contour", id, { field: "u" })`
   */
  solve(options = {}) {
    if (!this.solverConfig || !this.meshConfig || !this.boundaryConditions) {
      errorLog("solverConfig, meshConfig and boundaryConditions must be set before solving");
    }
    this.validateBoundaryConditions();
    /**
     * For consistency across both linear and nonlinear formulations,
     * we always refer 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/conductivity matrix and `residualVector`
     * corresponds to the traditional load (RHS) vector.
     */

    let jacobianMatrix = [];
    let residualVector = [];
    let solutionVector = new Float64Array(0);
    let initialSolution = [];
    let fields = null; // Named nodal fields, for models with more than one unknown per node

    // Prepare the mesh
    basicLog("Preparing mesh...");
    const meshData = prepareMesh(this.meshConfig);
    basicLog("Mesh preparation completed");

    // Extract node coordinates and nodal numbering from meshData
    const nodesCoordinates = {
      nodesXCoordinates: meshData.nodesXCoordinates,
      nodesYCoordinates: meshData.nodesYCoordinates,
    };

    // Select and execute the appropriate solver based on solverConfig
    basicLog("Beginning solving process...");
    console.time("totalSolvingTime");
    basicLog(`Using solver ${this.solverConfig}`);
    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(
          (elementArgs) =>
            assembleHeatConductionFront({
              ...elementArgs,
              coefficientFunctions: this.coefficientFunctions,
            }),
          meshData,
          thermalBoundaryConditions,
        );
        solutionVector = frontalResult.solutionVector;
      } else {
        // Use regular linear solver methods
        ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(
          meshData,
          this.boundaryConditions,
          this.coefficientFunctions,
        ));
        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
          maxIterations: options.maxIterations ?? this.maxIterations,
          tolerance: options.tolerance ?? this.tolerance,
        });
        solutionVector = linearSystemResult.solutionVector;
      }
    } else if (this.solverConfig === "frontPropagationScript") {
      // Initialize eikonalActivationFlag
      let eikonalActivationFlag = 0; // TODO: make activationFlag a generic variable (not only for eikonal)
      const eikonalExteralIterations = 5; // Number of incremental steps for the eikonal equation

      // Create context object with all necessary properties
      const context = {
        meshData: meshData,
        boundaryConditions: this.boundaryConditions,
        eikonalActivationFlag: eikonalActivationFlag,
        solverMethod: this.solverMethod,
        initialSolution,
        // 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
      // exactly (accumulating 1 / eikonalExteralIterations can overshoot 1 and skip it)
      for (let stepIndex = 0; stepIndex <= eikonalExteralIterations; stepIndex++) {
        eikonalActivationFlag = stepIndex / eikonalExteralIterations;

        // Update the context object with current eikonalActivationFlag
        context.eikonalActivationFlag = eikonalActivationFlag;

        // Pass the previous solution as initial guess
        if (solutionVector.length > 0) {
          context.initialSolution = [...solutionVector];
        }

        // Solve the assembled non-linear system
        const newtonRaphsonResult = newtonRaphson(assembleFrontPropagationMat, context);

        // Extract results
        jacobianMatrix = newtonRaphsonResult.jacobianMatrix;
        residualVector = newtonRaphsonResult.residualVector;
        solutionVector = newtonRaphsonResult.solutionVector;

        // Stop at the first step that does not converge, instead of continuing from a wrong solution
        if (!newtonRaphsonResult.converged) {
          errorLog(
            `Front propagation did not converge at eikonal activation ${eikonalActivationFlag} (step ${stepIndex} of ${eikonalExteralIterations}). The solution is not valid; a finer mesh may be required.`,
          );
          break;
        }
      }
    } else if (this.solverConfig === "generalFormPDEScript") {
      // Check if using frontal solver
      if (this.solverMethod === "frontal" && this.nonlinear) {
        errorLog(
          "Frontal solver is not yet supported for nonlinear generalFormPDEScript. Please use 'lusolve' or 'jacobi'.",
        );
      } else if (this.solverMethod === "frontal") {
        const genericBoundaryConditions = new GenericBoundaryConditions(
          this.boundaryConditions,
          meshData.boundaryElements,
          meshData.nop,
          meshData.meshDimension,
          meshData.elementOrder,
          meshData.elementShape,
        );
        const frontalResult = runFrontalSolver(
          (elementArgs) =>
            assembleGeneralFormPDEFront({
              ...elementArgs,
              coefficientFunctions: this.coefficientFunctions,
            }),
          meshData,
          genericBoundaryConditions,
        );
        solutionVector = frontalResult.solutionVector;
      } 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(
          meshData,
          this.boundaryConditions,
          this.coefficientFunctions,
        ));

        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
          maxIterations: options.maxIterations ?? this.maxIterations,
          tolerance: options.tolerance ?? this.tolerance,
        });
        solutionVector = linearSystemResult.solutionVector;
      }
    } else if (this.solverConfig === "creepingFlowScript") {
      // Use regular linear solver methods for steady creeping flow
      const creepingFlowResult = assembleCreepingFlowMatrix(
        meshData,
        this.boundaryConditions,
        this.coefficientFunctions,
      );
      jacobianMatrix = creepingFlowResult.jacobianMatrix;
      residualVector = creepingFlowResult.residualVector;

      const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
        maxIterations: options.maxIterations ?? this.maxIterations,
        tolerance: options.tolerance ?? this.tolerance,
      });
      solutionVector = linearSystemResult.solutionVector;

      // Split the solution vector into velocity and pressure fields
      fields = extractCreepingFlowFields(solutionVector, meshData, 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;

      // Split the solution vector into deflection and rotation fields (2 DOFs per node)
      fields = extractEulerBernoulliBeamFields(solutionVector, beamResult.dofsPerNode);
    }

    console.timeEnd("totalSolvingTime");
    basicLog("Solving process completed");

    return fields ? { solutionVector, nodesCoordinates, fields } : { solutionVector, nodesCoordinates };
  }

  /**
   * Method to solve the finite element problem asynchronously
   * @param {object} [computeEngine] - An initialized WebGPU compute engine to reuse across solves; when omitted
   *  or null, a temporary one is created for this solve and terminated afterwards. A provided engine is not
   *  terminated, so the caller remains responsible for it
   * @param {object} [options] - Additional parameters for the solver, such as `maxIterations` and `tolerance`
   * @returns {Promise<object>} A promise that resolves to an object containing the solution vector (Float64Array) and the coordinates of the mesh nodes
   */
  async solveAsync(computeEngine, options = {}) {
    if (!this.solverConfig || !this.meshConfig || !this.boundaryConditions) {
      errorLog("Solver config, mesh config, and boundary conditions must be set before solving.");
    }
    this.validateBoundaryConditions();

    let jacobianMatrix = [];
    let residualVector = [];
    let solutionVector = new Float64Array(0);

    basicLog("Preparing mesh...");
    const meshData = prepareMesh(this.meshConfig);
    basicLog("Mesh preparation completed");
    const nodesCoordinates = {
      nodesXCoordinates: meshData.nodesXCoordinates,
      nodesYCoordinates: meshData.nodesYCoordinates,
    };

    basicLog("Beginning solving process...");
    console.time("totalSolvingTime");

    basicLog(`Using solver: ${this.solverConfig}`);
    if (this.solverConfig === "heatConductionScript") {
      ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(
        meshData,
        this.boundaryConditions,
        this.coefficientFunctions,
      ));

      if (this.solverMethod === "jacobi-gpu") {
        const { solutionVector: x } = await solveLinearSystemAsync(
          "jacobi-gpu",
          jacobianMatrix,
          residualVector,
          {
            computeEngine,
            maxIterations: options.maxIterations ?? this.maxIterations,
            tolerance: options.tolerance ?? this.tolerance,
          },
        );
        solutionVector = x;
      } else {
        // Any other async solver
      }
    }
    console.timeEnd("totalSolvingTime");
    basicLog("Solving process completed");

    return { solutionVector, nodesCoordinates };
  }
}