From a1ef9e69b6e56aad1dfcb54d2e03eba30f2a12df Mon Sep 17 00:00:00 2001 From: nikoscham Date: Wed, 28 Oct 2026 20:14:56 +0200 Subject: [PATCH] Fix non-contiguous node tags in the Gmsh reader (#96) - Gmsh reader: when node tags have gaps (max tag > number of nodes), map each tag to a compact 1-based index and remap element and boundary node references; meshes with contiguous tags are read as before - Previously each gap added a disconnected node with an all-zero matrix row, so Jacobi returned NaN and the LU solver failed - Tests: add gmshReader unit test (test meshes copied with gaps in their node tags give the same nodes, connectivity and solution) --- src/readers/gmshReader.js | 23 +++- tests/unit/README.md | 1 + tests/unit/gmshReader.test.js | 223 ++++++++++++++++++++++++++++++++++ 3 files changed, 243 insertions(+), 4 deletions(-) create mode 100644 tests/unit/gmshReader.test.js diff --git a/src/readers/gmshReader.js b/src/readers/gmshReader.js index 6da405d..5f1a792 100644 --- a/src/readers/gmshReader.js +++ b/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; diff --git a/tests/unit/README.md b/tests/unit/README.md index 601ec9d..10ecb6b 100644 --- a/tests/unit/README.md +++ b/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: diff --git a/tests/unit/gmshReader.test.js b/tests/unit/gmshReader.test.js new file mode 100644 index 0000000..74476d4 --- /dev/null +++ b/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);