diff --git a/.prettierignore b/.prettierignore new file mode 100644 index 0000000..c72df9a --- /dev/null +++ b/.prettierignore @@ -0,0 +1,5 @@ +# Build output +dist/ + +# Third-party bundles +src/vendor/ diff --git a/.prettierrc b/.prettierrc new file mode 100644 index 0000000..3c9a3b2 --- /dev/null +++ b/.prettierrc @@ -0,0 +1,4 @@ +{ + "printWidth": 110, + "trailingComma": "all" +} diff --git a/src/FEAScript.js b/src/FEAScript.js index e32dfa3..fccc5d9 100644 --- a/src/FEAScript.js +++ b/src/FEAScript.js @@ -12,11 +12,13 @@ import { newtonRaphson } from "./solvers/newtonRaphson.js"; import { solveLinearSystem } from "./solvers/linearSystemSolver.js"; import { solveLinearSystemAsync } from "./solvers/linearSystemSolver.js"; import { prepareMesh } from "./mesh/meshUtils.js"; -import { assembleFrontPropagationMat } from "./models/frontPropagation.js"; -import { assembleGeneralFormPDEMat, assembleGeneralFormPDENonlinearMat } from "./models/generalFormPDE.js"; -import { assembleHeatConductionMat, assembleHeatConductionFront } from "./models/heatConduction.js"; -import { assembleCreepingFlowMatrix } from "./models/creepingFlow.js"; -import { assembleEulerBernoulliBeamMat } from "./models/eulerBernoulliBeam.js"; +import { assembleFrontPropagationMat, assembleFrontPropagationFront } from "./physics/frontPropagation.js"; +import { assembleGeneralFormPDEMat, 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 } from "./physics/creepingFlow.js"; +import { assembleEulerBernoulliBeamMat } from "./physics/eulerBernoulliBeam.js"; import { runFrontalSolver } from "./solvers/frontalSolver.js"; import { basicLog, debugLog, warnLog, errorLog } from "./utilities/logging.js"; @@ -141,11 +143,22 @@ export class FEAScriptModel { 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( - assembleHeatConductionFront, + (elementArgs) => + assembleHeatConductionFront({ + ...elementArgs, + coefficientFunctions: this.coefficientFunctions, + }), meshData, - this.boundaryConditions, - { coefficientFunctions: this.coefficientFunctions }, + thermalBoundaryConditions, ); solutionVector = frontalResult.solutionVector; } else { @@ -176,6 +189,20 @@ export class FEAScriptModel { // 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 diff --git a/src/mesh/basisFunctions.js b/src/elements/basisFunctions.js similarity index 84% rename from src/mesh/basisFunctions.js rename to src/elements/basisFunctions.js index cacbd54..712234a 100644 --- a/src/mesh/basisFunctions.js +++ b/src/elements/basisFunctions.js @@ -28,6 +28,53 @@ export const triangleBoundarySides = { }, }; +/** + * Local node indices on each side of the 1D reference element + */ +const lineBoundarySides = { + linear: { + 0: [0], // Node at the left side of the reference element + 1: [1], // Node at the right side of the reference element + }, + quadratic: { + 0: [0], // Node at the left side of the reference element + 1: [2], // Node at the right side of the reference element + }, +}; + +/** + * Local node indices on each side of the reference quadrilateral + * (side 0: bottom, eta = 0; side 1: left, ksi = 0; side 2: top, eta = 1; side 3: right, ksi = 1) + * Nodes are listed in order along each side + */ +const quadBoundarySides = { + linear: { + 0: [0, 2], // Nodes at the bottom side of the reference element + 1: [0, 1], // Nodes at the left side of the reference element + 2: [1, 3], // Nodes at the top side of the reference element + 3: [2, 3], // Nodes at the right side of the reference element + }, + quadratic: { + 0: [0, 3, 6], // Nodes at the bottom side of the reference element + 1: [0, 1, 2], // Nodes at the left side of the reference element + 2: [2, 5, 8], // Nodes at the top side of the reference element + 3: [6, 7, 8], // Nodes at the right side of the reference element + }, +}; + +/** + * Function to get the local node indices on each side of the reference element + * @param {string} meshDimension - The dimension of the mesh, either '1D' or '2D' + * @param {string} elementOrder - The order of elements, either 'linear' or 'quadratic' + * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' + * @returns {object} Map from side index to the local node indices on that side + */ +export function getBoundarySides(meshDimension, elementOrder, elementShape = "quad") { + if (meshDimension === "1D") return lineBoundarySides[elementOrder]; + if (elementShape === "triangle") return triangleBoundarySides[elementOrder]; + return quadBoundarySides[elementOrder]; +} + /** * Class to handle basis functions and their derivatives based on element configuration */ diff --git a/src/elements/elementGeometry.js b/src/elements/elementGeometry.js new file mode 100644 index 0000000..80dc124 --- /dev/null +++ b/src/elements/elementGeometry.js @@ -0,0 +1,129 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ +/** + * Function to test if a point is inside a triangle using barycentric coordinates, + * also returning the natural coordinates (ksi, eta) + * @param {number} x - X-coordinate of the point + * @param {number} y - Y-coordinate of the point + * @param {array} vertices - Triangle vertices [[x0,y0],[x1,y1],[x2,y2]] + * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} + */ +export function pointInsideTriangle(x, y, vertices) { + const tolerance = 1e-12; + const [v0, v1, v2] = vertices; + + const denom = (v1[1] - v2[1]) * (v0[0] - v2[0]) + (v2[0] - v1[0]) * (v0[1] - v2[1]); + + const ksi = ((v1[1] - v2[1]) * (x - v2[0]) + (v2[0] - v1[0]) * (y - v2[1])) / denom; + const eta = ((v2[1] - v0[1]) * (x - v2[0]) + (v0[0] - v2[0]) * (y - v2[1])) / denom; + const gamma = 1 - ksi - eta; + + const inside = ksi >= -tolerance && eta >= -tolerance && gamma >= -tolerance; + return { inside, ksi, eta }; +} + +/** + * Function to test if a point is inside a triangular element, also returning the natural coordinates + * (ksi, eta) of the FEAScript reference triangle + * @param {number} x - X-coordinate of the point + * @param {number} y - Y-coordinate of the point + * @param {array} vertices - Triangle vertices in FEAScript order [[x0,y0],[x1,y1],[x2,y2]] + * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} + */ +export function pointInsideTriangularElement(x, y, vertices) { + // Vertices order (reference coordinates): + // 1 (0,1) + // | \ + // | \ + // 0 --- 2 (1,0) + // The barycentric weights of vertices 2 and 1 are the natural coordinates ksi and eta + const { inside, ksi: weight0, eta: weight1 } = pointInsideTriangle(x, y, vertices); + return { inside, ksi: 1 - weight0 - weight1, eta: weight1 }; +} + +/** + * Function to test if a point is inside a quadrilateral by spliting it into triangles and using barycentric coordinates + * @param {number} x - X-coordinate of the point + * @param {number} y - Y-coordinate of the point + * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]] + * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} + */ +export function pointInsideQuadrilateral(x, y, vertices) { + const [firstTriangleVertices, secondTriangleVertices] = splitQuadrilateral(vertices); + const pointInsideFirstTriangle = pointInsideTriangle(x, y, firstTriangleVertices); + const pointInsideSecondTriangle = pointInsideTriangle(x, y, secondTriangleVertices); + + const inside = pointInsideFirstTriangle.inside || pointInsideSecondTriangle.inside; + let ksi = 0; + let eta = 0; + + if (inside) { + const [v0, v1, v2, v3] = vertices; + + // Function to calculate distance from point to line segment + const getDistanceFromLine = (p1, p2) => { + const num = Math.abs((p2[0] - p1[0]) * (p1[1] - y) - (p1[0] - x) * (p2[1] - p1[1])); + const den = Math.sqrt((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2); + return num / den; + }; + + // Calculate distances to edges based on vertex order: + // 1 (v1) --- 3 (v3) + // | | + // 0 (v0) --- 2 (v2) + + const distLeft = getDistanceFromLine(v0, v1); + const distRight = getDistanceFromLine(v2, v3); + const distBottom = getDistanceFromLine(v0, v2); + const distTop = getDistanceFromLine(v1, v3); + + ksi = distLeft / (distLeft + distRight); + eta = distBottom / (distBottom + distTop); + + // The distance ratios are exact only for parallelograms, so refine them with Newton iterations + // on the bilinear mapping x(ksi, eta) = v0 (1-ksi)(1-eta) + v1 (1-ksi) eta + v2 ksi (1-eta) + v3 ksi eta + for (let iteration = 0; iteration < 20; iteration++) { + const mappedX = + v0[0] * (1 - ksi) * (1 - eta) + v1[0] * (1 - ksi) * eta + v2[0] * ksi * (1 - eta) + v3[0] * ksi * eta; + const mappedY = + v0[1] * (1 - ksi) * (1 - eta) + v1[1] * (1 - ksi) * eta + v2[1] * ksi * (1 - eta) + v3[1] * ksi * eta; + const ksiDerivX = (v2[0] - v0[0]) * (1 - eta) + (v3[0] - v1[0]) * eta; + const ksiDerivY = (v2[1] - v0[1]) * (1 - eta) + (v3[1] - v1[1]) * eta; + const etaDerivX = (v1[0] - v0[0]) * (1 - ksi) + (v3[0] - v2[0]) * ksi; + const etaDerivY = (v1[1] - v0[1]) * (1 - ksi) + (v3[1] - v2[1]) * ksi; + const detJacobian = ksiDerivX * etaDerivY - etaDerivX * ksiDerivY; + const residualX = x - mappedX; + const residualY = y - mappedY; + const deltaKsi = (etaDerivY * residualX - etaDerivX * residualY) / detJacobian; + const deltaEta = (ksiDerivX * residualY - ksiDerivY * residualX) / detJacobian; + ksi += deltaKsi; + eta += deltaEta; + if (Math.abs(deltaKsi) + Math.abs(deltaEta) < 1e-14) break; + } + } + + return { inside, ksi, eta }; +} + +/** + * Function to split the quadrilateral elements into two triangles + * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]] + * @returns {array} Array of triangle vertices: [[v0,v1,v3], [v0,v2,v3]] + */ +export function splitQuadrilateral(vertices) { + const [v0, v1, v2, v3] = vertices; + // Vertices order: + // 1 --- 3 + // | | + // 0 --- 2 + return [ + [v0, v1, v3], + [v0, v2, v3], + ]; +} diff --git a/src/elements/feaInitialization.js b/src/elements/feaInitialization.js new file mode 100644 index 0000000..b30c657 --- /dev/null +++ b/src/elements/feaInitialization.js @@ -0,0 +1,71 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { BasisFunctions } from "./basisFunctions.js"; +import { NumericalIntegration } from "./numericalIntegration.js"; +import { getIntegrationPoints2D } from "./isoparametricMapping.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; +/** + * Function to initialize the FEA matrices and numerical tools + * @param {object} meshData - Object containing mesh data from prepareMesh() + * @returns {object} An object containing initialized matrices and numerical tools + */ +export function initializeFEA(meshData) { + const { totalNodes, nop, meshDimension, elementOrder, elementShape } = meshData; + + // Initialize variables for matrix assembly + let residualVector = new Array(totalNodes).fill(0); + let jacobianMatrix = new SparseMatrix(totalNodes); + let localToGlobalMap = []; + + // Initialize the BasisFunctions class + const basisFunctions = new BasisFunctions({ + meshDimension, + elementOrder, + elementShape, + }); + + // Initialize the NumericalIntegration class + const numericalIntegration = new NumericalIntegration({ + meshDimension, + elementOrder, + elementShape, + }); + + // Calculate Gauss points and weights + let gaussPointsAndWeights = numericalIntegration.getGaussPointsAndWeights(); + let gaussPoints = gaussPointsAndWeights.gaussPoints; + let gaussWeights = gaussPointsAndWeights.gaussWeights; + + // Integration points for 2D elements (tensor product for quads, collapsed rule for triangles) + const integrationPoints2D = + meshDimension === "2D" ? getIntegrationPoints2D(gaussPoints, gaussWeights, elementShape) : []; + + // Basis functions on the reference element at each integration point (the same for every element) + const referenceBasisFunctions = + meshDimension === "1D" + ? gaussPoints.map((gaussPoint) => basisFunctions.getBasisFunctions(gaussPoint)) + : integrationPoints2D.map(([ksi, eta]) => basisFunctions.getBasisFunctions(ksi, eta)); + + // Determine the number of nodes in the reference element based on the first element in the nop array + const nodesPerElement = nop[0].length; + + return { + residualVector, + jacobianMatrix, + localToGlobalMap, + basisFunctions, + gaussPoints, + gaussWeights, + integrationPoints2D, + referenceBasisFunctions, + nodesPerElement, + }; +} diff --git a/src/elements/isoparametricMapping.js b/src/elements/isoparametricMapping.js new file mode 100644 index 0000000..e8908f9 --- /dev/null +++ b/src/elements/isoparametricMapping.js @@ -0,0 +1,138 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ +/** + * Function to build the list of 2D integration points with their combined weights + * @param {array} gaussPoints - Gauss points from NumericalIntegration + * @param {array} gaussWeights - Gauss weights from NumericalIntegration + * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' + * @returns {array} Array of [ksi, eta, weight] entries + */ +export function getIntegrationPoints2D(gaussPoints, gaussWeights, elementShape = "quad") { + const integrationPoints = []; + + if (elementShape === "triangle") { + // Triangle rules already provide (ksi, eta) pairs with their weights + for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { + const [ksi, eta] = gaussPoints[gaussPointIndex]; + integrationPoints.push([ksi, eta, gaussWeights[gaussPointIndex]]); + } + } else { + // Quadrilaterals use the tensor product of the 1D Gauss rule + for (let gaussPointIndex1 = 0; gaussPointIndex1 < gaussPoints.length; gaussPointIndex1++) { + for (let gaussPointIndex2 = 0; gaussPointIndex2 < gaussPoints.length; gaussPointIndex2++) { + integrationPoints.push([ + gaussPoints[gaussPointIndex1], + gaussPoints[gaussPointIndex2], + gaussWeights[gaussPointIndex1] * gaussWeights[gaussPointIndex2], + ]); + } + } + } + + return integrationPoints; +} + +/** + * Function to perform isoparametric mapping for 1D elements + * @param {object} params - Parameters for the mapping. An optional `basisFunctionDerivX` array can be + * passed to be filled in place, to avoid allocating a new array at every integration point + * @returns {object} An object containing the mapped data + */ +export function performIsoparametricMapping1D(params) { + const { + basisFunction, + basisFunctionDerivKsi, + nodesXCoordinates, + localToGlobalMap, + nodesPerElement, + basisFunctionDerivX = [], + } = params; + + let xCoordinates = 0; + let ksiDerivX = 0; + + // Isoparametric mapping + for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { + xCoordinates += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; + ksiDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; + } + let detJacobian = ksiDerivX; + + // Compute x-derivative of basis functions + for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { + basisFunctionDerivX[localNodeIndex] = basisFunctionDerivKsi[localNodeIndex] / detJacobian; + } + + return { + xCoordinates, + detJacobian, + basisFunctionDerivX, + }; +} + +/** + * Function to perform isoparametric mapping for 2D elements + * @param {object} params - Parameters for the mapping. Optional `basisFunctionDerivX` and + * `basisFunctionDerivY` arrays can be passed to be filled in place, to avoid allocating new arrays at + * every integration point + * @returns {object} An object containing the mapped data + */ +export function performIsoparametricMapping2D(params) { + const { + basisFunction, + basisFunctionDerivKsi, + basisFunctionDerivEta, + nodesXCoordinates, + nodesYCoordinates, + localToGlobalMap, + nodesPerElement, + basisFunctionDerivX = [], + basisFunctionDerivY = [], + } = params; + + let xCoordinates = 0; + let yCoordinates = 0; + let ksiDerivX = 0; + let etaDerivX = 0; + let ksiDerivY = 0; + let etaDerivY = 0; + + // Isoparametric mapping + for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { + xCoordinates += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; + yCoordinates += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; + ksiDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; + etaDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex]; + ksiDerivY += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; + etaDerivY += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex]; + } + let detJacobian = ksiDerivX * etaDerivY - etaDerivX * ksiDerivY; + + // Compute x-derivative and y-derivative of basis functions + for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { + // The x-derivative of the n basis function + basisFunctionDerivX[localNodeIndex] = + (etaDerivY * basisFunctionDerivKsi[localNodeIndex] - + ksiDerivY * basisFunctionDerivEta[localNodeIndex]) / + detJacobian; + // The y-derivative of the n basis function + basisFunctionDerivY[localNodeIndex] = + (ksiDerivX * basisFunctionDerivEta[localNodeIndex] - + etaDerivX * basisFunctionDerivKsi[localNodeIndex]) / + detJacobian; + } + + return { + xCoordinates, + yCoordinates, + detJacobian, + basisFunctionDerivX, + basisFunctionDerivY, + }; +} diff --git a/src/numerics/numericalIntegration.js b/src/elements/numericalIntegration.js similarity index 100% rename from src/numerics/numericalIntegration.js rename to src/elements/numericalIntegration.js diff --git a/src/index.js b/src/index.js index cb07256..0ac75b8 100644 --- a/src/index.js +++ b/src/index.js @@ -8,7 +8,7 @@ */ export { FEAScriptModel } from "./FEAScript.js"; -export { importGmshMesh } from "./readers/gmshReader.js"; +export { importGmshMesh } from "./io/gmshReader.js"; export { logSystem } from "./utilities/logging.js"; export { plotSolution, plotInterpolatedSolution } from "./visualization/plotlyPlot.js"; export { diff --git a/src/io/gmshReader.js b/src/io/gmshReader.js new file mode 100644 index 0000000..8486eec --- /dev/null +++ b/src/io/gmshReader.js @@ -0,0 +1,695 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { triangleBoundarySides } from "../elements/basisFunctions.js"; +import { debugLog, errorLog } from "../utilities/logging.js"; + +/** + * Function to import mesh data from Gmsh (.msh v4.1) containing quadrilateral and triangular elements + * @param {File} file + * @returns {object} + */ +const importGmshMesh = async (file) => { + let result = { + nodesXCoordinates: [], + nodesYCoordinates: [], + nodalNumbering: { + quadElements: [], + triangleElements: [], + }, + boundaryElements: [], + boundaryConditions: [], + boundaryNodePairs: {}, + gmshV: 0, + ascii: false, + fltBytes: "8", + totalNodesX: 0, + totalNodesY: 0, + physicalPropMap: [], + elementTypes: {}, + }; + + // Entities to physical tags map + const entityPhysicalMap = { curves: {} }; + + let content = await file.text(); + let lines = content + .split("\n") + .map((l) => l.trim()) + .filter((l) => l !== ""); + + let section = ""; + let lineIndex = 0; + + // Nodes + let nodeEntityBlocks = 0; + let totalNodes = 0; + let nodeBlocksProcessed = 0; + let currentNodeBlock = { numNodes: 0 }; + let nodeTags = []; + 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; + let currentElementBlock = { numElements: 0 }; + let elementsProcessedInBlock = 0; + + let boundaryElementsByTag = {}; + + // Entities helpers + let entityCounts = null; + let entitiesPhase = null; + let processedPoints = 0; + let processedCurves = 0; + let processedSurfaces = 0; + + while (lineIndex < lines.length) { + const line = lines[lineIndex]; + + // Section switches + if (line === "$MeshFormat") { + section = "meshFormat"; + lineIndex++; + continue; + } + if (line === "$EndMeshFormat") { + section = ""; + lineIndex++; + continue; + } + if (line === "$PhysicalNames") { + section = "physicalNames"; + lineIndex++; + continue; + } + if (line === "$EndPhysicalNames") { + section = ""; + lineIndex++; + continue; + } + if (line === "$Entities") { + section = "entities"; + entitiesPhase = "counts"; + lineIndex++; + continue; + } + if (line === "$EndEntities") { + section = ""; + entityCounts = null; + entitiesPhase = null; + lineIndex++; + continue; + } + if (line === "$Nodes") { + section = "nodes"; + lineIndex++; + continue; + } + if (line === "$EndNodes") { + section = ""; + lineIndex++; + continue; + } + if (line === "$Elements") { + section = "elements"; + lineIndex++; + continue; + } + if (line === "$EndElements") { + section = ""; + lineIndex++; + continue; + } + + const parts = line.split(/\s+/); + + // Mesh format + if (section === "meshFormat") { + result.gmshV = parseFloat(parts[0]); + result.ascii = parts[1] === "0"; + result.fltBytes = parts[2]; + } + + // Physical names + else if (section === "physicalNames") { + const dimension = parseInt(parts[0], 10); + const tag = parseInt(parts[1], 10); + let name = parts.slice(2).join(" ").replace(/^"|"$/g, ""); + + result.physicalPropMap.push({ tag, dimension, name }); + } + + // Entities + else if (section === "entities") { + if (entitiesPhase === "counts") { + entityCounts = { + points: parseInt(parts[0], 10), + curves: parseInt(parts[1], 10), + surfaces: parseInt(parts[2], 10), + volumes: parseInt(parts[3], 10), + }; + entitiesPhase = "points"; + } else if (entitiesPhase === "points") { + processedPoints++; + if (processedPoints === entityCounts.points) entitiesPhase = "curves"; + } else if (entitiesPhase === "curves") { + const tag = parseInt(parts[0], 10); + const numPhysical = parseInt(parts[7], 10); + const physTags = parts.slice(8, 8 + numPhysical).map((p) => parseInt(p, 10)); + entityPhysicalMap.curves[tag] = physTags; + processedCurves++; + if (processedCurves === entityCounts.curves) entitiesPhase = "surfaces"; + } else if (entitiesPhase === "surfaces") { + processedSurfaces++; + if (processedSurfaces === entityCounts.surfaces) entitiesPhase = "volumes"; + } + } + + // Nodes + else if (section === "nodes") { + 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) { + currentNodeBlock = { + dim: parseInt(parts[0], 10), + tag: parseInt(parts[1], 10), + parametric: parseInt(parts[2], 10), + numNodes: parseInt(parts[3], 10), + }; + nodeTags = []; + nodeTagsCollected = 0; + nodeCoordinatesCollected = 0; + } else if (nodeTagsCollected < currentNodeBlock.numNodes) { + for (let p of parts) { + nodeTags.push(parseInt(p, 10)); + nodeTagsCollected++; + if (nodeTagsCollected === currentNodeBlock.numNodes) break; + } + } else if (nodeCoordinatesCollected < currentNodeBlock.numNodes) { + 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++; + if (nodeCoordinatesCollected === currentNodeBlock.numNodes) { + nodeBlocksProcessed++; + currentNodeBlock = { numNodes: 0 }; + } + } + } + + // Elements + else if (section === "elements") { + if (elementEntityBlocks === 0) { + elementEntityBlocks = parseInt(parts[0], 10); + } else if (elementBlocksProcessed < elementEntityBlocks && currentElementBlock.numElements === 0) { + currentElementBlock = { + dim: parseInt(parts[0], 10), + tag: parseInt(parts[1], 10), // entity tag + elementType: parseInt(parts[2], 10), + numElements: parseInt(parts[3], 10), + }; + + result.elementTypes[currentElementBlock.elementType] = + (result.elementTypes[currentElementBlock.elementType] || 0) + currentElementBlock.numElements; + + elementsProcessedInBlock = 0; + } else if (elementsProcessedInBlock < currentElementBlock.numElements) { + 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; + if (currentElementBlock.dim === 1) { + const mapped = entityPhysicalMap.curves[currentElementBlock.tag]; + if (mapped && mapped.length > 0) physicalTag = mapped[0]; + } + + if (currentElementBlock.elementType === 1 || currentElementBlock.elementType === 8) { + if (!boundaryElementsByTag[physicalTag]) boundaryElementsByTag[physicalTag] = []; + boundaryElementsByTag[physicalTag].push(nodeIndices); + + if (!result.boundaryNodePairs[physicalTag]) result.boundaryNodePairs[physicalTag] = []; + result.boundaryNodePairs[physicalTag].push(nodeIndices); + } else if (currentElementBlock.elementType === 2 || currentElementBlock.elementType === 9) { + result.nodalNumbering.triangleElements.push(nodeIndices); + } else if (currentElementBlock.elementType === 3 || currentElementBlock.elementType === 10) { + result.nodalNumbering.quadElements.push(nodeIndices); + } + + elementsProcessedInBlock++; + if (elementsProcessedInBlock === currentElementBlock.numElements) { + elementBlocksProcessed++; + currentElementBlock = { numElements: 0 }; + } + } + } + + lineIndex++; + } + + // Boundary conditions + result.physicalPropMap.forEach((prop) => { + if (prop.dimension === 1) { + const nodes = boundaryElementsByTag[prop.tag] || []; + if (nodes.length > 0) { + result.boundaryConditions.push({ + name: prop.name, + tag: prop.tag, + nodes, + }); + } + } + }); + + debugLog(`Parsed boundary node pairs by physical tag: ${JSON.stringify(result.boundaryNodePairs)}`); + + return result; +}; + +export { importGmshMesh }; + +/** + * Function to convert a mesh imported with importGmshMesh to the FEAScript format. The nodal numbering is + * mapped in place to the FEAScript node ordering, and boundary node pairs are converted to element sides + * @param {object} parsedMesh - Mesh data returned by importGmshMesh (modified in place) + * @param {string} elementOrder - The order of elements, either 'linear' or 'quadratic' + * @returns {boolean} True when the boundary elements have been processed + */ +export function mapGmshMeshToFEAScript(parsedMesh, elementOrder) { + let boundaryElementsProcessed = false; + + if (!parsedMesh.nodalNumbering) { + errorLog("No valid nodal numbering found in the parsed mesh."); + } + + // If this parsed mesh was already converted in a previous run, don't re-process it. + // Just mark this Mesh instance as ready so prepareMesh() doesn't fall back to generateMesh(). + if (Array.isArray(parsedMesh.nodalNumbering)) { + parsedMesh.boundaryElementsProcessed = true; + return true; + } + + if (typeof parsedMesh.nodalNumbering === "object" && !Array.isArray(parsedMesh.nodalNumbering)) { + // Store the nodal numbering structure before converting + const quadElements = parsedMesh.nodalNumbering.quadElements || []; + const triangleElements = parsedMesh.nodalNumbering.triangleElements || []; + + debugLog( + "Initial parsed mesh nodal numbering from Gmsh format: " + JSON.stringify(parsedMesh.nodalNumbering), + ); + + const hasQuadElements = quadElements.length > 0; + const hasTriangleElements = triangleElements.length > 0; + + if (hasQuadElements && hasTriangleElements) { + errorLog("Mixed quadrilateral and triangular meshes are not supported yet."); + } + + // Check if it has quadElements or triangleElements structure from gmshReader + if (hasQuadElements) { + parsedMesh.elementShape = "quad"; + // Map nodal numbering from Gmsh format to FEAScript format for quad elements + const mappedNodalNumbering = []; + + for (let elementIndex = 0; elementIndex < quadElements.length; elementIndex++) { + const gmshNodes = quadElements[elementIndex]; + const FEAScriptNodes = new Array(gmshNodes.length); + + // Check for element type based on number of nodes + if (gmshNodes.length === 4) { + // Simple mapping for linear quad elements (4 nodes) + // Gmsh: FEAScript: + // 3 --- 2 1 --- 3 + // | | --> | | + // 0 --- 1 0 --- 2 + + FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 + FEAScriptNodes[1] = gmshNodes[3]; // 3 -> 1 + FEAScriptNodes[2] = gmshNodes[1]; // 1 -> 2 + FEAScriptNodes[3] = gmshNodes[2]; // 2 -> 3 + } else if (gmshNodes.length === 9) { + // Mapping for quadratic quad elements (9 nodes) + // Gmsh: FEAScript: + // 3--6--2 2--5--8 + // | | | | + // 7 8 5 --> 1 4 7 + // | | | | + // 0--4--1 0--3--6 + + FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 + FEAScriptNodes[1] = gmshNodes[7]; // 7 -> 1 + FEAScriptNodes[2] = gmshNodes[3]; // 3 -> 2 + FEAScriptNodes[3] = gmshNodes[4]; // 4 -> 3 + FEAScriptNodes[4] = gmshNodes[8]; // 8 -> 4 + FEAScriptNodes[5] = gmshNodes[6]; // 6 -> 5 + FEAScriptNodes[6] = gmshNodes[1]; // 1 -> 6 + FEAScriptNodes[7] = gmshNodes[5]; // 5 -> 7 + FEAScriptNodes[8] = gmshNodes[2]; // 2 -> 8 + } + + mappedNodalNumbering.push(FEAScriptNodes); + } + + parsedMesh.nodalNumbering = mappedNodalNumbering; + } else if (hasTriangleElements) { + parsedMesh.elementShape = "triangle"; + + // Map nodal numbering from Gmsh format to FEAScript format for triangular elements + const mappedNodalNumbering = []; + + for (let elementIndex = 0; elementIndex < triangleElements.length; elementIndex++) { + const gmshNodes = triangleElements[elementIndex]; + const FEAScriptNodes = new Array(gmshNodes.length); + + // Check for element type based on number of nodes + if (gmshNodes.length === 3) { + // Mapping for linear triangular elements (3 nodes) + // Gmsh: FEAScript: + // 2 1 + // | \ | \ + // | \ --> | \ + // 0 --- 1 0 --- 2 + + FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 + FEAScriptNodes[1] = gmshNodes[2]; // 2 -> 1 + FEAScriptNodes[2] = gmshNodes[1]; // 1 -> 2 + } else if (gmshNodes.length === 6) { + // Mapping for quadratic triangular elements (6 nodes) + // Gmsh: FEAScript: + // 2 2 + // | \ | \ + // 5 4 --> 1 4 + // | \ | \ + // 0--3--1 0--3--5 + + FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 + FEAScriptNodes[1] = gmshNodes[5]; // 5 -> 1 + FEAScriptNodes[2] = gmshNodes[2]; // 2 -> 2 + FEAScriptNodes[3] = gmshNodes[3]; // 3 -> 3 + FEAScriptNodes[4] = gmshNodes[4]; // 4 -> 4 + FEAScriptNodes[5] = gmshNodes[1]; // 1 -> 5 + } + + mappedNodalNumbering.push(FEAScriptNodes); + } + + parsedMesh.nodalNumbering = mappedNodalNumbering; + + const expectedNodesPerTriangle = elementOrder === "quadratic" ? 6 : 3; + if (triangleElements[0].length !== expectedNodesPerTriangle) { + errorLog( + `The Gmsh mesh contains ${triangleElements[0].length}-node triangles, but elementOrder is '${elementOrder}' (expects ${expectedNodesPerTriangle}-node triangles). Use 'linear' for first-order and 'quadratic' for second-order meshes.`, + ); + } + } else { + errorLog("Element type is neither triangle nor quad; mapping for this type is not implemented yet."); + } + + debugLog( + "Nodal numbering after mapping from Gmsh to FEAScript format: " + + JSON.stringify(parsedMesh.nodalNumbering), + ); + + // Process boundary elements if they exist and if physical property mapping exists + if (parsedMesh.physicalPropMap && parsedMesh.boundaryElements) { + // Check if boundary elements need to be processed + if ( + Array.isArray(parsedMesh.boundaryElements) && + parsedMesh.boundaryElements.length > 0 && + parsedMesh.boundaryElements[0] === undefined + ) { + // Create a new array without the empty first element + const fixedBoundaryElements = []; + for (let i = 1; i < parsedMesh.boundaryElements.length; i++) { + if (parsedMesh.boundaryElements[i]) { + fixedBoundaryElements.push(parsedMesh.boundaryElements[i]); + } + } + parsedMesh.boundaryElements = fixedBoundaryElements; + } + + // If boundary node pairs exist but boundary elements haven't been processed + if (parsedMesh.boundaryNodePairs && !parsedMesh.boundaryElementsProcessed) { + // Reset boundary elements array + parsedMesh.boundaryElements = []; + + // Process each physical property from the Gmsh file + parsedMesh.physicalPropMap.forEach((prop) => { + // Only process 1D physical entities (boundary lines) + if (prop.dimension === 1) { + // Get all node pairs for this boundary + const boundaryNodePairs = parsedMesh.boundaryNodePairs[prop.tag] || []; + + if (boundaryNodePairs.length > 0) { + // Initialize array for this boundary tag + if (!parsedMesh.boundaryElements[prop.tag]) { + parsedMesh.boundaryElements[prop.tag] = []; + } + + // For each boundary line segment (defined by a pair of nodes) + boundaryNodePairs.forEach((nodesPair) => { + const node1 = nodesPair[0]; // First node in the pair + const node2 = nodesPair[1]; // Second node in the pair + + debugLog( + `Processing boundary node pair: [${node1}, ${node2}] for boundary ${prop.tag} (${ + prop.name || "unnamed" + })`, + ); + + // Search through all elements to find which one contains both nodes + let foundElement = false; + + // Loop through all elements in the mesh + for (let elementIndex = 0; elementIndex < parsedMesh.nodalNumbering.length; elementIndex++) { + const elementConnectivity = parsedMesh.nodalNumbering[elementIndex]; + + // For linear (3 nodes) and quadratic (6 nodes) triangular elements + if (parsedMesh.elementShape === "triangle") { + // Based on FEAScript triangle numbering: + // Linear: Quadratic: + // 1 2 + // | \ | \ + // | \ 1 4 + // | \ | \ + // 0 --- 2 0--3--5 + // Side 0: bottom, side 1: left, side 2: hypotenuse + const sides = + triangleBoundarySides[elementConnectivity.length === 6 ? "quadratic" : "linear"]; + + // Boundary node pairs are the end vertices of the boundary line, so they are + // compared against the first and last node of each side + const side = Object.keys(sides).find((sideKey) => { + const sideNodes = sides[sideKey]; + const firstNode = elementConnectivity[sideNodes[0]]; + const lastNode = elementConnectivity[sideNodes[sideNodes.length - 1]]; + return ( + (firstNode === node1 && lastNode === node2) || + (firstNode === node2 && lastNode === node1) + ); + }); + + if (side !== undefined) { + debugLog( + ` Found triangle ${elementIndex} containing boundary nodes on side ${side}. Element nodes: [${elementConnectivity.join( + ", ", + )}]`, + ); + + parsedMesh.boundaryElements[prop.tag].push([elementIndex, Number(side)]); + foundElement = true; + break; + } + } + // For linear quadrilateral linear elements (4 nodes) + else if (elementConnectivity.length === 4) { + // Check if both boundary nodes are in this element + if (elementConnectivity.includes(node1) && elementConnectivity.includes(node2)) { + // Find which side of the element these nodes form + let side; + + const node1Index = elementConnectivity.indexOf(node1); + const node2Index = elementConnectivity.indexOf(node2); + + debugLog( + ` Found element ${elementIndex} containing boundary nodes. Element nodes: [${elementConnectivity.join( + ", ", + )}]`, + ); + debugLog( + ` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`, + ); + + // Based on FEAScript linear quadrilateral numbering: + // 1 --- 3 + // | | + // 0 --- 2 + + if ((node1Index === 0 && node2Index === 2) || (node1Index === 2 && node2Index === 0)) { + side = 0; // Bottom side + debugLog(` These nodes form the BOTTOM side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 0 && node2Index === 1) || + (node1Index === 1 && node2Index === 0) + ) { + side = 1; // Left side + debugLog(` These nodes form the LEFT side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 1 && node2Index === 3) || + (node1Index === 3 && node2Index === 1) + ) { + side = 2; // Top side + debugLog(` These nodes form the TOP side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 2 && node2Index === 3) || + (node1Index === 3 && node2Index === 2) + ) { + side = 3; // Right side + debugLog(` These nodes form the RIGHT side (${side}) of element ${elementIndex}`); + } + + // Add the element and side to the boundary elements array + parsedMesh.boundaryElements[prop.tag].push([elementIndex, side]); + debugLog( + ` Added element-side pair [${elementIndex}, ${side}] to boundary tag ${prop.tag}`, + ); + foundElement = true; + break; + } + } else if (elementConnectivity.length === 9) { + // For quadratic quadrilateral elements (9 nodes) + // Check if both boundary nodes are in this element + if (elementConnectivity.includes(node1) && elementConnectivity.includes(node2)) { + // Find which side of the element these nodes form + let side; + + const node1Index = elementConnectivity.indexOf(node1); + const node2Index = elementConnectivity.indexOf(node2); + + debugLog( + ` Found element ${elementIndex} containing boundary nodes. Element nodes: [${elementConnectivity.join( + ", ", + )}]`, + ); + debugLog( + ` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`, + ); + + // Based on FEAScript quadratic quadrilateral numbering: + // 2--5--8 + // | | + // 1 4 7 + // | | + // 0--3--6 + + // TODO: Transform into dictionaries for better readability + if ( + (node1Index === 0 && node2Index === 6) || + (node1Index === 6 && node2Index === 0) || + (node1Index === 0 && node2Index === 3) || + (node1Index === 3 && node2Index === 0) || + (node1Index === 3 && node2Index === 6) || + (node1Index === 6 && node2Index === 3) + ) { + side = 0; // Bottom side (nodes 0, 3, 6) + debugLog(` These nodes form the BOTTOM side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 0 && node2Index === 2) || + (node1Index === 2 && node2Index === 0) || + (node1Index === 0 && node2Index === 1) || + (node1Index === 1 && node2Index === 0) || + (node1Index === 1 && node2Index === 2) || + (node1Index === 2 && node2Index === 1) + ) { + side = 1; // Left side (nodes 0, 1, 2) + debugLog(` These nodes form the LEFT side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 2 && node2Index === 8) || + (node1Index === 8 && node2Index === 2) || + (node1Index === 2 && node2Index === 5) || + (node1Index === 5 && node2Index === 2) || + (node1Index === 5 && node2Index === 8) || + (node1Index === 8 && node2Index === 5) + ) { + side = 2; // Top side (nodes 2, 5, 8) + debugLog(` These nodes form the TOP side (${side}) of element ${elementIndex}`); + } else if ( + (node1Index === 6 && node2Index === 8) || + (node1Index === 8 && node2Index === 6) || + (node1Index === 6 && node2Index === 7) || + (node1Index === 7 && node2Index === 6) || + (node1Index === 7 && node2Index === 8) || + (node1Index === 8 && node2Index === 7) + ) { + side = 3; // Right side (nodes 6, 7, 8) + debugLog(` These nodes form the RIGHT side (${side}) of element ${elementIndex}`); + } + + // Add the element and side to the boundary elements array + parsedMesh.boundaryElements[prop.tag].push([elementIndex, side]); + debugLog( + ` Added element-side pair [${elementIndex}, ${side}] to boundary tag ${prop.tag}`, + ); + foundElement = true; + break; + } + } + } + + if (!foundElement) { + errorLog( + `Could not find element containing boundary nodes ${node1} and ${node2}. Boundary may be incomplete.`, + ); + } + }); + } + } + }); + + // Mark as processed + boundaryElementsProcessed = true; + + // Fix boundary elements array - remove undefined entries + if (parsedMesh.boundaryElements.length > 0 && parsedMesh.boundaryElements[0] === undefined) { + const fixedBoundaryElements = []; + for (let i = 1; i < parsedMesh.boundaryElements.length; i++) { + if (parsedMesh.boundaryElements[i]) { + fixedBoundaryElements.push(parsedMesh.boundaryElements[i]); + } + } + parsedMesh.boundaryElements = fixedBoundaryElements; + } + } + } + } + + return boundaryElementsProcessed; +} diff --git a/src/numerics/blasUtilities.js b/src/linalg/blasUtilities.js similarity index 100% rename from src/numerics/blasUtilities.js rename to src/linalg/blasUtilities.js diff --git a/src/numerics/index.js b/src/linalg/index.js similarity index 100% rename from src/numerics/index.js rename to src/linalg/index.js diff --git a/src/numerics/sparseMatrix.js b/src/linalg/sparseMatrix.js similarity index 100% rename from src/numerics/sparseMatrix.js rename to src/linalg/sparseMatrix.js diff --git a/src/mesh/meshGeneration.js b/src/mesh/meshGeneration.js index d889325..d4a4d60 100644 --- a/src/mesh/meshGeneration.js +++ b/src/mesh/meshGeneration.js @@ -8,7 +8,7 @@ */ // Internal imports -import { triangleBoundarySides } from "./basisFunctions.js"; +import { mapGmshMeshToFEAScript } from "../io/gmshReader.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; /** @@ -55,406 +55,7 @@ export class Mesh { * Method to parse the mesh from the Gmsh format to the FEAScript format */ parseMeshFromGmsh() { - if (!this.parsedMesh.nodalNumbering) { - errorLog("No valid nodal numbering found in the parsed mesh."); - } - - // If this parsed mesh was already converted in a previous run, don't re-process it. - // Just mark this Mesh instance as ready so prepareMesh() doesn't fall back to generateMesh(). - if (Array.isArray(this.parsedMesh.nodalNumbering)) { - this.boundaryElementsProcessed = true; - this.parsedMesh.boundaryElementsProcessed = true; - return this.parsedMesh; - } - - if ( - typeof this.parsedMesh.nodalNumbering === "object" && - !Array.isArray(this.parsedMesh.nodalNumbering) - ) { - // Store the nodal numbering structure before converting - const quadElements = this.parsedMesh.nodalNumbering.quadElements || []; - const triangleElements = this.parsedMesh.nodalNumbering.triangleElements || []; - - debugLog( - "Initial parsed mesh nodal numbering from Gmsh format: " + - JSON.stringify(this.parsedMesh.nodalNumbering), - ); - - const hasQuadElements = quadElements.length > 0; - const hasTriangleElements = triangleElements.length > 0; - - if (hasQuadElements && hasTriangleElements) { - errorLog("Mixed quadrilateral and triangular meshes are not supported yet."); - } - - // Check if it has quadElements or triangleElements structure from gmshReader - if (hasQuadElements) { - this.parsedMesh.elementShape = "quad"; - // Map nodal numbering from Gmsh format to FEAScript format for quad elements - const mappedNodalNumbering = []; - - for (let elementIndex = 0; elementIndex < quadElements.length; elementIndex++) { - const gmshNodes = quadElements[elementIndex]; - const FEAScriptNodes = new Array(gmshNodes.length); - - // Check for element type based on number of nodes - if (gmshNodes.length === 4) { - // Simple mapping for linear quad elements (4 nodes) - // Gmsh: FEAScript: - // 3 --- 2 1 --- 3 - // | | --> | | - // 0 --- 1 0 --- 2 - - FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 - FEAScriptNodes[1] = gmshNodes[3]; // 3 -> 1 - FEAScriptNodes[2] = gmshNodes[1]; // 1 -> 2 - FEAScriptNodes[3] = gmshNodes[2]; // 2 -> 3 - } else if (gmshNodes.length === 9) { - // Mapping for quadratic quad elements (9 nodes) - // Gmsh: FEAScript: - // 3--6--2 2--5--8 - // | | | | - // 7 8 5 --> 1 4 7 - // | | | | - // 0--4--1 0--3--6 - - FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 - FEAScriptNodes[1] = gmshNodes[7]; // 7 -> 1 - FEAScriptNodes[2] = gmshNodes[3]; // 3 -> 2 - FEAScriptNodes[3] = gmshNodes[4]; // 4 -> 3 - FEAScriptNodes[4] = gmshNodes[8]; // 8 -> 4 - FEAScriptNodes[5] = gmshNodes[6]; // 6 -> 5 - FEAScriptNodes[6] = gmshNodes[1]; // 1 -> 6 - FEAScriptNodes[7] = gmshNodes[5]; // 5 -> 7 - FEAScriptNodes[8] = gmshNodes[2]; // 2 -> 8 - } - - mappedNodalNumbering.push(FEAScriptNodes); - } - - this.parsedMesh.nodalNumbering = mappedNodalNumbering; - } else if (hasTriangleElements) { - this.parsedMesh.elementShape = "triangle"; - - // Map nodal numbering from Gmsh format to FEAScript format for triangular elements - const mappedNodalNumbering = []; - - for (let elementIndex = 0; elementIndex < triangleElements.length; elementIndex++) { - const gmshNodes = triangleElements[elementIndex]; - const FEAScriptNodes = new Array(gmshNodes.length); - - // Check for element type based on number of nodes - if (gmshNodes.length === 3) { - // Mapping for linear triangular elements (3 nodes) - // Gmsh: FEAScript: - // 2 1 - // | \ | \ - // | \ --> | \ - // 0 --- 1 0 --- 2 - - FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 - FEAScriptNodes[1] = gmshNodes[2]; // 2 -> 1 - FEAScriptNodes[2] = gmshNodes[1]; // 1 -> 2 - } else if (gmshNodes.length === 6) { - // Mapping for quadratic triangular elements (6 nodes) - // Gmsh: FEAScript: - // 2 2 - // | \ | \ - // 5 4 --> 1 4 - // | \ | \ - // 0--3--1 0--3--5 - - FEAScriptNodes[0] = gmshNodes[0]; // 0 -> 0 - FEAScriptNodes[1] = gmshNodes[5]; // 5 -> 1 - FEAScriptNodes[2] = gmshNodes[2]; // 2 -> 2 - FEAScriptNodes[3] = gmshNodes[3]; // 3 -> 3 - FEAScriptNodes[4] = gmshNodes[4]; // 4 -> 4 - FEAScriptNodes[5] = gmshNodes[1]; // 1 -> 5 - } - - mappedNodalNumbering.push(FEAScriptNodes); - } - - this.parsedMesh.nodalNumbering = mappedNodalNumbering; - - const expectedNodesPerTriangle = this.elementOrder === "quadratic" ? 6 : 3; - if (triangleElements[0].length !== expectedNodesPerTriangle) { - errorLog( - `The Gmsh mesh contains ${triangleElements[0].length}-node triangles, but elementOrder is '${this.elementOrder}' (expects ${expectedNodesPerTriangle}-node triangles). Use 'linear' for first-order and 'quadratic' for second-order meshes.`, - ); - } - } else { - errorLog("Element type is neither triangle nor quad; mapping for this type is not implemented yet."); - } - - debugLog( - "Nodal numbering after mapping from Gmsh to FEAScript format: " + - JSON.stringify(this.parsedMesh.nodalNumbering), - ); - - // Process boundary elements if they exist and if physical property mapping exists - if (this.parsedMesh.physicalPropMap && this.parsedMesh.boundaryElements) { - // Check if boundary elements need to be processed - if ( - Array.isArray(this.parsedMesh.boundaryElements) && - this.parsedMesh.boundaryElements.length > 0 && - this.parsedMesh.boundaryElements[0] === undefined - ) { - // Create a new array without the empty first element - const fixedBoundaryElements = []; - for (let i = 1; i < this.parsedMesh.boundaryElements.length; i++) { - if (this.parsedMesh.boundaryElements[i]) { - fixedBoundaryElements.push(this.parsedMesh.boundaryElements[i]); - } - } - this.parsedMesh.boundaryElements = fixedBoundaryElements; - } - - // If boundary node pairs exist but boundary elements haven't been processed - if (this.parsedMesh.boundaryNodePairs && !this.parsedMesh.boundaryElementsProcessed) { - // Reset boundary elements array - this.parsedMesh.boundaryElements = []; - - // Process each physical property from the Gmsh file - this.parsedMesh.physicalPropMap.forEach((prop) => { - // Only process 1D physical entities (boundary lines) - if (prop.dimension === 1) { - // Get all node pairs for this boundary - const boundaryNodePairs = this.parsedMesh.boundaryNodePairs[prop.tag] || []; - - if (boundaryNodePairs.length > 0) { - // Initialize array for this boundary tag - if (!this.parsedMesh.boundaryElements[prop.tag]) { - this.parsedMesh.boundaryElements[prop.tag] = []; - } - - // For each boundary line segment (defined by a pair of nodes) - boundaryNodePairs.forEach((nodesPair) => { - const node1 = nodesPair[0]; // First node in the pair - const node2 = nodesPair[1]; // Second node in the pair - - debugLog( - `Processing boundary node pair: [${node1}, ${node2}] for boundary ${prop.tag} (${ - prop.name || "unnamed" - })`, - ); - - // Search through all elements to find which one contains both nodes - let foundElement = false; - - // Loop through all elements in the mesh - for ( - let elementIndex = 0; - elementIndex < this.parsedMesh.nodalNumbering.length; - elementIndex++ - ) { - const elementConnectivity = this.parsedMesh.nodalNumbering[elementIndex]; - - // For linear (3 nodes) and quadratic (6 nodes) triangular elements - if (this.parsedMesh.elementShape === "triangle") { - // Based on FEAScript triangle numbering: - // Linear: Quadratic: - // 1 2 - // | \ | \ - // | \ 1 4 - // | \ | \ - // 0 --- 2 0--3--5 - // Side 0: bottom, side 1: left, side 2: hypotenuse - const sides = - triangleBoundarySides[elementConnectivity.length === 6 ? "quadratic" : "linear"]; - - // Boundary node pairs are the end vertices of the boundary line, so they are - // compared against the first and last node of each side - const side = Object.keys(sides).find((sideKey) => { - const sideNodes = sides[sideKey]; - const firstNode = elementConnectivity[sideNodes[0]]; - const lastNode = elementConnectivity[sideNodes[sideNodes.length - 1]]; - return ( - (firstNode === node1 && lastNode === node2) || - (firstNode === node2 && lastNode === node1) - ); - }); - - if (side !== undefined) { - debugLog( - ` Found triangle ${elementIndex} containing boundary nodes on side ${side}. Element nodes: [${elementConnectivity.join( - ", ", - )}]`, - ); - - this.parsedMesh.boundaryElements[prop.tag].push([elementIndex, Number(side)]); - foundElement = true; - break; - } - } - // For linear quadrilateral linear elements (4 nodes) - else if (elementConnectivity.length === 4) { - // Check if both boundary nodes are in this element - if (elementConnectivity.includes(node1) && elementConnectivity.includes(node2)) { - // Find which side of the element these nodes form - let side; - - const node1Index = elementConnectivity.indexOf(node1); - const node2Index = elementConnectivity.indexOf(node2); - - debugLog( - ` Found element ${elementIndex} containing boundary nodes. Element nodes: [${elementConnectivity.join( - ", ", - )}]`, - ); - debugLog( - ` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`, - ); - - // Based on FEAScript linear quadrilateral numbering: - // 1 --- 3 - // | | - // 0 --- 2 - - if ( - (node1Index === 0 && node2Index === 2) || - (node1Index === 2 && node2Index === 0) - ) { - side = 0; // Bottom side - debugLog(` These nodes form the BOTTOM side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 0 && node2Index === 1) || - (node1Index === 1 && node2Index === 0) - ) { - side = 1; // Left side - debugLog(` These nodes form the LEFT side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 1 && node2Index === 3) || - (node1Index === 3 && node2Index === 1) - ) { - side = 2; // Top side - debugLog(` These nodes form the TOP side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 2 && node2Index === 3) || - (node1Index === 3 && node2Index === 2) - ) { - side = 3; // Right side - debugLog(` These nodes form the RIGHT side (${side}) of element ${elementIndex}`); - } - - // Add the element and side to the boundary elements array - this.parsedMesh.boundaryElements[prop.tag].push([elementIndex, side]); - debugLog( - ` Added element-side pair [${elementIndex}, ${side}] to boundary tag ${prop.tag}`, - ); - foundElement = true; - break; - } - } else if (elementConnectivity.length === 9) { - // For quadratic quadrilateral elements (9 nodes) - // Check if both boundary nodes are in this element - if (elementConnectivity.includes(node1) && elementConnectivity.includes(node2)) { - // Find which side of the element these nodes form - let side; - - const node1Index = elementConnectivity.indexOf(node1); - const node2Index = elementConnectivity.indexOf(node2); - - debugLog( - ` Found element ${elementIndex} containing boundary nodes. Element nodes: [${elementConnectivity.join( - ", ", - )}]`, - ); - debugLog( - ` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`, - ); - - // Based on FEAScript quadratic quadrilateral numbering: - // 2--5--8 - // | | - // 1 4 7 - // | | - // 0--3--6 - - // TODO: Transform into dictionaries for better readability - if ( - (node1Index === 0 && node2Index === 6) || - (node1Index === 6 && node2Index === 0) || - (node1Index === 0 && node2Index === 3) || - (node1Index === 3 && node2Index === 0) || - (node1Index === 3 && node2Index === 6) || - (node1Index === 6 && node2Index === 3) - ) { - side = 0; // Bottom side (nodes 0, 3, 6) - debugLog(` These nodes form the BOTTOM side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 0 && node2Index === 2) || - (node1Index === 2 && node2Index === 0) || - (node1Index === 0 && node2Index === 1) || - (node1Index === 1 && node2Index === 0) || - (node1Index === 1 && node2Index === 2) || - (node1Index === 2 && node2Index === 1) - ) { - side = 1; // Left side (nodes 0, 1, 2) - debugLog(` These nodes form the LEFT side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 2 && node2Index === 8) || - (node1Index === 8 && node2Index === 2) || - (node1Index === 2 && node2Index === 5) || - (node1Index === 5 && node2Index === 2) || - (node1Index === 5 && node2Index === 8) || - (node1Index === 8 && node2Index === 5) - ) { - side = 2; // Top side (nodes 2, 5, 8) - debugLog(` These nodes form the TOP side (${side}) of element ${elementIndex}`); - } else if ( - (node1Index === 6 && node2Index === 8) || - (node1Index === 8 && node2Index === 6) || - (node1Index === 6 && node2Index === 7) || - (node1Index === 7 && node2Index === 6) || - (node1Index === 7 && node2Index === 8) || - (node1Index === 8 && node2Index === 7) - ) { - side = 3; // Right side (nodes 6, 7, 8) - debugLog(` These nodes form the RIGHT side (${side}) of element ${elementIndex}`); - } - - // Add the element and side to the boundary elements array - this.parsedMesh.boundaryElements[prop.tag].push([elementIndex, side]); - debugLog( - ` Added element-side pair [${elementIndex}, ${side}] to boundary tag ${prop.tag}`, - ); - foundElement = true; - break; - } - } - } - - if (!foundElement) { - errorLog( - `Could not find element containing boundary nodes ${node1} and ${node2}. Boundary may be incomplete.`, - ); - } - }); - } - } - }); - - // Mark as processed - this.boundaryElementsProcessed = true; - - // Fix boundary elements array - remove undefined entries - if ( - this.parsedMesh.boundaryElements.length > 0 && - this.parsedMesh.boundaryElements[0] === undefined - ) { - const fixedBoundaryElements = []; - for (let i = 1; i < this.parsedMesh.boundaryElements.length; i++) { - if (this.parsedMesh.boundaryElements[i]) { - fixedBoundaryElements.push(this.parsedMesh.boundaryElements[i]); - } - } - this.parsedMesh.boundaryElements = fixedBoundaryElements; - } - } - } - } - + this.boundaryElementsProcessed = mapGmshMeshToFEAScript(this.parsedMesh, this.elementOrder); return this.parsedMesh; } } diff --git a/src/mesh/meshUtils.js b/src/mesh/meshUtils.js index 5356553..ef12b06 100644 --- a/src/mesh/meshUtils.js +++ b/src/mesh/meshUtils.js @@ -8,12 +8,9 @@ */ // Internal imports -import { BasisFunctions, triangleBoundarySides } from "./basisFunctions.js"; +import { getBoundarySides } from "../elements/basisFunctions.js"; import { Mesh1D, Mesh2D } from "./meshGeneration.js"; -import { NumericalIntegration } from "../numerics/numericalIntegration.js"; -import { SparseMatrix } from "../numerics/sparseMatrix.js"; import { debugLog, errorLog } from "../utilities/logging.js"; - /** * Function to prepare the mesh for finite element analysis * @param {object} meshConfig - Object containing computational mesh details @@ -96,317 +93,6 @@ export function prepareMesh(meshConfig) { }; } -/** - * Function to initialize the FEA matrices and numerical tools - * @param {object} meshData - Object containing mesh data from prepareMesh() - * @returns {object} An object containing initialized matrices and numerical tools - */ -export function initializeFEA(meshData) { - const { totalNodes, nop, meshDimension, elementOrder, elementShape } = meshData; - - // Initialize variables for matrix assembly - let residualVector = new Array(totalNodes).fill(0); - let jacobianMatrix = new SparseMatrix(totalNodes); - let localToGlobalMap = []; - - // Initialize the BasisFunctions class - const basisFunctions = new BasisFunctions({ - meshDimension, - elementOrder, - elementShape, - }); - - // Initialize the NumericalIntegration class - const numericalIntegration = new NumericalIntegration({ - meshDimension, - elementOrder, - elementShape, - }); - - // Calculate Gauss points and weights - let gaussPointsAndWeights = numericalIntegration.getGaussPointsAndWeights(); - let gaussPoints = gaussPointsAndWeights.gaussPoints; - let gaussWeights = gaussPointsAndWeights.gaussWeights; - - // Integration points for 2D elements (tensor product for quads, collapsed rule for triangles) - const integrationPoints2D = - meshDimension === "2D" ? getIntegrationPoints2D(gaussPoints, gaussWeights, elementShape) : []; - - // Basis functions on the reference element at each integration point (the same for every element) - const referenceBasisFunctions = - meshDimension === "1D" - ? gaussPoints.map((gaussPoint) => basisFunctions.getBasisFunctions(gaussPoint)) - : integrationPoints2D.map(([ksi, eta]) => basisFunctions.getBasisFunctions(ksi, eta)); - - // Determine the number of nodes in the reference element based on the first element in the nop array - const nodesPerElement = nop[0].length; - - return { - residualVector, - jacobianMatrix, - localToGlobalMap, - basisFunctions, - gaussPoints, - gaussWeights, - integrationPoints2D, - referenceBasisFunctions, - nodesPerElement, - }; -} - -/** - * Function to build the list of 2D integration points with their combined weights - * @param {array} gaussPoints - Gauss points from NumericalIntegration - * @param {array} gaussWeights - Gauss weights from NumericalIntegration - * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' - * @returns {array} Array of [ksi, eta, weight] entries - */ -export function getIntegrationPoints2D(gaussPoints, gaussWeights, elementShape = "quad") { - const integrationPoints = []; - - if (elementShape === "triangle") { - // Triangle rules already provide (ksi, eta) pairs with their weights - for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { - const [ksi, eta] = gaussPoints[gaussPointIndex]; - integrationPoints.push([ksi, eta, gaussWeights[gaussPointIndex]]); - } - } else { - // Quadrilaterals use the tensor product of the 1D Gauss rule - for (let gaussPointIndex1 = 0; gaussPointIndex1 < gaussPoints.length; gaussPointIndex1++) { - for (let gaussPointIndex2 = 0; gaussPointIndex2 < gaussPoints.length; gaussPointIndex2++) { - integrationPoints.push([ - gaussPoints[gaussPointIndex1], - gaussPoints[gaussPointIndex2], - gaussWeights[gaussPointIndex1] * gaussWeights[gaussPointIndex2], - ]); - } - } - } - - return integrationPoints; -} - -/** - * Function to perform isoparametric mapping for 1D elements - * @param {object} params - Parameters for the mapping. An optional `basisFunctionDerivX` array can be - * passed to be filled in place, to avoid allocating a new array at every integration point - * @returns {object} An object containing the mapped data - */ -export function performIsoparametricMapping1D(params) { - const { - basisFunction, - basisFunctionDerivKsi, - nodesXCoordinates, - localToGlobalMap, - nodesPerElement, - basisFunctionDerivX = [], - } = params; - - let xCoordinates = 0; - let ksiDerivX = 0; - - // Isoparametric mapping - for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { - xCoordinates += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; - ksiDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; - } - let detJacobian = ksiDerivX; - - // Compute x-derivative of basis functions - for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { - basisFunctionDerivX[localNodeIndex] = basisFunctionDerivKsi[localNodeIndex] / detJacobian; - } - - return { - xCoordinates, - detJacobian, - basisFunctionDerivX, - }; -} - -/** - * Function to perform isoparametric mapping for 2D elements - * @param {object} params - Parameters for the mapping. Optional `basisFunctionDerivX` and - * `basisFunctionDerivY` arrays can be passed to be filled in place, to avoid allocating new arrays at - * every integration point - * @returns {object} An object containing the mapped data - */ -export function performIsoparametricMapping2D(params) { - const { - basisFunction, - basisFunctionDerivKsi, - basisFunctionDerivEta, - nodesXCoordinates, - nodesYCoordinates, - localToGlobalMap, - nodesPerElement, - basisFunctionDerivX = [], - basisFunctionDerivY = [], - } = params; - - let xCoordinates = 0; - let yCoordinates = 0; - let ksiDerivX = 0; - let etaDerivX = 0; - let ksiDerivY = 0; - let etaDerivY = 0; - - // Isoparametric mapping - for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { - xCoordinates += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; - yCoordinates += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex]; - ksiDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; - etaDerivX += nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex]; - ksiDerivY += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex]; - etaDerivY += nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex]; - } - let detJacobian = ksiDerivX * etaDerivY - etaDerivX * ksiDerivY; - - // Compute x-derivative and y-derivative of basis functions - for (let localNodeIndex = 0; localNodeIndex < nodesPerElement; localNodeIndex++) { - // The x-derivative of the n basis function - basisFunctionDerivX[localNodeIndex] = - (etaDerivY * basisFunctionDerivKsi[localNodeIndex] - - ksiDerivY * basisFunctionDerivEta[localNodeIndex]) / - detJacobian; - // The y-derivative of the n basis function - basisFunctionDerivY[localNodeIndex] = - (ksiDerivX * basisFunctionDerivEta[localNodeIndex] - - etaDerivX * basisFunctionDerivKsi[localNodeIndex]) / - detJacobian; - } - - return { - xCoordinates, - yCoordinates, - detJacobian, - basisFunctionDerivX, - basisFunctionDerivY, - }; -} - -/** - * Function to test if a point is inside a triangle using barycentric coordinates, - * also returning the natural coordinates (ksi, eta) - * @param {number} x - X-coordinate of the point - * @param {number} y - Y-coordinate of the point - * @param {array} vertices - Triangle vertices [[x0,y0],[x1,y1],[x2,y2]] - * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} - */ -export function pointInsideTriangle(x, y, vertices) { - const tolerance = 1e-12; - const [v0, v1, v2] = vertices; - - const denom = (v1[1] - v2[1]) * (v0[0] - v2[0]) + (v2[0] - v1[0]) * (v0[1] - v2[1]); - - const ksi = ((v1[1] - v2[1]) * (x - v2[0]) + (v2[0] - v1[0]) * (y - v2[1])) / denom; - const eta = ((v2[1] - v0[1]) * (x - v2[0]) + (v0[0] - v2[0]) * (y - v2[1])) / denom; - const gamma = 1 - ksi - eta; - - const inside = ksi >= -tolerance && eta >= -tolerance && gamma >= -tolerance; - return { inside, ksi, eta }; -} - -/** - * Function to test if a point is inside a triangular element, also returning the natural coordinates - * (ksi, eta) of the FEAScript reference triangle - * @param {number} x - X-coordinate of the point - * @param {number} y - Y-coordinate of the point - * @param {array} vertices - Triangle vertices in FEAScript order [[x0,y0],[x1,y1],[x2,y2]] - * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} - */ -export function pointInsideTriangularElement(x, y, vertices) { - // Vertices order (reference coordinates): - // 1 (0,1) - // | \ - // | \ - // 0 --- 2 (1,0) - // The barycentric weights of vertices 2 and 1 are the natural coordinates ksi and eta - const { inside, ksi: weight0, eta: weight1 } = pointInsideTriangle(x, y, vertices); - return { inside, ksi: 1 - weight0 - weight1, eta: weight1 }; -} - -/** - * Function to test if a point is inside a quadrilateral by spliting it into triangles and using barycentric coordinates - * @param {number} x - X-coordinate of the point - * @param {number} y - Y-coordinate of the point - * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]] - * @returns {object} Object containing inside boolean and natural coordinates {inside, ksi, eta} - */ -export function pointInsideQuadrilateral(x, y, vertices) { - const [firstTriangleVertices, secondTriangleVertices] = splitQuadrilateral(vertices); - const pointInsideFirstTriangle = pointInsideTriangle(x, y, firstTriangleVertices); - const pointInsideSecondTriangle = pointInsideTriangle(x, y, secondTriangleVertices); - - const inside = pointInsideFirstTriangle.inside || pointInsideSecondTriangle.inside; - let ksi = 0; - let eta = 0; - - if (inside) { - const [v0, v1, v2, v3] = vertices; - - // Function to calculate distance from point to line segment - const getDistanceFromLine = (p1, p2) => { - const num = Math.abs((p2[0] - p1[0]) * (p1[1] - y) - (p1[0] - x) * (p2[1] - p1[1])); - const den = Math.sqrt((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2); - return num / den; - }; - - // Calculate distances to edges based on vertex order: - // 1 (v1) --- 3 (v3) - // | | - // 0 (v0) --- 2 (v2) - - const distLeft = getDistanceFromLine(v0, v1); - const distRight = getDistanceFromLine(v2, v3); - const distBottom = getDistanceFromLine(v0, v2); - const distTop = getDistanceFromLine(v1, v3); - - ksi = distLeft / (distLeft + distRight); - eta = distBottom / (distBottom + distTop); - - // The distance ratios are exact only for parallelograms, so refine them with Newton iterations - // on the bilinear mapping x(ksi, eta) = v0 (1-ksi)(1-eta) + v1 (1-ksi) eta + v2 ksi (1-eta) + v3 ksi eta - for (let iteration = 0; iteration < 20; iteration++) { - const mappedX = - v0[0] * (1 - ksi) * (1 - eta) + v1[0] * (1 - ksi) * eta + v2[0] * ksi * (1 - eta) + v3[0] * ksi * eta; - const mappedY = - v0[1] * (1 - ksi) * (1 - eta) + v1[1] * (1 - ksi) * eta + v2[1] * ksi * (1 - eta) + v3[1] * ksi * eta; - const ksiDerivX = (v2[0] - v0[0]) * (1 - eta) + (v3[0] - v1[0]) * eta; - const ksiDerivY = (v2[1] - v0[1]) * (1 - eta) + (v3[1] - v1[1]) * eta; - const etaDerivX = (v1[0] - v0[0]) * (1 - ksi) + (v3[0] - v2[0]) * ksi; - const etaDerivY = (v1[1] - v0[1]) * (1 - ksi) + (v3[1] - v2[1]) * ksi; - const detJacobian = ksiDerivX * etaDerivY - etaDerivX * ksiDerivY; - const residualX = x - mappedX; - const residualY = y - mappedY; - const deltaKsi = (etaDerivY * residualX - etaDerivX * residualY) / detJacobian; - const deltaEta = (ksiDerivX * residualY - ksiDerivY * residualX) / detJacobian; - ksi += deltaKsi; - eta += deltaEta; - if (Math.abs(deltaKsi) + Math.abs(deltaEta) < 1e-14) break; - } - } - - return { inside, ksi, eta }; -} - -/** - * Function to split the quadrilateral elements into two triangles - * @param {array} vertices - Quadrilateral vertices [[x0,y0],[x1,y1],[x2,y2],[x3,y3]] - * @returns {array} Array of triangle vertices: [[v0,v1,v3], [v0,v2,v3]] - */ -export function splitQuadrilateral(vertices) { - const [v0, v1, v2, v3] = vertices; - // Vertices order: - // 1 --- 3 - // | | - // 0 --- 2 - return [ - [v0, v1, v3], - [v0, v2, v3], - ]; -} - /** * Function that finds the list of adjacent elements for each node in the mesh * @param {object} meshData - Object containing nodal numbering (NOP) @@ -448,7 +134,6 @@ export function getBoundarySegments(meshData) { let boundaryLineElements = []; let boundaryNodesSegments = []; let boundaryGlobalElementIndex = 0; - let boundarySides; const { nodesXCoordinates, nodesYCoordinates, @@ -458,39 +143,7 @@ export function getBoundarySegments(meshData) { elementOrder, elementShape, } = meshData; - - if (meshDimension === "1D") { - if (elementOrder === "linear") { - boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - } else if (elementOrder === "quadratic") { - boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - } - } else if (meshDimension === "2D" && elementShape === "triangle") { - // Nodes at the bottom (0), left (1) and hypotenuse (2) sides of the reference triangle - boundarySides = triangleBoundarySides[elementOrder]; - } else if (meshDimension === "2D") { - if (elementOrder === "linear") { - boundarySides = { - 0: [0, 2], // Nodes at the bottom side of the reference element - 1: [0, 1], // Nodes at the left side of the reference element - 2: [1, 3], // Nodes at the top side of the reference element - 3: [2, 3], // Nodes at the right side of the reference element - }; - } else if (elementOrder === "quadratic") { - boundarySides = { - 0: [0, 3, 6], // Nodes at the bottom side of the reference element - 1: [0, 1, 2], // Nodes at the left side of the reference element - 2: [2, 5, 8], // Nodes at the top side of the reference element - 3: [6, 7, 8], // Nodes at the right side of the reference element - }; - } - } + const boundarySides = getBoundarySides(meshDimension, elementOrder, elementShape); // Iterate over all boundaries for (let boundaryIndex = 0; boundaryIndex < boundaryElements.length; boundaryIndex++) { diff --git a/src/models/genericBoundaryConditions.js b/src/models/genericBoundaryConditions.js deleted file mode 100644 index b017645..0000000 --- a/src/models/genericBoundaryConditions.js +++ /dev/null @@ -1,283 +0,0 @@ -/** - * ════════════════════════════════════════════════════════════════ - * FEAScript Core Library - * Finite Element Simulation in JavaScript - * Version: 0.3.0 (RC) | https://feascript.com - * MIT License © 2023–2026 FEAScript - * ════════════════════════════════════════════════════════════════ - */ - -// Internal imports -import { triangleBoundarySides } from "../mesh/basisFunctions.js"; -import { debugLog } from "../utilities/logging.js"; - -/** - * Class to handle generic boundary conditions application - */ -export class GenericBoundaryConditions { - /** - * Constructor to initialize the GenericBoundaryConditions class - * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis - * @param {array} boundaryElements - Array containing elements that belong to each boundary - * @param {array} nop - Nodal numbering (NOP) array representing the connectivity between elements and nodes - * @param {string} meshDimension - The dimension of the mesh (e.g., "2D") - * @param {string} elementOrder - The order of elements (e.g., "linear", "quadratic") - * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' - */ - constructor(boundaryConditions, boundaryElements, nop, meshDimension, elementOrder, elementShape = "quad") { - this.boundaryConditions = boundaryConditions; - this.boundaryElements = boundaryElements; - this.nop = nop; - this.meshDimension = meshDimension; - this.elementOrder = elementOrder; - this.elementShape = elementShape; - } - - /** - * Function to impose Dirichlet boundary conditions - * @param {array} residualVector - The residual vector to be modified - * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified - * @param {array} [solutionVector] - Current solution (Newton-Raphson iterate); when provided, the - * residual is set to the increment needed to reach the prescribed value instead of the value itself - * - * For consistency across both linear and nonlinear formulations, - * this project always refers 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. - */ - imposeDirichletBoundaryConditions(residualVector, jacobianMatrix, solutionVector) { - if (this.meshDimension === "1D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantValue") { - const value = this.boundaryConditions[boundaryKey][1]; - debugLog(`Boundary ${boundaryKey}: Applying constant value of ${value} (Dirichlet condition)`); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson - residualVector[globalNodeIndex] = solutionVector - ? value - solutionVector[globalNodeIndex] - : value; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [2], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson - residualVector[globalNodeIndex] = solutionVector - ? value - solutionVector[globalNodeIndex] - : value; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } - }); - } - }); - } else if (this.meshDimension === "2D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantValue") { - const value = this.boundaryConditions[boundaryKey][1]; - debugLog(`Boundary ${boundaryKey}: Applying constant value of ${value} (Dirichlet condition)`); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementShape === "triangle") { - triangleBoundarySides[this.elementOrder][side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson - residualVector[globalNodeIndex] = solutionVector - ? value - solutionVector[globalNodeIndex] - : value; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0, 2], // Nodes at the bottom side of the reference element - 1: [0, 1], // Nodes at the left side of the reference element - 2: [1, 3], // Nodes at the top side of the reference element - 3: [2, 3], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson - residualVector[globalNodeIndex] = solutionVector - ? value - solutionVector[globalNodeIndex] - : value; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0, 3, 6], // Nodes at the bottom side of the reference element - 1: [0, 1, 2], // Nodes at the left side of the reference element - 2: [2, 5, 8], // Nodes at the top side of the reference element - 3: [6, 7, 8], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson - residualVector[globalNodeIndex] = solutionVector - ? value - solutionVector[globalNodeIndex] - : value; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } - }); - } - }); - } - } - - /** - * Function to impose constant value (Dirichlet) boundary conditions for the frontal solver - * @param {array} nodeConstraintCode - Array indicating boundary condition code for each node - * @param {array} boundaryValues - Array containing boundary condition values - */ - imposeConstantValueBoundaryConditionsFront(nodeConstraintCode, boundaryValues) { - if (this.meshDimension === "1D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantValue") { - const value = this.boundaryConditions[boundaryKey][1]; - debugLog(`Boundary ${boundaryKey}: Applying constant value of ${value} (Dirichlet condition)`); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = value; - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [2], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = value; - }); - } - }); - } - }); - } else if (this.meshDimension === "2D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantValue") { - const value = this.boundaryConditions[boundaryKey][1]; - debugLog(`Boundary ${boundaryKey}: Applying constant value of ${value} (Dirichlet condition)`); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementShape === "triangle") { - triangleBoundarySides[this.elementOrder][side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = value; - }); - } else if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0, 2], // Nodes at the bottom side of the reference element - 1: [0, 1], // Nodes at the left side of the reference element - 2: [1, 3], // Nodes at the top side of the reference element - 3: [2, 3], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = value; - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0, 3, 6], // Nodes at the bottom side of the reference element - 1: [0, 1, 2], // Nodes at the left side of the reference element - 2: [2, 5, 8], // Nodes at the top side of the reference element - 3: [6, 7, 8], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = value; - }); - } - }); - } - }); - } - } -} diff --git a/src/models/thermalBoundaryConditions.js b/src/models/thermalBoundaryConditions.js deleted file mode 100644 index a1a021e..0000000 --- a/src/models/thermalBoundaryConditions.js +++ /dev/null @@ -1,989 +0,0 @@ -/** - * ════════════════════════════════════════════════════════════════ - * FEAScript Core Library - * Finite Element Simulation in JavaScript - * Version: 0.3.0 (RC) | https://feascript.com - * MIT License © 2023–2026 FEAScript - * ════════════════════════════════════════════════════════════════ - */ - -// Internal imports -import { triangleBoundarySides } from "../mesh/basisFunctions.js"; -import { NumericalIntegration } from "../numerics/numericalIntegration.js"; -import { debugLog } from "../utilities/logging.js"; - -/** - * Class to handle thermal boundary conditions application - */ -export class ThermalBoundaryConditions { - /** - * Constructor to initialize the ThermalBoundaryConditions class - * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis - * @param {array} boundaryElements - Array containing elements that belong to each boundary - * @param {array} nop - Nodal numbering (NOP) array representing the connectivity between elements and nodes - * @param {string} meshDimension - The dimension of the mesh (e.g., "2D") - * @param {string} elementOrder - The order of elements (e.g., "linear", "quadratic") - * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' - */ - constructor(boundaryConditions, boundaryElements, nop, meshDimension, elementOrder, elementShape = "quad") { - this.boundaryConditions = boundaryConditions; - this.boundaryElements = boundaryElements; - this.nop = nop; - this.meshDimension = meshDimension; - this.elementOrder = elementOrder; - this.elementShape = elementShape; - } - - /** - * Function to compute the convection contributions on one side of a triangular element - * @param {number} elementIndex - Index of the element - * @param {number} side - Side of the reference triangle (0: bottom, 1: left, 2: hypotenuse) - * @param {number} convectionHeatTransferCoefficient - Convection heat transfer coefficient - * @param {number} externalTemperature - External temperature - * @param {array} nodesXCoordinates - Array of x-coordinates of nodes - * @param {array} nodesYCoordinates - Array of y-coordinates of nodes - * @param {object} basisFunctions - Object containing basis functions and their derivatives - * @returns {object} An object containing: - * - sideNodeIndices: Local node indices on the side - * - localJacobianMatrix: Contributions indexed by local node indices - * - localResidualVector: Contributions indexed by local node indices - */ - computeTriangleConvectionContributions( - elementIndex, - side, - convectionHeatTransferCoefficient, - externalTemperature, - nodesXCoordinates, - nodesYCoordinates, - basisFunctions, - ) { - const nodesPerElement = this.nop[elementIndex].length; - const sideNodeIndices = triangleBoundarySides[this.elementOrder][side]; - const localJacobianMatrix = Array(nodesPerElement) - .fill() - .map(() => Array(nodesPerElement).fill(0)); - const localResidualVector = Array(nodesPerElement).fill(0); - - // Parametrize the side with t in [0, 1]: (ksi, eta) = start + t * direction - const sideParametrization = { - 0: { start: [0, 0], direction: [1, 0] }, // Bottom side (eta = 0) - 1: { start: [0, 0], direction: [0, 1] }, // Left side (ksi = 0) - 2: { start: [0, 1], direction: [1, -1] }, // Hypotenuse (ksi + eta = 1) - }; - const { start, direction } = sideParametrization[side]; - - // 3-point Gauss rule on [0, 1], exact for the boundary mass terms of both linear and quadratic triangles - const { gaussPoints, gaussWeights } = new NumericalIntegration({ - meshDimension: "1D", - elementOrder: "quadratic", - }).getGaussPointsAndWeights(); - - for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { - const t = gaussPoints[gaussPointIndex]; - const { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta } = - basisFunctions.getBasisFunctions(start[0] + t * direction[0], start[1] + t * direction[1]); - - // Tangent vector of the side in physical space (chain rule along the side parametrization) - let tangentX = 0; - let tangentY = 0; - for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - const basisFunctionDerivT = - basisFunctionDerivKsi[nodeIndex] * direction[0] + basisFunctionDerivEta[nodeIndex] * direction[1]; - tangentX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivT; - tangentY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivT; - } - const tangentVectorLength = Math.sqrt(tangentX ** 2 + tangentY ** 2); - - for (const localNodeIndex of sideNodeIndices) { - localResidualVector[localNodeIndex] += - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - convectionHeatTransferCoefficient * - externalTemperature; - - for (const localNodeIndex2 of sideNodeIndices) { - localJacobianMatrix[localNodeIndex][localNodeIndex2] += - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - basisFunction[localNodeIndex2] * - convectionHeatTransferCoefficient; - } - } - } - - return { sideNodeIndices, localJacobianMatrix, localResidualVector }; - } - - /** - * Function to impose constant temperature boundary conditions (Dirichlet type) - * @param {array} residualVector - The residual vector to be modified - * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified - * - * For consistency across both linear and nonlinear formulations, - * this project always refers 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. - */ - imposeConstantTempBoundaryConditions(residualVector, jacobianMatrix) { - if (this.meshDimension === "1D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantTemperature") { - const constantTemperatureValue = this.boundaryConditions[boundaryKey][1]; - debugLog( - `Boundary ${boundaryKey}: Applying constant temperature of ${constantTemperatureValue} K (Dirichlet condition)`, - ); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the ConstantTemp value - residualVector[globalNodeIndex] = constantTemperatureValue; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [2], // Node at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the ConstantTemp value - residualVector[globalNodeIndex] = constantTemperatureValue; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } - }); - } - }); - } else if (this.meshDimension === "2D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantTemperature") { - const constantTemperatureValue = this.boundaryConditions[boundaryKey][1]; - debugLog( - `Boundary ${boundaryKey}: Applying constant temperature of ${constantTemperatureValue} K (Dirichlet condition)`, - ); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementShape === "triangle") { - triangleBoundarySides[this.elementOrder][side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the ConstantTemp value - residualVector[globalNodeIndex] = constantTemperatureValue; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0, 2], // Nodes at the bottom side of the reference element - 1: [0, 1], // Nodes at the left side of the reference element - 2: [1, 3], // Nodes at the top side of the reference element - 3: [2, 3], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the ConstantTemp value - residualVector[globalNodeIndex] = constantTemperatureValue; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0, 3, 6], // Nodes at the bottom side of the reference element - 1: [0, 1, 2], // Nodes at the left side of the reference element - 2: [2, 5, 8], // Nodes at the top side of the reference element - 3: [6, 7, 8], // Nodes at the right side of the reference element - }; - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - // Set the residual vector to the ConstantTemp value - residualVector[globalNodeIndex] = constantTemperatureValue; - // Set the Jacobian matrix row to zero - jacobianMatrix.clearRow(globalNodeIndex); - // Set the diagonal entry of the Jacobian matrix to one - jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); - }); - } - }); - } - }); - } - } - - /** - * Function to impose constant temperature boundary conditions for the frontal solver - * @param {array} nodeConstraintCode - Array indicating boundary condition code for each node - * @param {array} boundaryValues - Array containing boundary condition values - */ - imposeConstantTempBoundaryConditionsFront(nodeConstraintCode, boundaryValues) { - if (this.meshDimension === "1D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantTemperature") { - const constantTemperatureValue = this.boundaryConditions[boundaryKey][1]; - debugLog( - `Boundary ${boundaryKey}: Applying constant temperature of ${constantTemperatureValue} K (Dirichlet condition)`, - ); - - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [1], // Node at the right side of the reference element - }; - - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - - // Set boundary condition code and value - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = constantTemperatureValue; - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0], // Node at the left side of the reference element - 1: [2], // Node at the right side of the reference element - }; - - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - - // Set boundary condition code and value - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = constantTemperatureValue; - }); - } - }); - } - }); - } else if (this.meshDimension === "2D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "constantTemperature") { - const constantTemperatureValue = this.boundaryConditions[boundaryKey][1]; - debugLog( - `Boundary ${boundaryKey}: Applying constant temperature of ${constantTemperatureValue} K (Dirichlet condition)`, - ); - - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementShape === "triangle") { - triangleBoundarySides[this.elementOrder][side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - - // Set boundary condition code and value - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = constantTemperatureValue; - }); - } else if (this.elementOrder === "linear") { - const boundarySides = { - 0: [0, 2], // Nodes at the bottom side of the reference element - 1: [0, 1], // Nodes at the left side of the reference element - 2: [1, 3], // Nodes at the top side of the reference element - 3: [2, 3], // Nodes at the right side of the reference element - }; - - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - - // Set boundary condition code and value - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = constantTemperatureValue; - }); - } else if (this.elementOrder === "quadratic") { - const boundarySides = { - 0: [0, 3, 6], // Nodes at the bottom side of the reference element - 1: [0, 1, 2], // Nodes at the left side of the reference element - 2: [2, 5, 8], // Nodes at the top side of the reference element - 3: [6, 7, 8], // Nodes at the right side of the reference element - }; - - boundarySides[side].forEach((nodeIndex) => { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - - // Set boundary condition code and value - nodeConstraintCode[globalNodeIndex] = 1; - boundaryValues[globalNodeIndex] = constantTemperatureValue; - }); - } - }); - } - }); - } - } - - /** - * Function to impose convection boundary conditions (Robin type) - * @param {array} residualVector - The residual vector to be modified - * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified - * @param {array} gaussPoints - Array of Gauss points for numerical integration - * @param {array} gaussWeights - Array of Gauss weights for numerical integration - * @param {array} nodesXCoordinates - Array of x-coordinates of nodes - * @param {array} nodesYCoordinates - Array of y-coordinates of nodes - * @param {object} basisFunctions - Object containing basis functions and their derivatives - */ - imposeConvectionBoundaryConditions( - residualVector, - jacobianMatrix, - gaussPoints, - gaussWeights, - nodesXCoordinates, - nodesYCoordinates, - basisFunctions, - ) { - // Extract convection parameters from boundary conditions - let convectionHeatTransferCoefficients = []; - let convectionExternalTemperatures = []; - Object.keys(this.boundaryConditions).forEach((key) => { - const boundaryCondition = this.boundaryConditions[key]; - if (boundaryCondition[0] === "convection") { - convectionHeatTransferCoefficients[key] = boundaryCondition[1]; - convectionExternalTemperatures[key] = boundaryCondition[2]; - } - }); - - if (this.meshDimension === "1D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "convection") { - const convectionHeatTransferCoefficient = convectionHeatTransferCoefficients[boundaryKey]; - const externalTemperature = convectionExternalTemperatures[boundaryKey]; - debugLog( - `Boundary ${boundaryKey}: Applying convection with heat transfer coefficient h=${convectionHeatTransferCoefficient} W/(m²·K) and external temperature T∞=${externalTemperature} K`, - ); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - let nodeIndex; - if (this.elementOrder === "linear") { - if (side === 0) { - // Node at the left side of the reference element - nodeIndex = 0; - } else { - // Node at the right side of the reference element - nodeIndex = 1; - } - } else if (this.elementOrder === "quadratic") { - if (side === 0) { - // Node at the left side of the reference element - nodeIndex = 0; - } else { - // Node at the right side of the reference element - nodeIndex = 2; - } - } - - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - debugLog( - ` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - residualVector[globalNodeIndex] += -convectionHeatTransferCoefficient * externalTemperature; - jacobianMatrix.addValue(globalNodeIndex, globalNodeIndex, convectionHeatTransferCoefficient); - }); - } - }); - } else if (this.meshDimension === "2D") { - Object.keys(this.boundaryConditions).forEach((boundaryKey) => { - if (this.boundaryConditions[boundaryKey][0] === "convection") { - const convectionHeatTransferCoefficient = convectionHeatTransferCoefficients[boundaryKey]; - const externalTemperature = convectionExternalTemperatures[boundaryKey]; - debugLog( - `Boundary ${boundaryKey}: Applying convection with heat transfer coefficient h=${convectionHeatTransferCoefficient} W/(m²·K) and external temperature T∞=${externalTemperature} K`, - ); - this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { - if (this.elementShape === "triangle") { - const { sideNodeIndices, localJacobianMatrix, localResidualVector } = - this.computeTriangleConvectionContributions( - elementIndex, - side, - convectionHeatTransferCoefficient, - externalTemperature, - nodesXCoordinates, - nodesYCoordinates, - basisFunctions, - ); - - for (const localNodeIndex of sideNodeIndices) { - const globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1; - debugLog( - ` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${localNodeIndex + 1})`, - ); - residualVector[globalNodeIndex] += localResidualVector[localNodeIndex]; - - for (const localNodeIndex2 of sideNodeIndices) { - const globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1; - jacobianMatrix.addValue( - globalNodeIndex, - globalNodeIndex2, - localJacobianMatrix[localNodeIndex][localNodeIndex2], - ); - } - } - } else if (this.elementOrder === "linear") { - let gaussPoint1, gaussPoint2, firstNodeIndex, lastNodeIndex, nodeIncrement; - if (side === 0) { - // Nodes at the bottom side of the reference element - gaussPoint1 = gaussPoints[0]; - gaussPoint2 = 0; - firstNodeIndex = 0; - lastNodeIndex = 3; - nodeIncrement = 2; - } else if (side === 1) { - // Nodes at the left side of the reference element - gaussPoint1 = 0; - gaussPoint2 = gaussPoints[0]; - firstNodeIndex = 0; - lastNodeIndex = 2; - nodeIncrement = 1; - } else if (side === 2) { - // Nodes at the top side of the reference element - gaussPoint1 = gaussPoints[0]; - gaussPoint2 = 1; - firstNodeIndex = 1; - lastNodeIndex = 4; - nodeIncrement = 2; - } else if (side === 3) { - // Nodes at the right side of the reference element - gaussPoint1 = 1; - gaussPoint2 = gaussPoints[0]; - firstNodeIndex = 2; - lastNodeIndex = 4; - nodeIncrement = 1; - } - - let basisFunctionsAndDerivatives = basisFunctions.getBasisFunctions(gaussPoint1, gaussPoint2); - let basisFunction = basisFunctionsAndDerivatives.basisFunction; - let basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi; - let basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta; - - let ksiDerivX = 0; - let ksiDerivY = 0; - let etaDerivX = 0; - let etaDerivY = 0; - const nodesPerElement = this.nop[elementIndex].length; - for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - - // For boundaries along Ksi (horizontal), use Ksi derivatives - if (side === 0 || side === 2) { - ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - } - // For boundaries along Eta (vertical), use Eta derivatives - else if (side === 1 || side === 3) { - etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - } - } - - // Compute the length of tangent vector - let tangentVectorLength; - if (side === 0 || side === 2) { - tangentVectorLength = Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2); - } else { - tangentVectorLength = Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2); - } - - for ( - let localNodeIndex = firstNodeIndex; - localNodeIndex < lastNodeIndex; - localNodeIndex += nodeIncrement - ) { - let globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1; - debugLog( - ` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${localNodeIndex + 1})`, - ); - - // Apply boundary condition with proper Jacobian for all sides - residualVector[globalNodeIndex] += - -gaussWeights[0] * - tangentVectorLength * - basisFunction[localNodeIndex] * - convectionHeatTransferCoefficient * - externalTemperature; - - for ( - let localNodeIndex2 = firstNodeIndex; - localNodeIndex2 < lastNodeIndex; - localNodeIndex2 += nodeIncrement - ) { - let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1; - jacobianMatrix.addValue( - globalNodeIndex, - globalNodeIndex2, - -gaussWeights[0] * - tangentVectorLength * - basisFunction[localNodeIndex] * - basisFunction[localNodeIndex2] * - convectionHeatTransferCoefficient, - ); - } - } - } else if (this.elementOrder === "quadratic") { - for (let gaussPointIndex = 0; gaussPointIndex < 3; gaussPointIndex++) { - let gaussPoint1, gaussPoint2, firstNodeIndex, lastNodeIndex, nodeIncrement; - if (side === 0) { - // Nodes at the bottom side of the reference element - gaussPoint1 = gaussPoints[gaussPointIndex]; - gaussPoint2 = 0; - firstNodeIndex = 0; - lastNodeIndex = 7; - nodeIncrement = 3; - } else if (side === 1) { - // Nodes at the left side of the reference element - gaussPoint1 = 0; - gaussPoint2 = gaussPoints[gaussPointIndex]; - firstNodeIndex = 0; - lastNodeIndex = 3; - nodeIncrement = 1; - } else if (side === 2) { - // Nodes at the top side of the reference element - gaussPoint1 = gaussPoints[gaussPointIndex]; - gaussPoint2 = 1; - firstNodeIndex = 2; - lastNodeIndex = 9; - nodeIncrement = 3; - } else if (side === 3) { - // Nodes at the right side of the reference element - gaussPoint1 = 1; - gaussPoint2 = gaussPoints[gaussPointIndex]; - firstNodeIndex = 6; - lastNodeIndex = 9; - nodeIncrement = 1; - } - let basisFunctionsAndDerivatives = basisFunctions.getBasisFunctions(gaussPoint1, gaussPoint2); - let basisFunction = basisFunctionsAndDerivatives.basisFunction; - let basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi; - let basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta; - - let ksiDerivX = 0; - let ksiDerivY = 0; - let etaDerivX = 0; - let etaDerivY = 0; - const nodesPerElement = this.nop[elementIndex].length; - for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - - // For boundaries along Ksi (horizontal), use Ksi derivatives - if (side === 0 || side === 2) { - ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - } - // For boundaries along Eta (vertical), use Eta derivatives - else if (side === 1 || side === 3) { - etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - } - } - - // Compute the length of tangent vector - let tangentVectorLength; - if (side === 0 || side === 2) { - tangentVectorLength = Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2); - } else { - tangentVectorLength = Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2); - } - - for ( - let localNodeIndex = firstNodeIndex; - localNodeIndex < lastNodeIndex; - localNodeIndex += nodeIncrement - ) { - let globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1; - debugLog( - ` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${localNodeIndex + 1})`, - ); - - // Apply boundary condition with proper Jacobian for all sides - residualVector[globalNodeIndex] += - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - convectionHeatTransferCoefficient * - externalTemperature; - - for ( - let localNodeIndex2 = firstNodeIndex; - localNodeIndex2 < lastNodeIndex; - localNodeIndex2 += nodeIncrement - ) { - let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1; - jacobianMatrix.addValue( - globalNodeIndex, - globalNodeIndex2, - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - basisFunction[localNodeIndex2] * - convectionHeatTransferCoefficient, - ); - } - } - } - } - }); - } - }); - } - } - - /** - * Function to impose convection boundary conditions for the frontal solver - * @param {number} elementIndex - Index of the element being processed - * @param {array} nodesXCoordinates - Array of x-coordinates of nodes - * @param {array} nodesYCoordinates - Array of y-coordinates of nodes - * @param {array} gaussPoints - Array of Gauss points for numerical integration - * @param {array} gaussWeights - Array of Gauss weights for numerical integration - * @param {object} basisFunctions - Object containing basis functions and their derivatives - * @returns {object} An object containing: - * - localJacobianMatrix: Local Jacobian matrix with convection contributions - * - localResidualVector: Residual vector with convection contributions - */ - imposeConvectionBoundaryConditionsFront( - elementIndex, - nodesXCoordinates, - nodesYCoordinates, - gaussPoints, - gaussWeights, - basisFunctions, - ) { - // Extract convection parameters from boundary conditions - let convectionHeatTransferCoefficients = []; - let convectionExternalTemperatures = []; - Object.keys(this.boundaryConditions).forEach((key) => { - const boundaryCondition = this.boundaryConditions[key]; - if (boundaryCondition[0] === "convection") { - convectionHeatTransferCoefficients[key] = boundaryCondition[1]; - convectionExternalTemperatures[key] = boundaryCondition[2]; - } - }); - - // Initialize local Jacobian matrix and local residual vector - const nodesPerElement = this.nop[elementIndex].length; - const localJacobianMatrix = Array(nodesPerElement) - .fill() - .map(() => Array(nodesPerElement).fill(0)); - const localResidualVector = Array(nodesPerElement).fill(0); - - // Check if this element is on a convection boundary - for (const boundaryKey in this.boundaryElements) { - if (this.boundaryConditions[boundaryKey]?.[0] === "convection") { - const convectionHeatTransferCoefficient = convectionHeatTransferCoefficients[boundaryKey]; - const externalTemperature = convectionExternalTemperatures[boundaryKey]; - debugLog( - `Boundary ${boundaryKey}: Applying convection with heat transfer coefficient h=${convectionHeatTransferCoefficient} W/(m²·K) and external temperature T∞=${externalTemperature} K`, - ); - - // Find if this element is on this boundary and which side - const boundaryElement = this.boundaryElements[boundaryKey].find( - ([boundaryElementIndex]) => boundaryElementIndex === elementIndex, - ); - - if (boundaryElement && this.meshDimension === "2D" && this.elementShape === "triangle") { - // A triangle may have more than one side on the same boundary, so all of its sides are processed - this.boundaryElements[boundaryKey] - .filter(([boundaryElementIndex]) => boundaryElementIndex === elementIndex) - .forEach(([, side]) => { - const triangleContributions = this.computeTriangleConvectionContributions( - elementIndex, - side, - convectionHeatTransferCoefficient, - externalTemperature, - nodesXCoordinates, - nodesYCoordinates, - basisFunctions, - ); - for (const localNodeIndex of triangleContributions.sideNodeIndices) { - localResidualVector[localNodeIndex] += - triangleContributions.localResidualVector[localNodeIndex]; - for (const localNodeIndex2 of triangleContributions.sideNodeIndices) { - localJacobianMatrix[localNodeIndex][localNodeIndex2] += - triangleContributions.localJacobianMatrix[localNodeIndex][localNodeIndex2]; - } - } - }); - } else if (boundaryElement) { - const side = boundaryElement[1]; - - if (this.meshDimension === "1D") { - // Handle 1D case - let nodeIndex; - if (this.elementOrder === "linear") { - nodeIndex = side === 0 ? 0 : 1; - } else if (this.elementOrder === "quadratic") { - nodeIndex = side === 0 ? 0 : 2; - } - - // Add contribution to local Jacobian matrix and local residual vector - debugLog( - ` - Applied convection boundary condition to node ${nodeIndex + 1} (element ${ - elementIndex + 1 - }, local node ${nodeIndex + 1})`, - ); - localResidualVector[nodeIndex] += -convectionHeatTransferCoefficient * externalTemperature; - localJacobianMatrix[nodeIndex][nodeIndex] += convectionHeatTransferCoefficient; - } else if (this.meshDimension === "2D") { - // Handle 2D case - if (this.elementOrder === "linear") { - let gaussPoint1, gaussPoint2, firstNodeIndex, lastNodeIndex, nodeIncrement; - - if (side === 0) { - // Nodes at the bottom side of the reference element - gaussPoint1 = gaussPoints[0]; - gaussPoint2 = 0; - firstNodeIndex = 0; - lastNodeIndex = 3; - nodeIncrement = 2; - } else if (side === 1) { - // Nodes at the left side of the reference element - gaussPoint1 = 0; - gaussPoint2 = gaussPoints[0]; - firstNodeIndex = 0; - lastNodeIndex = 2; - nodeIncrement = 1; - } else if (side === 2) { - // Nodes at the top side of the reference element - gaussPoint1 = gaussPoints[0]; - gaussPoint2 = 1; - firstNodeIndex = 1; - lastNodeIndex = 4; - nodeIncrement = 2; - } else if (side === 3) { - // Nodes at the right side of the reference element - gaussPoint1 = 1; - gaussPoint2 = gaussPoints[0]; - firstNodeIndex = 2; - lastNodeIndex = 4; - nodeIncrement = 1; - } - - // Get basis functions - const basisFunctionsAndDerivatives = basisFunctions.getBasisFunctions(gaussPoint1, gaussPoint2); - const basisFunction = basisFunctionsAndDerivatives.basisFunction; - const basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi; - const basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta; - - // Calculate tangent vector components - let ksiDerivX = 0, - ksiDerivY = 0, - etaDerivX = 0, - etaDerivY = 0; - for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - - if (side === 0 || side === 2) { - ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - } else if (side === 1 || side === 3) { - etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - } - } - - // Compute tangent vector length - let tangentVectorLength; - if (side === 0 || side === 2) { - tangentVectorLength = Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2); - } else { - tangentVectorLength = Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2); - } - - // Apply boundary conditions to local matrices - for ( - let localNodeIndex = firstNodeIndex; - localNodeIndex < lastNodeIndex; - localNodeIndex += nodeIncrement - ) { - localResidualVector[localNodeIndex] += - -gaussWeights[0] * - tangentVectorLength * - basisFunction[localNodeIndex] * - convectionHeatTransferCoefficient * - externalTemperature; - - for ( - let localNodeIndex2 = firstNodeIndex; - localNodeIndex2 < lastNodeIndex; - localNodeIndex2 += nodeIncrement - ) { - localJacobianMatrix[localNodeIndex][localNodeIndex2] += - -gaussWeights[0] * - tangentVectorLength * - basisFunction[localNodeIndex] * - basisFunction[localNodeIndex2] * - convectionHeatTransferCoefficient; - } - } - } else if (this.elementOrder === "quadratic") { - // Handle quadratic elements (similar pattern but with more Gauss points) - for (let gaussPointIndex = 0; gaussPointIndex < 3; gaussPointIndex++) { - let gaussPoint1, gaussPoint2, firstNodeIndex, lastNodeIndex, nodeIncrement; - - if (side === 0) { - // Nodes at the bottom side of the reference element - gaussPoint1 = gaussPoints[gaussPointIndex]; - gaussPoint2 = 0; - firstNodeIndex = 0; - lastNodeIndex = 7; - nodeIncrement = 3; - } else if (side === 1) { - // Nodes at the left side of the reference element - gaussPoint1 = 0; - gaussPoint2 = gaussPoints[gaussPointIndex]; - firstNodeIndex = 0; - lastNodeIndex = 3; - nodeIncrement = 1; - } else if (side === 2) { - // Nodes at the top side of the reference element - gaussPoint1 = gaussPoints[gaussPointIndex]; - gaussPoint2 = 1; - firstNodeIndex = 2; - lastNodeIndex = 9; - nodeIncrement = 3; - } else if (side === 3) { - // Nodes at the right side of the reference element - gaussPoint1 = 1; - gaussPoint2 = gaussPoints[gaussPointIndex]; - firstNodeIndex = 6; - lastNodeIndex = 9; - nodeIncrement = 1; - } - let basisFunctionsAndDerivatives = basisFunctions.getBasisFunctions(gaussPoint1, gaussPoint2); - let basisFunction = basisFunctionsAndDerivatives.basisFunction; - let basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi; - let basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta; - - let ksiDerivX = 0; - let ksiDerivY = 0; - let etaDerivX = 0; - let etaDerivY = 0; - const nodesPerElement = this.nop[elementIndex].length; - for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { - const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; - - // For boundaries along Ksi (horizontal), use Ksi derivatives - if (side === 0 || side === 2) { - ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex]; - } - // For boundaries along Eta (vertical), use Eta derivatives - else if (side === 1 || side === 3) { - etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex]; - } - } - - // Compute the length of tangent vector - let tangentVectorLength; - if (side === 0 || side === 2) { - tangentVectorLength = Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2); - } else { - tangentVectorLength = Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2); - } - - // Apply boundary conditions to local matrices - for ( - let localNodeIndex = firstNodeIndex; - localNodeIndex < lastNodeIndex; - localNodeIndex += nodeIncrement - ) { - localResidualVector[localNodeIndex] += - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - convectionHeatTransferCoefficient * - externalTemperature; - - for ( - let localNodeIndex2 = firstNodeIndex; - localNodeIndex2 < lastNodeIndex; - localNodeIndex2 += nodeIncrement - ) { - localJacobianMatrix[localNodeIndex][localNodeIndex2] += - -gaussWeights[gaussPointIndex] * - tangentVectorLength * - basisFunction[localNodeIndex] * - basisFunction[localNodeIndex2] * - convectionHeatTransferCoefficient; - } - } - } - } - } - } - } - } - - return { localJacobianMatrix, localResidualVector }; - } -} diff --git a/src/models/beamBoundaryConditions.js b/src/physics/beamBoundaryConditions.js similarity index 100% rename from src/models/beamBoundaryConditions.js rename to src/physics/beamBoundaryConditions.js diff --git a/src/models/creepingFlow.js b/src/physics/creepingFlow.js similarity index 98% rename from src/models/creepingFlow.js rename to src/physics/creepingFlow.js index 8e506c4..7f34bd3 100644 --- a/src/models/creepingFlow.js +++ b/src/physics/creepingFlow.js @@ -8,10 +8,10 @@ */ // Internal imports -import { getIntegrationPoints2D, performIsoparametricMapping2D } from "../mesh/meshUtils.js"; -import { BasisFunctions } from "../mesh/basisFunctions.js"; -import { NumericalIntegration } from "../numerics/numericalIntegration.js"; -import { SparseMatrix } from "../numerics/sparseMatrix.js"; +import { getIntegrationPoints2D, performIsoparametricMapping2D } from "../elements/isoparametricMapping.js"; +import { BasisFunctions } from "../elements/basisFunctions.js"; +import { NumericalIntegration } from "../elements/numericalIntegration.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; import { FlowBoundaryConditions } from "./flowBoundaryConditions.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; diff --git a/src/models/eulerBernoulliBeam.js b/src/physics/eulerBernoulliBeam.js similarity index 97% rename from src/models/eulerBernoulliBeam.js rename to src/physics/eulerBernoulliBeam.js index e74bcfa..f68de4e 100644 --- a/src/models/eulerBernoulliBeam.js +++ b/src/physics/eulerBernoulliBeam.js @@ -8,9 +8,9 @@ */ // Internal imports -import { BasisFunctions } from "../mesh/basisFunctions.js"; -import { NumericalIntegration } from "../numerics/numericalIntegration.js"; -import { SparseMatrix } from "../numerics/sparseMatrix.js"; +import { BasisFunctions } from "../elements/basisFunctions.js"; +import { NumericalIntegration } from "../elements/numericalIntegration.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; import { BeamBoundaryConditions } from "./beamBoundaryConditions.js"; import { basicLog, errorLog } from "../utilities/logging.js"; diff --git a/src/models/flowBoundaryConditions.js b/src/physics/flowBoundaryConditions.js similarity index 99% rename from src/models/flowBoundaryConditions.js rename to src/physics/flowBoundaryConditions.js index aa08301..940fd32 100644 --- a/src/models/flowBoundaryConditions.js +++ b/src/physics/flowBoundaryConditions.js @@ -8,7 +8,7 @@ */ // Internal imports -import { triangleBoundarySides } from "../mesh/basisFunctions.js"; +import { triangleBoundarySides } from "../elements/basisFunctions.js"; import { debugLog } from "../utilities/logging.js"; /** diff --git a/src/models/frontPropagation.js b/src/physics/frontPropagation.js similarity index 99% rename from src/models/frontPropagation.js rename to src/physics/frontPropagation.js index 4573e2d..58f92e5 100644 --- a/src/models/frontPropagation.js +++ b/src/physics/frontPropagation.js @@ -9,7 +9,8 @@ // Internal imports import { GenericBoundaryConditions } from "./genericBoundaryConditions.js"; -import { initializeFEA, performIsoparametricMapping2D } from "../mesh/meshUtils.js"; +import { initializeFEA } from "../elements/feaInitialization.js"; +import { performIsoparametricMapping2D } from "../elements/isoparametricMapping.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; // Base viscous term that remains when eikonal equation is fully activated diff --git a/src/models/generalFormPDE.js b/src/physics/generalFormPDE.js similarity index 99% rename from src/models/generalFormPDE.js rename to src/physics/generalFormPDE.js index a9e4ab5..0a22d78 100644 --- a/src/models/generalFormPDE.js +++ b/src/physics/generalFormPDE.js @@ -8,7 +8,8 @@ */ // Internal imports -import { initializeFEA, performIsoparametricMapping1D } from "../mesh/meshUtils.js"; +import { initializeFEA } from "../elements/feaInitialization.js"; +import { performIsoparametricMapping1D } from "../elements/isoparametricMapping.js"; import { GenericBoundaryConditions } from "./genericBoundaryConditions.js"; import { basicLog, errorLog } from "../utilities/logging.js"; diff --git a/src/physics/genericBoundaryConditions.js b/src/physics/genericBoundaryConditions.js new file mode 100644 index 0000000..f378cbd --- /dev/null +++ b/src/physics/genericBoundaryConditions.js @@ -0,0 +1,99 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { getBoundarySides } from "../elements/basisFunctions.js"; +import { debugLog } from "../utilities/logging.js"; + +/** + * Class to handle generic boundary conditions application + */ +export class GenericBoundaryConditions { + /** + * Constructor to initialize the GenericBoundaryConditions class + * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis + * @param {array} boundaryElements - Array containing elements that belong to each boundary + * @param {array} nop - Nodal numbering (NOP) array representing the connectivity between elements and nodes + * @param {string} meshDimension - The dimension of the mesh (e.g., "2D") + * @param {string} elementOrder - The order of elements (e.g., "linear", "quadratic") + * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' + */ + constructor(boundaryConditions, boundaryElements, nop, meshDimension, elementOrder, elementShape = "quad") { + this.boundaryConditions = boundaryConditions; + this.boundaryElements = boundaryElements; + this.nop = nop; + this.meshDimension = meshDimension; + this.elementOrder = elementOrder; + this.elementShape = elementShape; + } + + /** + * Function to visit every node on a constant value boundary + * @param {function} callback - Called with (globalNodeIndex, value) for each node + */ + forEachConstantValueNode(callback) { + const boundarySides = getBoundarySides(this.meshDimension, this.elementOrder, this.elementShape); + Object.keys(this.boundaryConditions).forEach((boundaryKey) => { + if (this.boundaryConditions[boundaryKey][0] === "constantValue") { + const value = this.boundaryConditions[boundaryKey][1]; + debugLog(`Boundary ${boundaryKey}: Applying constant value of ${value} (Dirichlet condition)`); + + this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { + boundarySides[side].forEach((nodeIndex) => { + const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; + debugLog( + ` - Applied constant value to node ${globalNodeIndex + 1} (element ${ + elementIndex + 1 + }, local node ${nodeIndex + 1})`, + ); + callback(globalNodeIndex, value); + }); + }); + } + }); + } + + /** + * Function to impose Dirichlet boundary conditions + * @param {array} residualVector - The residual vector to be modified + * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified + * @param {array} [solutionVector] - Current solution (Newton-Raphson iterate); when provided, the + * residual is set to the increment needed to reach the prescribed value instead of the value itself + * + * For consistency across both linear and nonlinear formulations, + * this project always refers 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. + */ + imposeDirichletBoundaryConditions(residualVector, jacobianMatrix, solutionVector) { + this.forEachConstantValueNode((globalNodeIndex, value) => { + // Set the residual vector to the value, or the increment needed to reach it for Newton-Raphson + residualVector[globalNodeIndex] = solutionVector ? value - solutionVector[globalNodeIndex] : value; + // Set the Jacobian matrix row to zero + jacobianMatrix.clearRow(globalNodeIndex); + // Set the diagonal entry of the Jacobian matrix to one + jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); + }); + } + + /** + * Function to impose constant value (Dirichlet) boundary conditions for the frontal solver + * @param {array} nodeConstraintCode - Array indicating boundary condition code for each node + * @param {array} boundaryValues - Array containing boundary condition values + */ + imposeDirichletBoundaryConditionsFront(nodeConstraintCode, boundaryValues) { + this.forEachConstantValueNode((globalNodeIndex, value) => { + nodeConstraintCode[globalNodeIndex] = 1; + boundaryValues[globalNodeIndex] = value; + }); + } +} diff --git a/src/models/heatConduction.js b/src/physics/heatConduction.js similarity index 98% rename from src/models/heatConduction.js rename to src/physics/heatConduction.js index a8a2734..6d0bd47 100644 --- a/src/models/heatConduction.js +++ b/src/physics/heatConduction.js @@ -8,11 +8,8 @@ */ // Internal imports -import { - initializeFEA, - performIsoparametricMapping1D, - performIsoparametricMapping2D, -} from "../mesh/meshUtils.js"; +import { initializeFEA } from "../elements/feaInitialization.js"; +import { performIsoparametricMapping1D, performIsoparametricMapping2D } from "../elements/isoparametricMapping.js"; import { ThermalBoundaryConditions } from "./thermalBoundaryConditions.js"; import { basicLog } from "../utilities/logging.js"; diff --git a/src/physics/thermalBoundaryConditions.js b/src/physics/thermalBoundaryConditions.js new file mode 100644 index 0000000..d75016c --- /dev/null +++ b/src/physics/thermalBoundaryConditions.js @@ -0,0 +1,425 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +// Internal imports +import { triangleBoundarySides, getBoundarySides } from "../elements/basisFunctions.js"; +import { NumericalIntegration } from "../elements/numericalIntegration.js"; +import { debugLog } from "../utilities/logging.js"; + +/** + * Class to handle thermal boundary conditions application + */ +export class ThermalBoundaryConditions { + /** + * Constructor to initialize the ThermalBoundaryConditions class + * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis + * @param {array} boundaryElements - Array containing elements that belong to each boundary + * @param {array} nop - Nodal numbering (NOP) array representing the connectivity between elements and nodes + * @param {string} meshDimension - The dimension of the mesh (e.g., "2D") + * @param {string} elementOrder - The order of elements (e.g., "linear", "quadratic") + * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle' + */ + constructor(boundaryConditions, boundaryElements, nop, meshDimension, elementOrder, elementShape = "quad") { + this.boundaryConditions = boundaryConditions; + this.boundaryElements = boundaryElements; + this.nop = nop; + this.meshDimension = meshDimension; + this.elementOrder = elementOrder; + this.elementShape = elementShape; + } + + /** + * Function to compute the convection contributions on one side of a triangular element + * @param {number} elementIndex - Index of the element + * @param {number} side - Side of the reference triangle (0: bottom, 1: left, 2: hypotenuse) + * @param {number} convectionHeatTransferCoefficient - Convection heat transfer coefficient + * @param {number} externalTemperature - External temperature + * @param {array} nodesXCoordinates - Array of x-coordinates of nodes + * @param {array} nodesYCoordinates - Array of y-coordinates of nodes + * @param {object} basisFunctions - Object containing basis functions and their derivatives + * @returns {object} An object containing: + * - sideNodeIndices: Local node indices on the side + * - localJacobianMatrix: Contributions indexed by local node indices + * - localResidualVector: Contributions indexed by local node indices + */ + computeTriangleConvectionContributions( + elementIndex, + side, + convectionHeatTransferCoefficient, + externalTemperature, + nodesXCoordinates, + nodesYCoordinates, + basisFunctions, + ) { + const nodesPerElement = this.nop[elementIndex].length; + const sideNodeIndices = triangleBoundarySides[this.elementOrder][side]; + const localJacobianMatrix = Array(nodesPerElement) + .fill() + .map(() => Array(nodesPerElement).fill(0)); + const localResidualVector = Array(nodesPerElement).fill(0); + + // Parametrize the side with t in [0, 1]: (ksi, eta) = start + t * direction + const sideParametrization = { + 0: { start: [0, 0], direction: [1, 0] }, // Bottom side (eta = 0) + 1: { start: [0, 0], direction: [0, 1] }, // Left side (ksi = 0) + 2: { start: [0, 1], direction: [1, -1] }, // Hypotenuse (ksi + eta = 1) + }; + const { start, direction } = sideParametrization[side]; + + // 3-point Gauss rule on [0, 1], exact for the boundary mass terms of both linear and quadratic triangles + const { gaussPoints, gaussWeights } = new NumericalIntegration({ + meshDimension: "1D", + elementOrder: "quadratic", + }).getGaussPointsAndWeights(); + + for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { + const t = gaussPoints[gaussPointIndex]; + const { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta } = + basisFunctions.getBasisFunctions(start[0] + t * direction[0], start[1] + t * direction[1]); + + // Tangent vector of the side in physical space (chain rule along the side parametrization) + let tangentX = 0; + let tangentY = 0; + for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { + const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; + const basisFunctionDerivT = + basisFunctionDerivKsi[nodeIndex] * direction[0] + basisFunctionDerivEta[nodeIndex] * direction[1]; + tangentX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivT; + tangentY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivT; + } + const tangentVectorLength = Math.sqrt(tangentX ** 2 + tangentY ** 2); + + for (const localNodeIndex of sideNodeIndices) { + localResidualVector[localNodeIndex] += + -gaussWeights[gaussPointIndex] * + tangentVectorLength * + basisFunction[localNodeIndex] * + convectionHeatTransferCoefficient * + externalTemperature; + + for (const localNodeIndex2 of sideNodeIndices) { + localJacobianMatrix[localNodeIndex][localNodeIndex2] += + -gaussWeights[gaussPointIndex] * + tangentVectorLength * + basisFunction[localNodeIndex] * + basisFunction[localNodeIndex2] * + convectionHeatTransferCoefficient; + } + } + } + + return { sideNodeIndices, localJacobianMatrix, localResidualVector }; + } + + /** + * Function to compute the convection contributions on one side of an element (1D, quadrilateral or triangle) + * @param {number} elementIndex - Index of the element + * @param {number} side - Side of the reference element + * @param {number} convectionHeatTransferCoefficient - Convection heat transfer coefficient + * @param {number} externalTemperature - External temperature + * @param {array} nodesXCoordinates - Array of x-coordinates of nodes + * @param {array} nodesYCoordinates - Array of y-coordinates of nodes + * @param {array} gaussPoints - Array of Gauss points for numerical integration + * @param {array} gaussWeights - Array of Gauss weights for numerical integration + * @param {object} basisFunctions - Object containing basis functions and their derivatives + * @returns {object} An object containing: + * - sideNodeIndices: Local node indices on the side + * - localJacobianMatrix: Contributions indexed by local node indices + * - localResidualVector: Contributions indexed by local node indices + */ + computeConvectionContributions( + elementIndex, + side, + convectionHeatTransferCoefficient, + externalTemperature, + nodesXCoordinates, + nodesYCoordinates, + gaussPoints, + gaussWeights, + basisFunctions, + ) { + if (this.meshDimension === "2D" && this.elementShape === "triangle") { + return this.computeTriangleConvectionContributions( + elementIndex, + side, + convectionHeatTransferCoefficient, + externalTemperature, + nodesXCoordinates, + nodesYCoordinates, + basisFunctions, + ); + } + + const nodesPerElement = this.nop[elementIndex].length; + const sideNodeIndices = getBoundarySides(this.meshDimension, this.elementOrder, this.elementShape)[side]; + const localJacobianMatrix = Array(nodesPerElement) + .fill() + .map(() => Array(nodesPerElement).fill(0)); + const localResidualVector = Array(nodesPerElement).fill(0); + + if (this.meshDimension === "1D") { + // The boundary is a single node + const nodeIndex = sideNodeIndices[0]; + localResidualVector[nodeIndex] += -convectionHeatTransferCoefficient * externalTemperature; + localJacobianMatrix[nodeIndex][nodeIndex] += convectionHeatTransferCoefficient; + return { sideNodeIndices, localJacobianMatrix, localResidualVector }; + } + + // Quadrilateral: integrate along the side with the 1D Gauss rule + for (let gaussPointIndex = 0; gaussPointIndex < gaussPoints.length; gaussPointIndex++) { + // Reference coordinates of the Gauss point on the side + const gaussPoint = gaussPoints[gaussPointIndex]; + const sideCoordinates = { + 0: [gaussPoint, 0], // Bottom side (eta = 0) + 1: [0, gaussPoint], // Left side (ksi = 0) + 2: [gaussPoint, 1], // Top side (eta = 1) + 3: [1, gaussPoint], // Right side (ksi = 1) + }; + const [ksi, eta] = sideCoordinates[side]; + const { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta } = + basisFunctions.getBasisFunctions(ksi, eta); + + // Tangent vector of the side: along ksi for the bottom/top sides, along eta for the left/right sides + const isKsiSide = side === 0 || side === 2; + let tangentX = 0; + let tangentY = 0; + for (let nodeIndex = 0; nodeIndex < nodesPerElement; nodeIndex++) { + const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; + const basisFunctionDerivT = isKsiSide + ? basisFunctionDerivKsi[nodeIndex] + : basisFunctionDerivEta[nodeIndex]; + tangentX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivT; + tangentY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivT; + } + const tangentVectorLength = Math.sqrt(tangentX ** 2 + tangentY ** 2); + + for (const localNodeIndex of sideNodeIndices) { + localResidualVector[localNodeIndex] += + -gaussWeights[gaussPointIndex] * + tangentVectorLength * + basisFunction[localNodeIndex] * + convectionHeatTransferCoefficient * + externalTemperature; + + for (const localNodeIndex2 of sideNodeIndices) { + localJacobianMatrix[localNodeIndex][localNodeIndex2] += + -gaussWeights[gaussPointIndex] * + tangentVectorLength * + basisFunction[localNodeIndex] * + basisFunction[localNodeIndex2] * + convectionHeatTransferCoefficient; + } + } + } + + return { sideNodeIndices, localJacobianMatrix, localResidualVector }; + } + + /** + * Function to visit every node on a constant temperature boundary + * @param {function} callback - Called with (globalNodeIndex, constantTemperatureValue) for each node + */ + forEachConstantTemperatureNode(callback) { + const boundarySides = getBoundarySides(this.meshDimension, this.elementOrder, this.elementShape); + Object.keys(this.boundaryConditions).forEach((boundaryKey) => { + if (this.boundaryConditions[boundaryKey][0] === "constantTemperature") { + const constantTemperatureValue = this.boundaryConditions[boundaryKey][1]; + debugLog( + `Boundary ${boundaryKey}: Applying constant temperature of ${constantTemperatureValue} K (Dirichlet condition)`, + ); + + this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { + boundarySides[side].forEach((nodeIndex) => { + const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1; + debugLog( + ` - Applied constant temperature to node ${globalNodeIndex + 1} (element ${ + elementIndex + 1 + }, local node ${nodeIndex + 1})`, + ); + callback(globalNodeIndex, constantTemperatureValue); + }); + }); + } + }); + } + + /** + * Function to impose constant temperature boundary conditions (Dirichlet type) + * @param {array} residualVector - The residual vector to be modified + * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified + * + * For consistency across both linear and nonlinear formulations, + * this project always refers 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. + */ + imposeConstantTempBoundaryConditions(residualVector, jacobianMatrix) { + this.forEachConstantTemperatureNode((globalNodeIndex, constantTemperatureValue) => { + // Set the residual vector to the constant temperature value + residualVector[globalNodeIndex] = constantTemperatureValue; + // Set the Jacobian matrix row to zero + jacobianMatrix.clearRow(globalNodeIndex); + // Set the diagonal entry of the Jacobian matrix to one + jacobianMatrix.setValue(globalNodeIndex, globalNodeIndex, 1); + }); + } + + /** + * Function to impose constant temperature boundary conditions for the frontal solver + * @param {array} nodeConstraintCode - Array indicating boundary condition code for each node + * @param {array} boundaryValues - Array containing boundary condition values + */ + imposeDirichletBoundaryConditionsFront(nodeConstraintCode, boundaryValues) { + this.forEachConstantTemperatureNode((globalNodeIndex, constantTemperatureValue) => { + // Set boundary condition code and value + nodeConstraintCode[globalNodeIndex] = 1; + boundaryValues[globalNodeIndex] = constantTemperatureValue; + }); + } + + /** + * Function to impose convection boundary conditions (Robin type) + * @param {array} residualVector - The residual vector to be modified + * @param {SparseMatrix} jacobianMatrix - The Jacobian matrix to be modified + * @param {array} gaussPoints - Array of Gauss points for numerical integration + * @param {array} gaussWeights - Array of Gauss weights for numerical integration + * @param {array} nodesXCoordinates - Array of x-coordinates of nodes + * @param {array} nodesYCoordinates - Array of y-coordinates of nodes + * @param {object} basisFunctions - Object containing basis functions and their derivatives + */ + imposeConvectionBoundaryConditions( + residualVector, + jacobianMatrix, + gaussPoints, + gaussWeights, + nodesXCoordinates, + nodesYCoordinates, + basisFunctions, + ) { + Object.keys(this.boundaryConditions).forEach((boundaryKey) => { + if (this.boundaryConditions[boundaryKey][0] === "convection") { + const [, convectionHeatTransferCoefficient, externalTemperature] = + this.boundaryConditions[boundaryKey]; + debugLog( + `Boundary ${boundaryKey}: Applying convection with heat transfer coefficient h=${convectionHeatTransferCoefficient} W/(m²·K) and external temperature T∞=${externalTemperature} K`, + ); + + this.boundaryElements[boundaryKey].forEach(([elementIndex, side]) => { + const { sideNodeIndices, localJacobianMatrix, localResidualVector } = + this.computeConvectionContributions( + elementIndex, + side, + convectionHeatTransferCoefficient, + externalTemperature, + nodesXCoordinates, + nodesYCoordinates, + gaussPoints, + gaussWeights, + basisFunctions, + ); + + // Add the side contributions to the global system + for (const localNodeIndex of sideNodeIndices) { + const globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1; + debugLog( + ` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${ + elementIndex + 1 + }, local node ${localNodeIndex + 1})`, + ); + residualVector[globalNodeIndex] += localResidualVector[localNodeIndex]; + + for (const localNodeIndex2 of sideNodeIndices) { + const globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1; + jacobianMatrix.addValue( + globalNodeIndex, + globalNodeIndex2, + localJacobianMatrix[localNodeIndex][localNodeIndex2], + ); + } + } + }); + } + }); + } + + /** + * Function to impose convection boundary conditions for the frontal solver + * @param {number} elementIndex - Index of the element being processed + * @param {array} nodesXCoordinates - Array of x-coordinates of nodes + * @param {array} nodesYCoordinates - Array of y-coordinates of nodes + * @param {array} gaussPoints - Array of Gauss points for numerical integration + * @param {array} gaussWeights - Array of Gauss weights for numerical integration + * @param {object} basisFunctions - Object containing basis functions and their derivatives + * @returns {object|null} An object containing the following, or null when the element is not on a convection boundary: + * - localJacobianMatrix: Local Jacobian matrix with convection contributions + * - localResidualVector: Residual vector with convection contributions + */ + imposeRobinBoundaryConditionsFront( + elementIndex, + nodesXCoordinates, + nodesYCoordinates, + gaussPoints, + gaussWeights, + basisFunctions, + ) { + let localJacobianMatrix = null; + let localResidualVector = null; + + for (const boundaryKey in this.boundaryElements) { + if (this.boundaryConditions[boundaryKey]?.[0] !== "convection") continue; + const [, convectionHeatTransferCoefficient, externalTemperature] = this.boundaryConditions[boundaryKey]; + + // An element may have more than one side on the same boundary, so all of its sides are processed + this.boundaryElements[boundaryKey] + .filter(([boundaryElementIndex]) => boundaryElementIndex === elementIndex) + .forEach(([, side]) => { + debugLog( + `Boundary ${boundaryKey}: Applying convection to element ${ + elementIndex + 1 + }, side ${side} with heat transfer coefficient h=${convectionHeatTransferCoefficient} W/(m²·K) and external temperature T∞=${externalTemperature} K`, + ); + + // Initialize local Jacobian matrix and local residual vector on the first convection side + if (!localJacobianMatrix) { + const nodesPerElement = this.nop[elementIndex].length; + localJacobianMatrix = Array(nodesPerElement) + .fill() + .map(() => Array(nodesPerElement).fill(0)); + localResidualVector = Array(nodesPerElement).fill(0); + } + + const sideContributions = this.computeConvectionContributions( + elementIndex, + side, + convectionHeatTransferCoefficient, + externalTemperature, + nodesXCoordinates, + nodesYCoordinates, + gaussPoints, + gaussWeights, + basisFunctions, + ); + for (const localNodeIndex of sideContributions.sideNodeIndices) { + localResidualVector[localNodeIndex] += sideContributions.localResidualVector[localNodeIndex]; + for (const localNodeIndex2 of sideContributions.sideNodeIndices) { + localJacobianMatrix[localNodeIndex][localNodeIndex2] += + sideContributions.localJacobianMatrix[localNodeIndex][localNodeIndex2]; + } + } + }); + } + + // Elements that are not on a convection boundary have no Robin-type contributions + if (!localJacobianMatrix) return null; + return { localJacobianMatrix, localResidualVector }; + } +} diff --git a/src/readers/gmshReader.js b/src/readers/gmshReader.js deleted file mode 100644 index 5f1a792..0000000 --- a/src/readers/gmshReader.js +++ /dev/null @@ -1,296 +0,0 @@ -/** - * ════════════════════════════════════════════════════════════════ - * FEAScript Core Library - * Finite Element Simulation in JavaScript - * Version: 0.3.0 (RC) | https://feascript.com - * MIT License © 2023–2026 FEAScript - * ════════════════════════════════════════════════════════════════ - */ - -// Internal imports -import { debugLog } from "../utilities/logging.js"; - -/** - * Function to import mesh data from Gmsh (.msh v4.1) containing quadrilateral and triangular elements - * @param {File} file - * @returns {object} - */ -const importGmshMesh = async (file) => { - let result = { - nodesXCoordinates: [], - nodesYCoordinates: [], - nodalNumbering: { - quadElements: [], - triangleElements: [], - }, - boundaryElements: [], - boundaryConditions: [], - boundaryNodePairs: {}, - gmshV: 0, - ascii: false, - fltBytes: "8", - totalNodesX: 0, - totalNodesY: 0, - physicalPropMap: [], - elementTypes: {}, - }; - - // Entities to physical tags map - const entityPhysicalMap = { curves: {} }; - - let content = await file.text(); - let lines = content - .split("\n") - .map((l) => l.trim()) - .filter((l) => l !== ""); - - let section = ""; - let lineIndex = 0; - - // Nodes - let nodeEntityBlocks = 0; - let totalNodes = 0; - let nodeBlocksProcessed = 0; - let currentNodeBlock = { numNodes: 0 }; - let nodeTags = []; - 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; - let currentElementBlock = { numElements: 0 }; - let elementsProcessedInBlock = 0; - - let boundaryElementsByTag = {}; - - // Entities helpers - let entityCounts = null; - let entitiesPhase = null; - let processedPoints = 0; - let processedCurves = 0; - let processedSurfaces = 0; - - while (lineIndex < lines.length) { - const line = lines[lineIndex]; - - // Section switches - if (line === "$MeshFormat") { - section = "meshFormat"; - lineIndex++; - continue; - } - if (line === "$EndMeshFormat") { - section = ""; - lineIndex++; - continue; - } - if (line === "$PhysicalNames") { - section = "physicalNames"; - lineIndex++; - continue; - } - if (line === "$EndPhysicalNames") { - section = ""; - lineIndex++; - continue; - } - if (line === "$Entities") { - section = "entities"; - entitiesPhase = "counts"; - lineIndex++; - continue; - } - if (line === "$EndEntities") { - section = ""; - entityCounts = null; - entitiesPhase = null; - lineIndex++; - continue; - } - if (line === "$Nodes") { - section = "nodes"; - lineIndex++; - continue; - } - if (line === "$EndNodes") { - section = ""; - lineIndex++; - continue; - } - if (line === "$Elements") { - section = "elements"; - lineIndex++; - continue; - } - if (line === "$EndElements") { - section = ""; - lineIndex++; - continue; - } - - const parts = line.split(/\s+/); - - // Mesh format - if (section === "meshFormat") { - result.gmshV = parseFloat(parts[0]); - result.ascii = parts[1] === "0"; - result.fltBytes = parts[2]; - } - - // Physical names - else if (section === "physicalNames") { - const dimension = parseInt(parts[0], 10); - const tag = parseInt(parts[1], 10); - let name = parts.slice(2).join(" ").replace(/^"|"$/g, ""); - - result.physicalPropMap.push({ tag, dimension, name }); - } - - // Entities - else if (section === "entities") { - if (entitiesPhase === "counts") { - entityCounts = { - points: parseInt(parts[0], 10), - curves: parseInt(parts[1], 10), - surfaces: parseInt(parts[2], 10), - volumes: parseInt(parts[3], 10), - }; - entitiesPhase = "points"; - } else if (entitiesPhase === "points") { - processedPoints++; - if (processedPoints === entityCounts.points) entitiesPhase = "curves"; - } else if (entitiesPhase === "curves") { - const tag = parseInt(parts[0], 10); - const numPhysical = parseInt(parts[7], 10); - const physTags = parts.slice(8, 8 + numPhysical).map((p) => parseInt(p, 10)); - entityPhysicalMap.curves[tag] = physTags; - processedCurves++; - if (processedCurves === entityCounts.curves) entitiesPhase = "surfaces"; - } else if (entitiesPhase === "surfaces") { - processedSurfaces++; - if (processedSurfaces === entityCounts.surfaces) entitiesPhase = "volumes"; - } - } - - // Nodes - else if (section === "nodes") { - 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) { - currentNodeBlock = { - dim: parseInt(parts[0], 10), - tag: parseInt(parts[1], 10), - parametric: parseInt(parts[2], 10), - numNodes: parseInt(parts[3], 10), - }; - nodeTags = []; - nodeTagsCollected = 0; - nodeCoordinatesCollected = 0; - } else if (nodeTagsCollected < currentNodeBlock.numNodes) { - for (let p of parts) { - nodeTags.push(parseInt(p, 10)); - nodeTagsCollected++; - if (nodeTagsCollected === currentNodeBlock.numNodes) break; - } - } else if (nodeCoordinatesCollected < currentNodeBlock.numNodes) { - 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++; - if (nodeCoordinatesCollected === currentNodeBlock.numNodes) { - nodeBlocksProcessed++; - currentNodeBlock = { numNodes: 0 }; - } - } - } - - // Elements - else if (section === "elements") { - if (elementEntityBlocks === 0) { - elementEntityBlocks = parseInt(parts[0], 10); - } else if (elementBlocksProcessed < elementEntityBlocks && currentElementBlock.numElements === 0) { - currentElementBlock = { - dim: parseInt(parts[0], 10), - tag: parseInt(parts[1], 10), // entity tag - elementType: parseInt(parts[2], 10), - numElements: parseInt(parts[3], 10), - }; - - result.elementTypes[currentElementBlock.elementType] = - (result.elementTypes[currentElementBlock.elementType] || 0) + currentElementBlock.numElements; - - elementsProcessedInBlock = 0; - } else if (elementsProcessedInBlock < currentElementBlock.numElements) { - 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; - if (currentElementBlock.dim === 1) { - const mapped = entityPhysicalMap.curves[currentElementBlock.tag]; - if (mapped && mapped.length > 0) physicalTag = mapped[0]; - } - - if (currentElementBlock.elementType === 1 || currentElementBlock.elementType === 8) { - if (!boundaryElementsByTag[physicalTag]) boundaryElementsByTag[physicalTag] = []; - boundaryElementsByTag[physicalTag].push(nodeIndices); - - if (!result.boundaryNodePairs[physicalTag]) result.boundaryNodePairs[physicalTag] = []; - result.boundaryNodePairs[physicalTag].push(nodeIndices); - } else if (currentElementBlock.elementType === 2 || currentElementBlock.elementType === 9) { - result.nodalNumbering.triangleElements.push(nodeIndices); - } else if (currentElementBlock.elementType === 3 || currentElementBlock.elementType === 10) { - result.nodalNumbering.quadElements.push(nodeIndices); - } - - elementsProcessedInBlock++; - if (elementsProcessedInBlock === currentElementBlock.numElements) { - elementBlocksProcessed++; - currentElementBlock = { numElements: 0 }; - } - } - } - - lineIndex++; - } - - // Boundary conditions - result.physicalPropMap.forEach((prop) => { - if (prop.dimension === 1) { - const nodes = boundaryElementsByTag[prop.tag] || []; - if (nodes.length > 0) { - result.boundaryConditions.push({ - name: prop.name, - tag: prop.tag, - nodes, - }); - } - } - }); - - debugLog(`Parsed boundary node pairs by physical tag: ${JSON.stringify(result.boundaryNodePairs)}`); - - return result; -}; - -export { importGmshMesh }; diff --git a/src/solvers/frontalSolver.js b/src/solvers/frontalSolver.js index 7888553..8c38fe2 100644 --- a/src/solvers/frontalSolver.js +++ b/src/solvers/frontalSolver.js @@ -8,12 +8,8 @@ */ // Internal imports -import { BasisFunctions } from "../mesh/basisFunctions.js"; -import { initializeFEA } from "../mesh/meshUtils.js"; -import { assembleHeatConductionFront } from "../models/heatConduction.js"; -import { ThermalBoundaryConditions } from "../models/thermalBoundaryConditions.js"; -import { assembleFrontPropagationFront } from "../models/frontPropagation.js"; -import { GenericBoundaryConditions } from "../models/genericBoundaryConditions.js"; +import { BasisFunctions } from "../elements/basisFunctions.js"; +import { initializeFEA } from "../elements/feaInitialization.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; // Create object templates @@ -25,13 +21,23 @@ let basisFunctions; /** * Function to run the frontal solver and obtain results for plotting - * @param {function} assembleFront - Matrix assembler based on the physical model + * + * The solver is independent of the physical model: the model supplies its element assembler and its + * boundary conditions handler. Model-specific parameters (e.g. coefficient functions, the current + * Newton-Raphson iterate) are bound into `assembleFront` by the caller. + * + * @param {function} assembleFront - Element assembler of the physical model. It receives + * `{ elementIndex, nop, meshData, basisFunctions, FEAData }` and returns + * `{ localJacobianMatrix, localResidualVector, ngl }` * @param {object} meshData - Object containing mesh data - * @param {object} boundaryConditions - Object containing boundary conditions - * @param {object} [options] - Additional options for the solver + * @param {object} boundaryConditionsHandler - Boundary conditions object of the physical model, implementing: + * - imposeDirichletBoundaryConditionsFront(nodeConstraintCode, boundaryValues) + * - imposeRobinBoundaryConditionsFront(elementIndex, nodesXCoordinates, nodesYCoordinates, gaussPoints, + * gaussWeights, basisFunctions) (optional) - returns `{ localJacobianMatrix, localResidualVector }`, + * or null when the element is not on a Robin-type boundary * @returns {object} An object containing the solution vector (Float64Array) and node coordinates */ -export function runFrontalSolver(assembleFront, meshData, boundaryConditions, options = {}) { +export function runFrontalSolver(assembleFront, meshData, boundaryConditionsHandler) { // Initialize FEA components const FEAData = initializeFEA(meshData); const totalNodes = meshData.nodesXCoordinates.length; @@ -66,39 +72,11 @@ export function runFrontalSolver(assembleFront, meshData, boundaryConditions, op frontalData.boundaryValues[nodeIndex] = 0; } - // Handle Dirichlet-type boundary conditions differently based on which solver is being used - let dirichletBoundaryConditionsHandler; - // Solid heat transfer model (heatConduction solver) - if (assembleFront === assembleHeatConductionFront) { - dirichletBoundaryConditionsHandler = new ThermalBoundaryConditions( - boundaryConditions, - meshData.boundaryElements, - meshData.nop, - meshData.meshDimension, - meshData.elementOrder, - meshData.elementShape, - ); - - dirichletBoundaryConditionsHandler.imposeConstantTempBoundaryConditionsFront( - frontalData.nodeConstraintCode, - frontalData.boundaryValues, - ); - // Front propagation model (frontPropagation solver) - } else if (assembleFront === assembleFrontPropagationFront) { - dirichletBoundaryConditionsHandler = new GenericBoundaryConditions( - boundaryConditions, - meshData.boundaryElements, - meshData.nop, - meshData.meshDimension, - meshData.elementOrder, - meshData.elementShape, - ); - - dirichletBoundaryConditionsHandler.imposeConstantValueBoundaryConditionsFront( - frontalData.nodeConstraintCode, - frontalData.boundaryValues, - ); - } + boundaryConditionsHandler.imposeDirichletBoundaryConditionsFront( + frontalData.nodeConstraintCode, + frontalData.boundaryValues, + ); + // Initialize global residual vector for (let nodeIndex = 0; nodeIndex < meshData.nodesXCoordinates.length; nodeIndex++) { frontalData.globalResidualVector[nodeIndex] = 0; @@ -113,13 +91,7 @@ export function runFrontalSolver(assembleFront, meshData, boundaryConditions, op frontalState.nodesPerElement[elementIndex] = FEAData.nodesPerElement; } - // 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); + runFrontalAlgorithm(meshData, FEAData, boundaryConditionsHandler, assembleFront); // Copy solution for (let nodeIndex = 0; nodeIndex < meshData.nodesXCoordinates.length; nodeIndex++) { @@ -220,10 +192,10 @@ function estimateFrontSize(nodesPerElement, numElements) { * Function to compute local Jacobian matrix and local residual vector * @param {object} meshData - Object containing mesh data * @param {object} FEAData - Object containing FEA-related data - * @param {object} thermalBoundaryConditions - Object containing thermal boundary conditions - * @param {function} assembleFront - Matrix assembler based on the physical model + * @param {object} boundaryConditionsHandler - Boundary conditions object of the physical model + * @param {function} assembleFront - Element assembler of the physical model */ -function assembleElementContribution(meshData, FEAData, thermalBoundaryConditions, assembleFront) { +function assembleElementContribution(meshData, FEAData, boundaryConditionsHandler, assembleFront) { const elementIndex = elementData.currentElementIndex - 1; // Guard against out-of-range indices @@ -239,57 +211,24 @@ function assembleElementContribution(meshData, FEAData, thermalBoundaryCondition meshData, basisFunctions: basisFunctions, FEAData, - // 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 - let boundaryLocalJacobianMatrix = Array(FEAData.nodesPerElement) - .fill() - .map(() => Array(FEAData.nodesPerElement).fill(0)); - let boundaryResidualVector = Array(FEAData.nodesPerElement).fill(0); - - // heatConduction solver - if (assembleFront === assembleHeatConductionFront) { - // Check if this element is on a Robin-type boundary - let isOnRobinTypeBoundary = false; - for (const boundaryKey in meshData.boundaryElements) { - if ( - thermalBoundaryConditions.boundaryConditions[boundaryKey]?.[0] === "convection" && - meshData.boundaryElements[boundaryKey].some( - ([boundaryElementIndex]) => boundaryElementIndex === elementIndex, - ) - ) { - isOnRobinTypeBoundary = true; - break; - } - } - - // Only calculate Robin-type for elements when required - if (isOnRobinTypeBoundary) { - const { gaussPoints, gaussWeights } = FEAData; - const result = thermalBoundaryConditions.imposeConvectionBoundaryConditionsFront( - elementIndex, - meshData.nodesXCoordinates, - meshData.nodesYCoordinates, - gaussPoints, - gaussWeights, - basisFunctions, - ); - boundaryLocalJacobianMatrix = result.localJacobianMatrix; - boundaryResidualVector = result.localResidualVector; - } - } else if (assembleFront === assembleFrontPropagationFront) { - // For now, no Robin-type boundary conditions exist for any other solver - } + // Robin-type boundary terms (only for models whose handler supports them) + const robinContribution = boundaryConditionsHandler.imposeRobinBoundaryConditionsFront?.( + elementIndex, + meshData.nodesXCoordinates, + meshData.nodesYCoordinates, + FEAData.gaussPoints, + FEAData.gaussWeights, + basisFunctions, + ); // Combine domain and boundary contributions for (let localNodeI = 0; localNodeI < FEAData.nodesPerElement; localNodeI++) { for (let localNodeJ = 0; localNodeJ < FEAData.nodesPerElement; localNodeJ++) { elementData.localJacobianMatrix[localNodeI][localNodeJ] = - localJacobianMatrix[localNodeI][localNodeJ] + boundaryLocalJacobianMatrix[localNodeI][localNodeJ]; + localJacobianMatrix[localNodeI][localNodeJ] + + (robinContribution ? robinContribution.localJacobianMatrix[localNodeI][localNodeJ] : 0); } } @@ -297,7 +236,8 @@ function assembleElementContribution(meshData, FEAData, thermalBoundaryCondition for (let localNodeIndex = 0; localNodeIndex < FEAData.nodesPerElement; localNodeIndex++) { const globalNodeIndex = ngl[localNodeIndex] - 1; frontalData.globalResidualVector[globalNodeIndex] += - localResidualVector[localNodeIndex] + boundaryResidualVector[localNodeIndex]; + localResidualVector[localNodeIndex] + + (robinContribution ? robinContribution.localResidualVector[localNodeIndex] : 0); } return true; @@ -307,10 +247,10 @@ function assembleElementContribution(meshData, FEAData, thermalBoundaryCondition * Function to implement the frontal solver algorithm * @param {object} meshData - Object containing mesh data * @param {object} FEAData - Object containing FEA-related data - * @param {object} thermalBoundaryConditions - Object containing thermal boundary conditions - * @param {function} assembleFront - Matrix assembler based on the physical model + * @param {object} boundaryConditionsHandler - Boundary conditions object of the physical model + * @param {function} assembleFront - Element assembler of the physical model */ -function runFrontalAlgorithm(meshData, FEAData, thermalBoundaryConditions, assembleFront) { +function runFrontalAlgorithm(meshData, FEAData, boundaryConditionsHandler, assembleFront) { // Allocate local arrays dynamically const totalElements = meshData.totalElements; const nodesPerElement = meshData.nodesXCoordinates.length; @@ -382,7 +322,7 @@ function runFrontalAlgorithm(meshData, FEAData, thermalBoundaryConditions, assem if (elementData.currentElementIndex < totalElements) { elementData.currentElementIndex++; - assembled = assembleElementContribution(meshData, FEAData, thermalBoundaryConditions, assembleFront); + assembled = assembleElementContribution(meshData, FEAData, boundaryConditionsHandler, assembleFront); } if (assembled) { diff --git a/src/solvers/index.js b/src/solvers/index.js index 1e63eb8..3c92d3a 100644 --- a/src/solvers/index.js +++ b/src/solvers/index.js @@ -9,3 +9,8 @@ export { jacobiSolver } from "./jacobiSolver.js"; export { mathjsLUSolver } from "./mathjsLUSolver.js"; +export { solveLinearSystem, solveLinearSystemAsync } from "./linearSystemSolver.js"; +export { runFrontalSolver } from "./frontalSolver.js"; +export { newtonRaphson } from "./newtonRaphson.js"; +// WebGPUComputeEngine (webgpuJacobiSolver.js) is not re-exported, since importing it loads taichi.js; +// it is used only by the WebGPU worker diff --git a/src/solvers/jacobiSolver.js b/src/solvers/jacobiSolver.js index b2d196b..51bf902 100644 --- a/src/solvers/jacobiSolver.js +++ b/src/solvers/jacobiSolver.js @@ -8,8 +8,8 @@ */ // Internal imports -import { copyVector, euclideanNorm } from "../numerics/blasUtilities.js"; -import { SparseMatrix } from "../numerics/sparseMatrix.js"; +import { copyVector, euclideanNorm } from "../linalg/blasUtilities.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; /** * Function to solve a system of linear equations using the Jacobi iterative method (CPU synchronous version) diff --git a/src/solvers/linearSystemSolver.js b/src/solvers/linearSystemSolver.js index 7e869b1..af87c7d 100644 --- a/src/solvers/linearSystemSolver.js +++ b/src/solvers/linearSystemSolver.js @@ -8,8 +8,9 @@ */ // Internal imports -import { jacobiSolver, mathjsLUSolver } from "./index.js"; -import { SparseMatrix } from "../numerics/sparseMatrix.js"; +import { jacobiSolver } from "./jacobiSolver.js"; +import { mathjsLUSolver } from "./mathjsLUSolver.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; import { basicLog, debugLog, errorLog } from "../utilities/logging.js"; import * as Comlink from "../vendor/comlink.mjs"; diff --git a/src/solvers/mathjsLUSolver.js b/src/solvers/mathjsLUSolver.js index eee0759..dbc8cbb 100644 --- a/src/solvers/mathjsLUSolver.js +++ b/src/solvers/mathjsLUSolver.js @@ -12,7 +12,7 @@ import { create, all } from "mathjs"; const math = create(all); // Internal imports -import { SparseMatrix } from "../numerics/sparseMatrix.js"; +import { SparseMatrix } from "../linalg/sparseMatrix.js"; /** * Function to solve a system of linear equations using LU decomposition (mathjs) diff --git a/src/solvers/newtonRaphson.js b/src/solvers/newtonRaphson.js index b653ef7..faf701e 100644 --- a/src/solvers/newtonRaphson.js +++ b/src/solvers/newtonRaphson.js @@ -8,16 +8,18 @@ */ // Internal imports -import { euclideanNorm, axpy } from "../numerics/index.js"; +import { euclideanNorm, axpy } from "../linalg/index.js"; import { solveLinearSystem } from "./linearSystemSolver.js"; import { basicLog, errorLog } from "../utilities/logging.js"; import { runFrontalSolver } from "./frontalSolver.js"; -import { assembleFrontPropagationFront } from "../models/frontPropagation.js"; /** * Function to solve a system of non-linear equations using the Newton-Raphson method * @param {function} assembleMat - Matrix assembler based on the physical model - * @param {object} context - Context object containing simulation data and options + * @param {object} context - Context object containing simulation data and options. When + * `solverMethod` is "frontal", it must also contain `assembleFront` (element assembler, which + * receives the current iterate as `solutionVector`) and `boundaryConditionsHandler` + * (see runFrontalSolver) * @returns {object} An object containing: * - solutionVector: The solution vector (Float64Array) * - iterations: The number of iterations performed @@ -53,10 +55,9 @@ export function newtonRaphson(assembleMat, context = {}) { // Check if using frontal solver if (context.solverMethod === "frontal") { const frontalResult = runFrontalSolver( - assembleFrontPropagationFront, + (elementArgs) => context.assembleFront({ ...elementArgs, solutionVector }), context.meshData, - context.boundaryConditions, - { solutionVector, eikonalActivationFlag: context.eikonalActivationFlag }, + context.boundaryConditionsHandler, ); deltaX = frontalResult.solutionVector; } else { diff --git a/src/visualization/plotlyPlot.js b/src/visualization/plotlyPlot.js index 4cb382d..9e80e97 100644 --- a/src/visualization/plotlyPlot.js +++ b/src/visualization/plotlyPlot.js @@ -8,14 +8,9 @@ */ // Internal imports -import { - prepareMesh, - pointInsideQuadrilateral, - pointInsideTriangularElement, - computeNodeNeighbors, - getBoundarySegments, -} from "../mesh/meshUtils.js"; -import { BasisFunctions } from "../mesh/basisFunctions.js"; +import { prepareMesh, computeNodeNeighbors, getBoundarySegments } from "../mesh/meshUtils.js"; +import { pointInsideQuadrilateral, pointInsideTriangularElement } from "../elements/elementGeometry.js"; +import { BasisFunctions } from "../elements/basisFunctions.js"; /** * Function to create plots of the solution vector diff --git a/src/visualization/vtkPlot.js b/src/visualization/vtkPlot.js index 33396a4..8ef7a82 100644 --- a/src/visualization/vtkPlot.js +++ b/src/visualization/vtkPlot.js @@ -8,14 +8,9 @@ */ // Internal imports -import { - prepareMesh, - pointInsideQuadrilateral, - pointInsideTriangularElement, - computeNodeNeighbors, - getBoundarySegments, -} from "../mesh/meshUtils.js"; -import { BasisFunctions } from "../mesh/basisFunctions.js"; +import { prepareMesh, computeNodeNeighbors, getBoundarySegments } from "../mesh/meshUtils.js"; +import { pointInsideQuadrilateral, pointInsideTriangularElement } from "../elements/elementGeometry.js"; +import { BasisFunctions } from "../elements/basisFunctions.js"; import { errorLog } from "../utilities/logging.js"; let _vtkModules = null; diff --git a/tests/regression/EulerBernoulliBeam/regression.test.js b/tests/regression/EulerBernoulliBeam/regression.test.js index 2143b8b..9759057 100644 --- a/tests/regression/EulerBernoulliBeam/regression.test.js +++ b/tests/regression/EulerBernoulliBeam/regression.test.js @@ -24,7 +24,7 @@ */ import { FEAScriptModel } from "../../../src/FEAScript.js"; -import { assembleEulerBernoulliBeamMat } from "../../../src/models/eulerBernoulliBeam.js"; +import { assembleEulerBernoulliBeamMat } from "../../../src/physics/eulerBernoulliBeam.js"; import { basicLog, errorLog } from "../../../src/utilities/logging.js"; // Baseline values (see README.md for derivation and verification) diff --git a/tests/regression/HeatConduction2DFin/regression.test.js b/tests/regression/HeatConduction2DFin/regression.test.js index f02eabc..0f90dca 100644 --- a/tests/regression/HeatConduction2DFin/regression.test.js +++ b/tests/regression/HeatConduction2DFin/regression.test.js @@ -26,7 +26,7 @@ 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 { importGmshMesh } from "../../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url)); diff --git a/tests/regression/LidDrivenCavity2DCreepingFlow/regression.test.js b/tests/regression/LidDrivenCavity2DCreepingFlow/regression.test.js index f9ca7df..2983a40 100644 --- a/tests/regression/LidDrivenCavity2DCreepingFlow/regression.test.js +++ b/tests/regression/LidDrivenCavity2DCreepingFlow/regression.test.js @@ -22,7 +22,7 @@ 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 { importGmshMesh } from "../../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url)); diff --git a/tests/regression/SolidificationFront2D/regression.test.js b/tests/regression/SolidificationFront2D/regression.test.js index 8aa148b..5c5ce35 100644 --- a/tests/regression/SolidificationFront2D/regression.test.js +++ b/tests/regression/SolidificationFront2D/regression.test.js @@ -22,7 +22,7 @@ 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 { importGmshMesh } from "../../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url)); diff --git a/tests/unit/elements2D.test.js b/tests/unit/elements2D.test.js index 718bad1..80a2893 100644 --- a/tests/unit/elements2D.test.js +++ b/tests/unit/elements2D.test.js @@ -19,9 +19,9 @@ * Run: node tests/unit/elements2D.test.js (or npm test) */ -import { BasisFunctions } from "../../src/mesh/basisFunctions.js"; -import { getIntegrationPoints2D } from "../../src/mesh/meshUtils.js"; -import { NumericalIntegration } from "../../src/numerics/numericalIntegration.js"; +import { BasisFunctions } from "../../src/elements/basisFunctions.js"; +import { getIntegrationPoints2D } from "../../src/elements/isoparametricMapping.js"; +import { NumericalIntegration } from "../../src/elements/numericalIntegration.js"; import { basicLog, errorLog } from "../../src/utilities/logging.js"; basicLog(""); diff --git a/tests/unit/gmshReader.test.js b/tests/unit/gmshReader.test.js index 74476d4..c7d71a6 100644 --- a/tests/unit/gmshReader.test.js +++ b/tests/unit/gmshReader.test.js @@ -26,7 +26,7 @@ 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 { importGmshMesh } from "../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url)); diff --git a/tests/unit/plotInterpolation.test.js b/tests/unit/plotInterpolation.test.js index c139edf..628a9b6 100644 --- a/tests/unit/plotInterpolation.test.js +++ b/tests/unit/plotInterpolation.test.js @@ -25,7 +25,7 @@ import fs from "fs"; import path from "path"; import { fileURLToPath } from "url"; import { prepareMesh } from "../../src/mesh/meshUtils.js"; -import { importGmshMesh } from "../../src/readers/gmshReader.js"; +import { importGmshMesh } from "../../src/io/gmshReader.js"; import { plotInterpolatedSolution } from "../../src/visualization/plotlyPlot.js"; import { basicLog, errorLog } from "../../src/utilities/logging.js"; diff --git a/tests/unit/sparseMatrix.test.js b/tests/unit/sparseMatrix.test.js index 7bf0fe9..81b2148 100644 --- a/tests/unit/sparseMatrix.test.js +++ b/tests/unit/sparseMatrix.test.js @@ -19,7 +19,7 @@ * Run: node tests/unit/sparseMatrix.test.js (or npm test) */ -import { SparseMatrix } from "../../src/numerics/sparseMatrix.js"; +import { SparseMatrix } from "../../src/linalg/sparseMatrix.js"; import { jacobiSolver, mathjsLUSolver } from "../../src/solvers/index.js"; import { basicLog, errorLog } from "../../src/utilities/logging.js"; diff --git a/tests/verification/analytical/creepingFlow2DCouette/analytical.test.js b/tests/verification/analytical/creepingFlow2DCouette/analytical.test.js index 4ee9fe1..578db7c 100644 --- a/tests/verification/analytical/creepingFlow2DCouette/analytical.test.js +++ b/tests/verification/analytical/creepingFlow2DCouette/analytical.test.js @@ -21,7 +21,7 @@ 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 { importGmshMesh } from "../../../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url)); diff --git a/tests/verification/analytical/eulerBernoulliBeam1DCantilever/analytical.test.js b/tests/verification/analytical/eulerBernoulliBeam1DCantilever/analytical.test.js index 5e08afc..6986732 100644 --- a/tests/verification/analytical/eulerBernoulliBeam1DCantilever/analytical.test.js +++ b/tests/verification/analytical/eulerBernoulliBeam1DCantilever/analytical.test.js @@ -23,7 +23,7 @@ * Run: node tests/verification/analytical/eulerBernoulliBeam1DCantilever/analytical.test.js (or npm test) */ -import { assembleEulerBernoulliBeamMat } from "../../../../src/models/eulerBernoulliBeam.js"; +import { assembleEulerBernoulliBeamMat } from "../../../../src/physics/eulerBernoulliBeam.js"; import { prepareMesh } from "../../../../src/mesh/meshUtils.js"; import { solveLinearSystem } from "../../../../src/solvers/linearSystemSolver.js"; import { basicLog, errorLog } from "../../../../src/utilities/logging.js"; diff --git a/tests/verification/analytical/frontPropagation2DPlanarFront/analytical.test.js b/tests/verification/analytical/frontPropagation2DPlanarFront/analytical.test.js index c1d9595..7371e77 100644 --- a/tests/verification/analytical/frontPropagation2DPlanarFront/analytical.test.js +++ b/tests/verification/analytical/frontPropagation2DPlanarFront/analytical.test.js @@ -23,7 +23,7 @@ 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 { importGmshMesh } from "../../../../src/io/gmshReader.js"; import { basicLog, errorLog } from "../../../../src/utilities/logging.js"; const __dirname = path.dirname(fileURLToPath(import.meta.url));