Return named result fields for creeping flow and the beam (#106)

examples/creepingFlowScript/lidDrivenCavity2DCreepingFlow/lidDrivenCavity2DCreepingFlow.js

@@ -34,9 +34,11 @@ model.addBoundaryCondition("3", ["constantVelocity", 0, 0]); // Right
 
 // 4. Solve
 model.setSolverMethod("lusolve");
-const { solutionVector, nodesCoordinates } = model.solve();
+const { nodesCoordinates, fields } = model.solve();
 
-// 5. Print the results
+// 5. Print the results (one value per mesh node for each field)
 console.log(`Number of nodes in mesh: ${nodesCoordinates.nodesXCoordinates.length}`);
 console.log("Node coordinates:", nodesCoordinates);
-console.log("Solution vector:", solutionVector);
+console.log("x-velocity u:", fields.u);
+console.log("y-velocity v:", fields.v);
+console.log("Pressure p:", fields.p);

examples/eulerBernoulliBeamScript/clampedSpringSupportedBeam1D/clampedSpringSupportedBeam1D.js

@@ -39,20 +39,16 @@ model.addBoundaryCondition("3", [
 
 // 4. Solve
 model.setSolverMethod("lusolve");
-const { solutionVector } = model.solve();
+const { nodesCoordinates, fields } = model.solve();
 
 // 5. Print the results
-const flatSolution = solutionVector.map((entry) => (Array.isArray(entry) ? entry[0] : entry));
-
-const nodeXCoordinates = [0, 5, 10];
+const { nodesXCoordinates } = nodesCoordinates;
 console.log("\nNode |    x (m) | Deflection w (m) | Rotation theta (rad)");
 console.log("-----|----------|-------------------|----------------------");
-for (let nodeIndex = 0; nodeIndex < nodeXCoordinates.length; nodeIndex++) {
-  const w = flatSolution[2 * nodeIndex];
-  const theta = flatSolution[2 * nodeIndex + 1];
+for (let nodeIndex = 0; nodeIndex < nodesXCoordinates.length; nodeIndex++) {
   console.log(
-    `  ${nodeIndex + 1}  | ${nodeXCoordinates[nodeIndex].toFixed(2).padStart(8)} | ${w
+    `  ${nodeIndex + 1}  | ${nodesXCoordinates[nodeIndex].toFixed(2).padStart(8)} | ${fields.w[nodeIndex]
       .toExponential(4)
-      .padStart(17)} | ${theta.toExponential(4).padStart(20)}`,
+      .padStart(17)} | ${fields.theta[nodeIndex].toExponential(4).padStart(20)}`,
   );
 }

src/FEAScript.js

@@ -21,8 +21,11 @@ import {
 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 { 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";
 
@@ -162,7 +165,9 @@ export class FEAScriptModel {
   /**
    * 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
+   * @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) {
@@ -183,6 +188,7 @@ export class FEAScriptModel {
     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...");
@@ -368,12 +374,8 @@ export class FEAScriptModel {
       });
       solutionVector = linearSystemResult.solutionVector;
 
-      // Store creeping-flow-specific metadata for solution extraction
-      this._creepingFlowMetadata = {
-        totalNodesVelocity: creepingFlowResult.totalNodesVelocity,
-        totalNodesPressure: creepingFlowResult.totalNodesPressure,
-        pressureNodeIndices: creepingFlowResult.pressureNodeIndices,
-      };
+      // 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(
@@ -390,17 +392,14 @@ export class FEAScriptModel {
       });
       solutionVector = linearSystemResult.solutionVector;
 
-      // Store beam-specific metadata for solution extraction (2 DOFs per node: deflection, rotation)
-      this._eulerBernoulliBeamMetadata = {
-        dofsPerNode: beamResult.dofsPerNode,
-        totalNodesX: meshData.totalNodesX,
-      };
+      // 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 { solutionVector, nodesCoordinates };
+    return fields ? { solutionVector, nodesCoordinates, fields } : { solutionVector, nodesCoordinates };
   }
 
   /**

src/physics/creepingFlow.js

@@ -274,3 +274,65 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
     pressureNodeIndices,
   };
 }
+
+/**
+ * Function to split the creeping flow solution vector into named nodal fields
+ *
+ * The pressure is solved only at the corner nodes (Q1/P1). It is interpolated linearly to the
+ * mid-side (and, for quadrilaterals, centre) nodes, so that every field has one value per mesh node
+ *
+ * @param {Float64Array} solutionVector - The solution vector [u_0 … u_{N_2−1}, v_0 … v_{N_2−1}, p_0 … p_{N_1−1}]
+ * @param {object} meshData - Object containing prepared mesh data
+ * @param {number[]} pressureNodeIndices - Array mapping pressure DOF index to global Q2 node index
+ * @returns {object} An object containing one Float64Array per field, each with one value per mesh node:
+ *  - velocityMagnitude: The velocity magnitude sqrt(u² + v²)
+ *  - u: The x-velocity
+ *  - v: The y-velocity
+ *  - p: The pressure
+ */
+export function extractCreepingFlowFields(solutionVector, meshData, pressureNodeIndices) {
+  const { nop, totalNodes, elementShape } = meshData;
+
+  const u = Float64Array.from(solutionVector.slice(0, totalNodes));
+  const v = Float64Array.from(solutionVector.slice(totalNodes, 2 * totalNodes));
+  const velocityMagnitude = u.map((uValue, nodeIndex) => Math.hypot(uValue, v[nodeIndex]));
+
+  // Place the pressure at the corner nodes
+  const p = new Float64Array(totalNodes);
+  pressureNodeIndices.forEach((globalNodeIndex, pressureIndex) => {
+    p[globalNodeIndex] = solutionVector[2 * totalNodes + pressureIndex];
+  });
+
+  // Interpolate the pressure to the remaining nodes as the average of the corner nodes they lie between
+  // Q2 local numbering:    P2 local numbering:
+  //   2 - 5 - 8            2
+  //   |       |            | \
+  //   1   4   7            1   4
+  //   |       |            |     \
+  //   0 - 3 - 6            0 - 3 - 5
+  const interpolatedNodes =
+    elementShape === "triangle"
+      ? [
+          [1, [0, 2]],
+          [3, [0, 5]],
+          [4, [2, 5]],
+        ]
+      : [
+          [1, [0, 2]],
+          [3, [0, 6]],
+          [5, [2, 8]],
+          [7, [6, 8]],
+          [4, [0, 2, 6, 8]],
+        ];
+  for (const elementNodes of nop) {
+    for (const [localNodeIndex, cornerLocalIndices] of interpolatedNodes) {
+      const cornerSum = cornerLocalIndices.reduce(
+        (sum, cornerLocalIndex) => sum + p[elementNodes[cornerLocalIndex] - 1],
+        0,
+      );
+      p[elementNodes[localNodeIndex] - 1] = cornerSum / cornerLocalIndices.length;
+    }
+  }
+
+  return { velocityMagnitude, u, v, p };
+}

src/physics/eulerBernoulliBeam.js

@@ -174,3 +174,17 @@ export function assembleEulerBernoulliBeamMat(meshData, boundaryConditions, coef
     totalDOFs,
   };
 }
+
+/**
+ * Function to split the Euler-Bernoulli beam solution vector into named nodal fields
+ * @param {Float64Array} solutionVector - The solution vector [w_0, theta_0, w_1, theta_1, …]
+ * @param {number} dofsPerNode - Number of degrees of freedom per node (2)
+ * @returns {object} An object containing one Float64Array per field, each with one value per mesh node:
+ *  - w: The deflection
+ *  - theta: The rotation (dw/dx)
+ */
+export function extractEulerBernoulliBeamFields(solutionVector, dofsPerNode) {
+  const w = Float64Array.from(solutionVector.filter((_, dofIndex) => dofIndex % dofsPerNode === 0));
+  const theta = Float64Array.from(solutionVector.filter((_, dofIndex) => dofIndex % dofsPerNode === 1));
+  return { w, theta };
+}

src/visualization/plotlyPlot.js

@@ -11,18 +11,23 @@
 import { prepareMesh, computeNodeNeighbors, getBoundarySegments } from "../mesh/meshUtils.js";
 import { pointInsideQuadrilateral, pointInsideTriangularElement } from "../elements/elementGeometry.js";
 import { BasisFunctions } from "../elements/basisFunctions.js";
+import { selectSolutionField } from "./solutionFields.js";
 
 /**
  * Function to create plots of the solution vector
  * @param {object} result - Object containing solution vector and mesh information
  * @param {object} model - Object containing model properties
  * @param {string} plotType - The type of plot
  * @param {string} plotDivId - The id of the div where the plot will be rendered
+ * @param {object} [options] - Plot options:
+ *  - field: The named field to plot, for models with more than one unknown per node (e.g. "u", "v", "p" for
+ *    creeping flow, "w", "theta" for the Euler-Bernoulli beam); defaults to the main field of the model
  */
-export function plotSolution(model, result, plotType, plotDivId) {
+export function plotSolution(model, result, plotType, plotDivId, options = {}) {
   console.time("plottingTime");
   const { nodesXCoordinates, nodesYCoordinates } = result.nodesCoordinates;
-  const solutionVector = result.solutionVector;
+  const { values: solutionVector, fieldName } = selectSolutionField(model, result, options.field);
+  const solutionLabel = fieldName ?? "Solution";
   const solverConfig = model.solverConfig;
   const meshDimension = model.meshConfig.meshDimension;
 
@@ -35,7 +40,7 @@ export function plotSolution(model, result, plotType, plotDivId) {
       mode: "lines",
       type: "scatter",
       line: { color: "rgb(219, 64, 82)", width: 2 },
-      name: "Solution",
+      name: solutionLabel,
     };
 
     let maxWindowWidth = Math.min(window.innerWidth, 700);
@@ -47,7 +52,7 @@ export function plotSolution(model, result, plotType, plotDivId) {
       width: plotWidth,
       height: plotHeight,
       xaxis: { title: "x" },
-      yaxis: { title: "Solution" },
+      yaxis: { title: solutionLabel },
       margin: { l: 50, r: 50, t: 50, b: 50 },
     };
 
@@ -98,7 +103,7 @@ export function plotSolution(model, result, plotType, plotDivId) {
       },
       //colorscale: 'Viridis',
       colorbar: {
-        title: "Solution",
+        title: solutionLabel,
       },
       name: "Solution Field",
     };
@@ -114,9 +119,13 @@ export function plotSolution(model, result, plotType, plotDivId) {
  * @param {object} model - Object containing model properties
  * @param {string} plotType - The type of plot
  * @param {string} plotDivId - The id of the div where the plot will be rendered
+ * @param {object} [options] - Plot options:
+ *  - field: The named field to plot, for models with more than one unknown per node; defaults to the main field of the model
  */
-export function plotInterpolatedSolution(model, result, plotType, plotDivId) {
+export function plotInterpolatedSolution(model, result, plotType, plotDivId, options = {}) {
   console.time("plottingTime");
+  const { values, fieldName } = selectSolutionField(model, result, options.field);
+  result = { ...result, solutionVector: values, fields: undefined }; // Interpolate the selected field
   const { nodesXCoordinates, nodesYCoordinates } = result.nodesCoordinates; // TODO: Check if we should place it inside the 2D block
   const meshDimension = model.meshConfig.meshDimension;
   const meshData = prepareMesh(model.meshConfig); // Retrieve mesh connectivity details
@@ -275,7 +284,7 @@ export function plotInterpolatedSolution(model, result, plotType, plotDivId) {
       },
       //colorscale: 'Viridis',
       colorbar: {
-        title: "Solution",
+        title: fieldName ?? "Solution",
       },
       name: "Interpolated Solution Field",
     };

src/visualization/solutionFields.js

@@ -0,0 +1,50 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+// Internal imports
+import { errorLog } from "../utilities/logging.js";
+
+// Field plotted when none is requested, for models that return named fields
+const defaultFields = {
+  creepingFlowScript: "velocityMagnitude",
+  eulerBernoulliBeamScript: "w",
+};
+
+/**
+ * Function to select the nodal values to plot from a solve result
+ * @param {object} model - Object containing model properties
+ * @param {object} result - Object containing the solution vector, mesh information and, for models with more than
+ *  one unknown per node, the named `fields`
+ * @param {string} [field] - The name of the field to plot (e.g. "u", "v", "p", "velocityMagnitude", "w", "theta");
+ *  when omitted, the default field of the model is used
+ * @returns {object} An object containing:
+ *  - values: The nodal values to plot (one value per mesh node)
+ *  - fieldName: The name of the plotted field, or null when the solution vector itself is plotted
+ */
+export function selectSolutionField(model, result, field) {
+  if (!result.fields) {
+    if (field) {
+      errorLog(
+        `Field "${field}" was requested, but this result has no named fields; the solution vector is plotted.`,
+      );
+    }
+    return { values: result.solutionVector, fieldName: null };
+  }
+
+  const defaultField = defaultFields[model.solverConfig] ?? Object.keys(result.fields)[0];
+  if (field && !(field in result.fields)) {
+    errorLog(
+      `Unknown field "${field}"; plotting "${defaultField}" instead. Available fields: ${Object.keys(
+        result.fields,
+      ).join(", ")}.`,
+    );
+  }
+  const fieldName = field in result.fields ? field : defaultField;
+  return { values: result.fields[fieldName], fieldName };
+}

src/visualization/vtkPlot.js

@@ -12,6 +12,7 @@ import { prepareMesh, computeNodeNeighbors, getBoundarySegments } from "../mesh/
 import { pointInsideQuadrilateral, pointInsideTriangularElement } from "../elements/elementGeometry.js";
 import { BasisFunctions } from "../elements/basisFunctions.js";
 import { errorLog } from "../utilities/logging.js";
+import { selectSolutionField } from "./solutionFields.js";
 
 let _vtkModules = null;
 
@@ -104,14 +105,16 @@ export function createContourLineOptions(options = {}) {
  * @param {object} result - Object containing solution vector and node coordinates
  * @param {string} plotType - The type of plot (e.g., "line", "contour")
  * @param {string} plotDivId - The id of the HTML element for rendering
- * @param {object} [renderOptions={}] - Rendering options including colorScale, contourLines
+ * @param {object} [renderOptions={}] - Rendering options including colorScale, contourLines and field (the named
+ *  field to plot, for models with more than one unknown per node)
  */
 export async function plotSolutionVtk(model, result, plotType, plotDivId, renderOptions = {}) {
   console.time("plottingTime");
   const meshDimension = model.meshConfig.meshDimension;
   const meshData = prepareMesh(model.meshConfig);
   const vtkData = await transformSolverOutputToVtkData(model, result, meshData, {
     mode: meshDimension === "1D" && plotType === "line" ? "line" : "surface",
+    field: renderOptions.field,
   });
   await renderVtkScene(vtkData, plotDivId, model.solverConfig, plotType, renderOptions);
   console.timeEnd("plottingTime");
@@ -123,10 +126,13 @@ export async function plotSolutionVtk(model, result, plotType, plotDivId, render
  * @param {object} result - Object containing solution vector and node coordinates
  * @param {string} plotType - The type of plot (e.g., "contour")
  * @param {string} plotDivId - The id of the HTML element for rendering
- * @param {object} [renderOptions={}] - Rendering options including colorScale, contourLines
+ * @param {object} [renderOptions={}] - Rendering options including colorScale, contourLines and field (the named
+ *  field to plot, for models with more than one unknown per node)
  */
 export async function plotInterpolatedSolutionVtk(model, result, plotType, plotDivId, renderOptions = {}) {
   console.time("plottingTime");
+  const { values } = selectSolutionField(model, result, renderOptions.field);
+  result = { ...result, solutionVector: values, fields: undefined }; // Plot the selected field
   if (model.meshConfig.meshDimension !== "2D" || plotType !== "contour") {
     // Fallback to non-interpolated for non-2D contour cases
     const meshData = prepareMesh(model.meshConfig);
@@ -154,13 +160,15 @@ export async function plotInterpolatedSolutionVtk(model, result, plotType, plotD
  * @param {object} model - FEAScript model containing mesh and solver configuration
  * @param {object} result - Object containing solution vector and node coordinates
  * @param {object|null} [meshData=null] - Pre-prepared mesh data, or null to auto-prepare
- * @param {object} [options={}] - Options including mode ("surface" or "line")
+ * @param {object} [options={}] - Options including mode ("surface" or "line") and field (the named field to export,
+ *  for models with more than one unknown per node)
  * @returns {Promise<object>} VTK data object with points, scalars, cells, polyData, and metadata
  */
 export async function transformSolverOutputToVtkData(model, result, meshData = null, options = {}) {
   const preparedMesh = meshData ?? prepareMesh(model.meshConfig);
   const { nodesXCoordinates, nodesYCoordinates } = result.nodesCoordinates;
-  const solutionArray = extractScalarSolution(result.solutionVector, nodesXCoordinates.length);
+  const { values } = selectSolutionField(model, result, options.field);
+  const solutionArray = extractScalarSolution(values, nodesXCoordinates.length);
   const points = buildPointsArray(nodesXCoordinates, nodesYCoordinates);
 
   const mode = options.mode ?? "surface";
@@ -191,7 +199,7 @@ export async function transformSolverOutputToVtkData(model, result, meshData = n
  * @param {object} model - FEAScript model containing mesh and solver configuration
  * @param {object} result - Object containing solution vector and node coordinates
  * @param {object|null} [meshData=null] - Pre-prepared mesh data, or null to auto-prepare
- * @param {object} [options={}] - Options including mode ("surface" or "line")
+ * @param {object} [options={}] - Options including mode ("surface" or "line") and field (the named field to export)
  * @returns {Promise<string>} VTP XML string
  */
 export async function transformSolverOutputToVTP(model, result, meshData = null, options = {}) {