Fix non-contiguous node tags in the Gmsh reader (#96)

src/readers/gmshReader.js

@@ -56,6 +56,9 @@ const importGmshMesh = async (file) => {
   let nodeTagsCollected = 0;
   let nodeCoordinatesCollected = 0;
 
+  // Node tags with gaps are mapped to compact 1-based indices (null when tags are contiguous)
+  let nodeTagToIndex = null;
+
   // Elements
   let elementEntityBlocks = 0;
   let elementBlocksProcessed = 0;
@@ -178,6 +181,13 @@ const importGmshMesh = async (file) => {
       if (nodeEntityBlocks === 0) {
         nodeEntityBlocks = parseInt(parts[0], 10);
         totalNodes = parseInt(parts[1], 10);
+        const maxNodeTag = parseInt(parts[3], 10);
+        if (maxNodeTag !== totalNodes) {
+          nodeTagToIndex = new Map();
+          debugLog(
+            `Gmsh node tags are not contiguous (max tag ${maxNodeTag}, ${totalNodes} nodes); renumbering`,
+          );
+        }
         result.nodesXCoordinates = new Array(totalNodes).fill(0);
         result.nodesYCoordinates = new Array(totalNodes).fill(0);
       } else if (nodeBlocksProcessed < nodeEntityBlocks && currentNodeBlock.numNodes === 0) {
@@ -197,9 +207,13 @@ const importGmshMesh = async (file) => {
           if (nodeTagsCollected === currentNodeBlock.numNodes) break;
         }
       } else if (nodeCoordinatesCollected < currentNodeBlock.numNodes) {
-        const nodeTag = nodeTags[nodeCoordinatesCollected] - 1;
-        result.nodesXCoordinates[nodeTag] = parseFloat(parts[0]);
-        result.nodesYCoordinates[nodeTag] = parseFloat(parts[1]);
+        let nodeIndex = nodeTags[nodeCoordinatesCollected] - 1;
+        if (nodeTagToIndex) {
+          nodeIndex = nodeTagToIndex.size;
+          nodeTagToIndex.set(nodeTags[nodeCoordinatesCollected], nodeIndex + 1);
+        }
+        result.nodesXCoordinates[nodeIndex] = parseFloat(parts[0]);
+        result.nodesYCoordinates[nodeIndex] = parseFloat(parts[1]);
         result.totalNodesX++;
         result.totalNodesY++;
         nodeCoordinatesCollected++;
@@ -227,7 +241,8 @@ const importGmshMesh = async (file) => {
 
         elementsProcessedInBlock = 0;
       } else if (elementsProcessedInBlock < currentElementBlock.numElements) {
-        const nodeIndices = parts.slice(1).map((n) => parseInt(n, 10));
+        let nodeIndices = parts.slice(1).map((n) => parseInt(n, 10));
+        if (nodeTagToIndex) nodeIndices = nodeIndices.map((tag) => nodeTagToIndex.get(tag));
 
         // Resolve physical tag via entity → physical map
         let physicalTag = currentElementBlock.tag;

tests/unit/README.md

@@ -17,6 +17,7 @@ and utilities.
 | `elements2D.test.js`         | Basis functions and quadrature rules of linear and quadratic quadrilaterals and triangles           |
 | `plotInterpolation.test.js`  | Point location and solution interpolation of the contour plot, for all 2D element types             |
 | `solutionFormat.test.js`     | Solution vector is a Float64Array for every model and solver method, with the same values and plots |
+| `gmshReader.test.js`         | Gmsh reader, including meshes with non-contiguous node tags (same nodes, connectivity, solution)    |
 
 Run all unit, regression, and verification tests from the repository root with:
 

tests/unit/gmshReader.test.js

@@ -0,0 +1,223 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 (RC) | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+/**
+ * Unit tests for the Gmsh reader (importGmshMesh)
+ *
+ * Uses the test meshes of linear and quadratic quadrilaterals and triangles
+ *
+ * Covers:
+ *  - Non-contiguous node tags: each mesh is copied with gaps in its node tags (same geometry and
+ *    connectivity) and compared against the original
+ *    - No extra (disconnected) nodes are created by the gaps
+ *    - Element and boundary node references point to the same coordinates
+ *    - The heat conduction solution (lusolve and jacobi) is the same at every element node
+ *
+ * Run: node tests/unit/gmshReader.test.js (or npm test)
+ */
+
+import fs from "fs";
+import path from "path";
+import { fileURLToPath } from "url";
+import { FEAScriptModel } from "../../src/FEAScript.js";
+import { importGmshMesh } from "../../src/readers/gmshReader.js";
+import { basicLog, errorLog } from "../../src/utilities/logging.js";
+
+const __dirname = path.dirname(fileURLToPath(import.meta.url));
+const meshesDirectory = path.join(__dirname, "../meshes");
+
+const TOLERANCE = 1e-8;
+
+basicLog("");
+basicLog("================================");
+basicLog("Unit tests: Gmsh reader");
+
+let passed = 0;
+let failed = 0;
+
+function assert(condition, message) {
+  if (!condition) {
+    errorLog(`FAIL: ${message}`);
+    failed++;
+  } else {
+    basicLog(`PASS: ${message}`);
+    passed++;
+  }
+}
+
+const meshCases = [
+  {
+    meshFile: "rect_quad_linear.msh",
+    elementOrder: "linear",
+    solverMethods: ["lusolve"],
+  },
+  {
+    meshFile: "rect_quad_quadratic.msh",
+    elementOrder: "quadratic",
+    solverMethods: ["lusolve"],
+  },
+  {
+    meshFile: "rect_triang_linear.msh",
+    elementOrder: "linear",
+    solverMethods: ["lusolve", "jacobi"],
+  },
+  {
+    meshFile: "rect_triang_quadratic.msh",
+    elementOrder: "quadratic",
+    solverMethods: ["lusolve"],
+  },
+];
+
+// Gmsh physical tags: 1 bottom, 2 right, 3 top, 4 left (stored as boundary keys 0-3)
+const boundaryConditions = {
+  0: ["constantTemperature", 200],
+  1: ["constantTemperature", 200],
+  2: ["convection", 1, 20],
+  3: ["symmetry"],
+};
+
+const readMeshFile = (meshFile) => fs.readFileSync(path.join(meshesDirectory, meshFile), "utf8");
+
+const parseMesh = (meshContent) => importGmshMesh({ text: async () => meshContent });
+
+// Coordinates of the nodes referenced by a list of 1-based node numbers
+const nodeCoordinates = (mesh, nodeNumbers) =>
+  nodeNumbers.map((node) => [mesh.nodesXCoordinates[node - 1], mesh.nodesYCoordinates[node - 1]]);
+
+const largestDifference = (a, b) => {
+  const flatB = b.flat(Infinity);
+  return Math.max(0, ...a.flat(Infinity).map((value, i) => Math.abs(value - flatB[i])));
+};
+
+function runSimulation(parsedMesh, elementOrder, solverMethod) {
+  const model = new FEAScriptModel();
+  model.setModelConfig("heatConductionScript");
+  model.setMeshConfig({ parsedMesh, meshDimension: "2D", elementOrder });
+  Object.entries(boundaryConditions).forEach(([boundaryKey, condition]) => {
+    model.addBoundaryCondition(boundaryKey, condition);
+  });
+  model.setSolverMethod(solverMethod);
+  return model.solve({ maxIterations: 200000, tolerance: 1e-10 });
+}
+
+// ─── Non-contiguous node tags ───────────────────────────────────
+
+// Node tags above TAG_GAP_START are shifted by TAG_GAP, leaving tags that no node uses
+const TAG_GAP_START = 10;
+const TAG_GAP = 100;
+
+/**
+ * Function to copy a Gmsh (.msh v4.1 ASCII) file with gaps in its node tags
+ * @param {string} meshContent - Content of the original mesh file
+ * @returns {string} Content of the mesh file with shifted node tags
+ */
+function addNodeTagGaps(meshContent) {
+  const shiftTag = (tag) => (tag > TAG_GAP_START ? tag + TAG_GAP : tag);
+  const lines = meshContent.split("\n");
+  const output = [];
+  let lineIndex = 0;
+  const nextParts = () => lines[lineIndex++].trim().split(/\s+/).map(Number);
+
+  while (lineIndex < lines.length) {
+    const line = lines[lineIndex++];
+    output.push(line);
+
+    if (line.trim() === "$Nodes") {
+      const [numBlocks, numNodes, minTag, maxTag] = nextParts();
+      output.push(`${numBlocks} ${numNodes} ${minTag} ${shiftTag(maxTag)}`);
+      for (let block = 0; block < numBlocks; block++) {
+        const blockHeader = nextParts();
+        output.push(blockHeader.join(" "));
+        const numBlockNodes = blockHeader[3];
+        for (let i = 0; i < numBlockNodes; i++) output.push(String(shiftTag(nextParts()[0])));
+        for (let i = 0; i < numBlockNodes; i++) output.push(lines[lineIndex++]);
+      }
+    } else if (line.trim() === "$Elements") {
+      const elementsHeader = nextParts();
+      output.push(elementsHeader.join(" "));
+      for (let block = 0; block < elementsHeader[0]; block++) {
+        const blockHeader = nextParts();
+        output.push(blockHeader.join(" "));
+        for (let i = 0; i < blockHeader[3]; i++) {
+          const [elementTag, ...nodeTags] = nextParts();
+          output.push([elementTag, ...nodeTags.map(shiftTag)].join(" "));
+        }
+      }
+    }
+  }
+
+  return output.join("\n");
+}
+
+for (const { meshFile, elementOrder, solverMethods } of meshCases) {
+  const meshContent = readMeshFile(meshFile);
+  const original = await parseMesh(meshContent);
+  const gapped = await parseMesh(addNodeTagGaps(meshContent));
+
+  const originalElements = [
+    ...original.nodalNumbering.quadElements,
+    ...original.nodalNumbering.triangleElements,
+  ];
+  const gappedElements = [...gapped.nodalNumbering.quadElements, ...gapped.nodalNumbering.triangleElements];
+
+  const numNodes = original.nodesXCoordinates.length;
+  const definedNodes = gapped.nodesXCoordinates.filter((x) => x !== undefined).length;
+  assert(
+    gapped.nodesXCoordinates.length === numNodes && definedNodes === numNodes,
+    `${meshFile}: ${numNodes} nodes with gapped tags (got ${gapped.nodesXCoordinates.length}, ${definedNodes} defined)`,
+  );
+
+  const elementDifference = largestDifference(
+    gappedElements.map((element) => nodeCoordinates(gapped, element)),
+    originalElements.map((element) => nodeCoordinates(original, element)),
+  );
+  assert(
+    elementDifference === 0,
+    `${meshFile}: element nodes have the same coordinates (largest difference ${elementDifference})`,
+  );
+
+  const boundaryTags = Object.keys(original.boundaryNodePairs);
+  const boundaryDifference = largestDifference(
+    boundaryTags.map((tag) => gapped.boundaryNodePairs[tag].map((pair) => nodeCoordinates(gapped, pair))),
+    boundaryTags.map((tag) => original.boundaryNodePairs[tag].map((pair) => nodeCoordinates(original, pair))),
+  );
+  assert(
+    boundaryDifference === 0,
+    `${meshFile}: boundary nodes have the same coordinates (largest difference ${boundaryDifference})`,
+  );
+
+  // The simulation modifies the parsed mesh, so each one parses its own copy
+  for (const solverMethod of solverMethods) {
+    const originalSolution = runSimulation(
+      await parseMesh(meshContent),
+      elementOrder,
+      solverMethod,
+    ).solutionVector;
+    const gappedSolution = runSimulation(
+      await parseMesh(addNodeTagGaps(meshContent)),
+      elementOrder,
+      solverMethod,
+    ).solutionVector;
+
+    const allFinite = Array.from(gappedSolution).every(Number.isFinite);
+    const solutionDifference = largestDifference(
+      gappedElements.map((element) => element.map((node) => gappedSolution[node - 1])),
+      originalElements.map((element) => element.map((node) => originalSolution[node - 1])),
+    );
+    assert(
+      allFinite && gappedSolution.length === numNodes && solutionDifference < TOLERANCE,
+      `${meshFile}, ${solverMethod}: same temperatures at the element nodes, largest difference ${solutionDifference.toExponential(
+        3,
+      )} (tolerance ${TOLERANCE})`,
+    );
+  }
+}
+
+basicLog(`Gmsh reader tests: ${passed} passed, ${failed} failed`);
+if (failed > 0) process.exit(1);