Feature/spatially varying coefficients

README.md

@@ -20,6 +20,7 @@
 - <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>Spatially varying coefficients:</b> `thermalConductivity(x, y)` and `heatSource(x, y)` can be scalars or functions, evaluated at each Gauss point
 - <b>Performance:</b> web worker support for multi-threaded computation
 - <b>Visualization:</b> interactive rendering with vtk.js and Plotly
 

examples/heatConductionScript/README.md

@@ -28,6 +28,37 @@ Implementation using a Gmsh-generated mesh for a rhomboid domain (the mesh file,
 
 For detailed information on the model setup, refer to the corresponding [tutorial](https://feascript.com/tutorials/heat-conduction-2d-rhom-fin-gmsh.html) in the FEAScript website.
 
+#### 5. Heat Conduction in a 1D Bi-Material Wall with Spatially Varying k(x) (`heatConduction1DVaryingK.js`)
+
+Demonstrates passing `thermalConductivity` as a function of position `k(x)`. The wall consists of two layers with different conductivities separated at mid-length. This example exercises the `coefficientFunctions` API with the standard linear solver.
+
+#### 6. Heat Conduction in a 2D Fin with Spatially Varying k(x,y) and Q(x,y) (`heatConduction2DVaryingKQ.js`)
+
+Demonstrates both `thermalConductivity(x, y)` and `heatSource(x, y)` as functions of position on a 2D structured mesh. The domain is split into a high-conductivity metal half and a low-conductivity ceramic half, with a localised volumetric heat source in the upper strip.
+
+## Spatially Varying Coefficients
+
+Both `thermalConductivity` and `heatSource` can be provided either as constants (scalars) or as functions of the physical coordinates. They are evaluated at each Gauss point during the isoparametric mapping loop, so any piecewise or smooth spatial variation is fully supported.
+
+```javascript
+model.setModelConfig("heatConductionScript", {
+  coefficientFunctions: {
+    // Scalar (uniform)
+    thermalConductivity: 10,
+    // Function of x only (1D or 2D)
+    // thermalConductivity: (x) => x < 0.5 ? 10 : 1,
+    // Function of x and y (2D)
+    // thermalConductivity: (x, y) => x < 2.0 ? 10 : 1,
+    // Uniform heat source
+    // heatSource: 500,
+    // Localised heat source (active only in upper strip)
+    heatSource: (x, y) => (y > 1.5 ? 500 : 0),
+  },
+});
+```
+
+When `coefficientFunctions` is omitted or a coefficient is not provided, the defaults `thermalConductivity = 1` and `heatSource = 0` are used. Both the standard matrix assembler and the frontal solver assembler support this feature.
+
 ## Running the Node.js Examples
 
 #### 1. Create package.json with ES module support:

examples/heatConductionScript/heatConduction1DWall/heatConduction1DVaryingK.js

@@ -0,0 +1,61 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+/**
+ * Heat Conduction in a 1D Bi-Material Wall with Spatially Varying k(x)
+ *
+ * Domain: 0 ≤ x ≤ 0.15 m
+ * Material interface at x = 0.075 m:
+ *   - Left half (x < 0.075 m):  k = 10 W/(m·K)  (e.g. concrete)
+ *   - Right half (x ≥ 0.075 m): k =  1 W/(m·K)  (e.g. insulation)
+ * Boundary conditions:
+ *   - Left  (x = 0):    convection, h = 25 W/(m²·K), T_inf = 5 °C
+ *   - 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";
+
+console.log("FEAScript Version:", printVersion);
+
+// Create a new FEAScript model
+const model = new FEAScriptModel();
+
+// Select physics/PDE
+model.setModelConfig("heatConductionScript", {
+  coefficientFunctions: {
+    // Bi-material: high-conductivity concrete | low-conductivity insulation
+    thermalConductivity: (x) => (x < 0.075 ? 10 : 1),
+    heatSource: 0,
+  },
+});
+
+// Define mesh configuration
+model.setMeshConfig({
+  meshDimension: "1D",
+  elementOrder: "linear",
+  numElementsX: 20,
+  maxX: 0.15,
+});
+
+// Define boundary conditions
+model.addBoundaryCondition("0", ["convection", 25, 5]); // Left boundary: convection
+model.addBoundaryCondition("1", ["constantTemp", 20]); // Right boundary: fixed temperature
+
+// Solve the problem
+const { solutionVector, nodesCoordinates } = model.solve();
+
+// Print results
+console.log(`Number of nodes: ${nodesCoordinates.nodesXCoordinates.length}`);
+console.log("Node coordinates (x):", nodesCoordinates.nodesXCoordinates);
+console.log("Temperature solution:", solutionVector);

examples/heatConductionScript/heatConduction2DFin/heatConduction2DVaryingKQ.js

@@ -0,0 +1,74 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Lightweight Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+/**
+ * Heat Conduction in a 2D Fin with Spatially Varying k(x,y) and Q(x,y)
+ *
+ * Domain: 0 ≤ x ≤ 4 m, 0 ≤ y ≤ 2 m
+ *
+ * Thermal conductivity varies by material region:
+ *   - Left half  (x < 2): k = 10 W/(m·K)  (high-conductivity metal)
+ *   - Right half (x ≥ 2): k =  1 W/(m·K)  (low-conductivity ceramic)
+ *
+ * Heat source active only in the upper strip (y > 1.5):
+ *   - Q = 500 W/m³  for y > 1.5
+ *   - Q = 0         otherwise
+ *
+ * Boundary conditions:
+ *   - Bottom (y = 0):  constant temperature, T = 200 °C  (heated base)
+ *   - Left   (x = 0):  symmetry
+ *   - Top    (y = 2):  convection, h = 1 W/(m²·K), T_inf = 20 °C
+ *   - 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";
+
+console.log("FEAScript Version:", printVersion);
+
+// Create a new FEAScript model
+const model = new FEAScriptModel();
+
+// Select physics/PDE with spatially varying coefficients
+model.setModelConfig("heatConductionScript", {
+  coefficientFunctions: {
+    // Bi-material fin: high-k metal on left, low-k ceramic on right
+    thermalConductivity: (x, y) => (x < 2.0 ? 10 : 1),
+    // Localised heat source in the upper strip
+    heatSource: (x, y) => (y > 1.5 ? 500 : 0),
+  },
+});
+
+// Define mesh configuration
+model.setMeshConfig({
+  meshDimension: "2D",
+  elementOrder: "quadratic",
+  numElementsX: 8,
+  numElementsY: 4,
+  maxX: 4,
+  maxY: 2,
+});
+
+// Define boundary conditions
+model.addBoundaryCondition("0", ["constantTemp", 200]); // Bottom boundary
+model.addBoundaryCondition("1", ["symmetry"]); // Left boundary
+model.addBoundaryCondition("2", ["convection", 1, 20]); // Top boundary
+model.addBoundaryCondition("3", ["constantTemp", 200]); // Right boundary
+
+// Solve the problem
+const { solutionVector, nodesCoordinates } = model.solve();
+
+// Print results
+console.log(`Number of nodes in mesh: ${nodesCoordinates.nodesXCoordinates.length}`);
+console.log("Node coordinates:", nodesCoordinates);
+console.log("Temperature solution:", solutionVector);

src/FEAScript.js

@@ -81,6 +81,19 @@ export class FEAScriptModel {
   }
 
   addBoundaryCondition(boundaryKey, condition) {
+    // 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]}`);
   }
@@ -136,11 +149,16 @@ export class FEAScriptModel {
           assembleHeatConductionFront,
           meshData,
           this.boundaryConditions,
+          { coefficientFunctions: this.coefficientFunctions },
         );
         solutionVector = frontalResult.solutionVector;
       } else {
         // Use regular linear solver methods
-        ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(meshData, this.boundaryConditions));
+        ({ 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,
@@ -303,7 +321,11 @@ export class FEAScriptModel {
 
     basicLog(`Using solver: ${this.solverConfig}`);
     if (this.solverConfig === "heatConductionScript") {
-      ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(meshData, this.boundaryConditions));
+      ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(
+        meshData,
+        this.boundaryConditions,
+        this.coefficientFunctions,
+      ));
 
       if (this.solverMethod === "jacobi-gpu") {
         const { solutionVector: x } = await solveLinearSystemAsync(

src/methods/frontalSolver.js

@@ -113,6 +113,7 @@ export function runFrontalSolver(assembleFront, meshData, boundaryConditions, op
   // Parameters for non-linear assemblers
   frontalState.currentSolutionVector = options.solutionVector;
   frontalState.eikonalActivationFlag = options.eikonalActivationFlag;
+  frontalState.coefficientFunctions = options.coefficientFunctions;
 
   // Pass assembleFront and dirichletBoundaryConditionsHandler to runFrontalAlgorithm
   runFrontalAlgorithm(meshData, FEAData, dirichletBoundaryConditionsHandler, assembleFront);
@@ -244,6 +245,7 @@ function assembleElementContribution(meshData, FEAData, thermalBoundaryCondition
     // These are ignored by linear assemblers
     solutionVector: frontalState.currentSolutionVector,
     eikonalActivationFlag: frontalState.eikonalActivationFlag,
+    coefficientFunctions: frontalState.coefficientFunctions,
   });
 
   // Handle Robin-type boundary conditions differently based on which solver is being used

src/models/creepingFlow.js

@@ -183,13 +183,13 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
         // Assemble viscous stiffness terms (K block)
         for (let localNodeIndex1 = 0; localNodeIndex1 < nodesPerVelocityElement; localNodeIndex1++) {
           let globalNode1 = velLocalToGlobalMap[localNodeIndex1];
-          let uDOF1 = globalNode1; // u-velocity DOF
-          let vDOF1 = totalNodesVelocity + globalNode1; // v-velocity DOF
+          let xVelocityDegreeOfFreedom1 = globalNode1; // u-velocity DOF
+          let yVelocityDegreeOfFreedom1 = totalNodesVelocity + globalNode1; // v-velocity DOF
 
           for (let localNodeIndex2 = 0; localNodeIndex2 < nodesPerVelocityElement; localNodeIndex2++) {
             let globalNode2 = velLocalToGlobalMap[localNodeIndex2];
-            let uDOF2 = globalNode2; // u-velocity DOF
-            let vDOF2 = totalNodesVelocity + globalNode2; // v-velocity DOF
+            let xVelocityDegreeOfFreedom2 = globalNode2; // u-velocity DOF
+            let yVelocityDegreeOfFreedom2 = totalNodesVelocity + globalNode2; // v-velocity DOF
 
             // Viscous stiffness
             let viscousContribution =
@@ -199,13 +199,13 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
                 basisFunctionDerivY[localNodeIndex1] * basisFunctionDerivY[localNodeIndex2]);
 
             // K appears in both u-u and v-v blocks
-            jacobianMatrix[uDOF1][uDOF2] += viscousContribution;
-            jacobianMatrix[vDOF1][vDOF2] += viscousContribution;
+            jacobianMatrix[xVelocityDegreeOfFreedom1][xVelocityDegreeOfFreedom2] += viscousContribution;
+            jacobianMatrix[yVelocityDegreeOfFreedom1][yVelocityDegreeOfFreedom2] += viscousContribution;
           }
 
           // Assemble pressure-velocity coupling terms
           for (let localPresIndex = 0; localPresIndex < nodesPerPressureElement; localPresIndex++) {
-            let pDOF = 2 * totalNodesVelocity + presLocalToGlobalMap[localPresIndex];
+            let pressureDegreeOfFreedom = 2 * totalNodesVelocity + presLocalToGlobalMap[localPresIndex];
 
             let bxContribution =
               weightFactor *
@@ -218,14 +218,14 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
               basisFunctionDerivY[localNodeIndex1];
 
             // Pressure gradient in x-momentum
-            jacobianMatrix[uDOF1][pDOF] += bxContribution;
+            jacobianMatrix[xVelocityDegreeOfFreedom1][pressureDegreeOfFreedom] += bxContribution;
 
             // Pressure gradient in y-momentum
-            jacobianMatrix[vDOF1][pDOF] += byContribution;
+            jacobianMatrix[yVelocityDegreeOfFreedom1][pressureDegreeOfFreedom] += byContribution;
 
             // Continuity equation
-            jacobianMatrix[pDOF][uDOF1] += -bxContribution;
-            jacobianMatrix[pDOF][vDOF1] += -byContribution;
+            jacobianMatrix[pressureDegreeOfFreedom][xVelocityDegreeOfFreedom1] += -bxContribution;
+            jacobianMatrix[pressureDegreeOfFreedom][yVelocityDegreeOfFreedom1] += -byContribution;
           }
         }
       }