Skip to content

Commit 0812dd7

Browse files
nikoschamnikoscham
committed
Update nodal numbering handling in mesh generation and enhance boundary element processing in GMSH reader - Still in debugging phase
1 parent ac09141 commit 0812dd7

File tree

3 files changed

+176
-114
lines changed

3 files changed

+176
-114
lines changed

‎src/mesh/meshGenerationScript.js‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,155 @@ export class meshGeneration {
5151
generateMesh() {
5252
// If pre-parsed mesh data is provided, use it directly
5353
if (this.parsedMesh) {
54+
// Process the nodalNumbering from gmshReader format to the format expected by the solver
55+
if (this.parsedMesh.nodalNumbering) {
56+
if (
57+
typeof this.parsedMesh.nodalNumbering === "object" &&
58+
!Array.isArray(this.parsedMesh.nodalNumbering)
59+
) {
60+
// Store the nodal numbering structure before converting
61+
const quadElements = this.parsedMesh.nodalNumbering.quadElements || [];
62+
const triangleElements = this.parsedMesh.nodalNumbering.triangleElements || [];
63+
64+
console.log("Parsed mesh nodal numbering:", this.parsedMesh.nodalNumbering);
65+
66+
// Check if it has quadElements or triangleElements structure from gmshReader
67+
if (quadElements && quadElements.length > 0) {
68+
// Map nodal numbering from GMSH format to FEAScript format for quad elements
69+
const mappedNodalNumbering = [];
70+
71+
for (let elemIdx = 0; elemIdx < quadElements.length; elemIdx++) {
72+
const gmshNodes = quadElements[elemIdx];
73+
const feaScriptNodes = new Array(gmshNodes.length);
74+
75+
// Simple mapping for linear quad elements (4 nodes)
76+
// GMSH: FEAScript:
77+
// 3 --- 2 1 --- 3
78+
// | | | |
79+
// 0 --- 1 0 --- 2
80+
81+
feaScriptNodes[0] = gmshNodes[0]; // 0 -> 0
82+
feaScriptNodes[1] = gmshNodes[3]; // 3 -> 1
83+
feaScriptNodes[2] = gmshNodes[1]; // 1 -> 2
84+
feaScriptNodes[3] = gmshNodes[2]; // 2 -> 3
85+
86+
mappedNodalNumbering.push(feaScriptNodes);
87+
}
88+
89+
this.parsedMesh.nodalNumbering = mappedNodalNumbering;
90+
} else if (triangleElements && triangleElements.length > 0) {
91+
this.parsedMesh.nodalNumbering = triangleElements;
92+
}
93+
94+
console.log("Mapped nodal numbering:", this.parsedMesh.nodalNumbering);
95+
96+
// Process boundary elements if they exist and if physical property mapping exists
97+
if (this.parsedMesh.physicalPropMap && this.parsedMesh.boundaryElements) {
98+
// Check if boundary elements need to be processed
99+
if (
100+
Array.isArray(this.parsedMesh.boundaryElements) &&
101+
this.parsedMesh.boundaryElements.length > 0 &&
102+
this.parsedMesh.boundaryElements[0] === undefined
103+
) {
104+
// Create a new array without the empty first element
105+
const fixedBoundaryElements = [];
106+
for (let i = 1; i < this.parsedMesh.boundaryElements.length; i++) {
107+
if (this.parsedMesh.boundaryElements[i]) {
108+
fixedBoundaryElements.push(this.parsedMesh.boundaryElements[i]);
109+
}
110+
}
111+
this.parsedMesh.boundaryElements = fixedBoundaryElements;
112+
}
113+
114+
// If boundary node pairs exist but boundary elements haven't been processed
115+
if (this.parsedMesh.boundaryNodePairs && !this.parsedMesh.boundaryElementsProcessed) {
116+
this.parsedMesh.boundaryElements = [];
117+
118+
this.parsedMesh.physicalPropMap.forEach((prop) => {
119+
if (prop.dimension === 1) {
120+
const boundaryNodes = this.parsedMesh.boundaryNodePairs[prop.tag] || [];
121+
122+
if (boundaryNodes.length > 0) {
123+
// Initialize boundary element array for this tag
124+
if (!this.parsedMesh.boundaryElements[prop.tag]) {
125+
this.parsedMesh.boundaryElements[prop.tag] = [];
126+
}
127+
128+
// For each boundary line segment (2 nodes)
129+
boundaryNodes.forEach((nodesPair) => {
130+
// Find which element contains both nodes of this boundary
131+
let foundElement = false;
132+
for (let elemIdx = 0; elemIdx < this.parsedMesh.nodalNumbering.length; elemIdx++) {
133+
const elemNodes = this.parsedMesh.nodalNumbering[elemIdx];
134+
135+
// Check if both boundary nodes are in this element
136+
if (elemNodes.includes(nodesPair[0]) && elemNodes.includes(nodesPair[1])) {
137+
// Find which side of the element these nodes form
138+
let side;
139+
140+
// Using FEAScript orientation for sides
141+
// Side 0 (Bottom): nodes 0-2
142+
// Side 1 (Left): nodes 0-1
143+
// Side 2 (Top): nodes 1-3
144+
// Side 3 (Right): nodes 2-3
145+
146+
if (
147+
(elemNodes.indexOf(nodesPair[0]) === 0 &&
148+
elemNodes.indexOf(nodesPair[1]) === 2) ||
149+
(elemNodes.indexOf(nodesPair[1]) === 0 && elemNodes.indexOf(nodesPair[0]) === 2)
150+
) {
151+
side = 0; // Bottom side
152+
} else if (
153+
(elemNodes.indexOf(nodesPair[0]) === 0 &&
154+
elemNodes.indexOf(nodesPair[1]) === 1) ||
155+
(elemNodes.indexOf(nodesPair[1]) === 0 && elemNodes.indexOf(nodesPair[0]) === 1)
156+
) {
157+
side = 1; // Left side
158+
} else if (
159+
(elemNodes.indexOf(nodesPair[0]) === 1 &&
160+
elemNodes.indexOf(nodesPair[1]) === 3) ||
161+
(elemNodes.indexOf(nodesPair[1]) === 1 && elemNodes.indexOf(nodesPair[0]) === 3)
162+
) {
163+
side = 2; // Top side
164+
} else if (
165+
(elemNodes.indexOf(nodesPair[0]) === 2 &&
166+
elemNodes.indexOf(nodesPair[1]) === 3) ||
167+
(elemNodes.indexOf(nodesPair[1]) === 2 && elemNodes.indexOf(nodesPair[0]) === 3)
168+
) {
169+
side = 3; // Right side
170+
}
171+
172+
// Add to boundary elements
173+
this.parsedMesh.boundaryElements[prop.tag].push([elemIdx, side]);
174+
foundElement = true;
175+
break;
176+
}
177+
}
178+
});
179+
}
180+
}
181+
});
182+
183+
// Mark as processed
184+
this.parsedMesh.boundaryElementsProcessed = true;
185+
186+
// Fix boundary elements array - remove undefined entries
187+
if (
188+
this.parsedMesh.boundaryElements.length > 0 &&
189+
this.parsedMesh.boundaryElements[0] === undefined
190+
) {
191+
const fixedBoundaryElements = [];
192+
for (let i = 1; i < this.parsedMesh.boundaryElements.length; i++) {
193+
if (this.parsedMesh.boundaryElements[i]) {
194+
fixedBoundaryElements.push(this.parsedMesh.boundaryElements[i]);
195+
}
196+
}
197+
this.parsedMesh.boundaryElements = fixedBoundaryElements;
198+
}
199+
}
200+
}
201+
}
202+
}
54203
return this.parsedMesh;
55204
} else {
56205
// Validate required geometry parameters based on mesh dimension
@@ -179,6 +328,9 @@ export class meshGeneration {
179328
// Find boundary elements
180329
const boundaryElements = this.findBoundaryElements();
181330

331+
console.log("x coordinates of nodes:", nodesXCoordinates);
332+
console.log("y coordinates of nodes:", nodesYCoordinates);
333+
182334
// Return x and y coordinates of nodes, total nodes, NOP array, and boundary elements
183335
return {
184336
nodesXCoordinates,

‎src/readers/gmshReaderScript.js‎

Lines changed: 8 additions & 98 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ const importGmshQuadTri = async (file) => {
2626
},
2727
boundaryElements: [],
2828
boundaryConditions: [],
29+
boundaryNodePairs: {}, // Store boundary node pairs for processing in meshGenerationScript
2930
gmshV: 0,
3031
ascii: false,
3132
fltBytes: "8",
@@ -228,6 +229,12 @@ const importGmshQuadTri = async (file) => {
228229
}
229230

230231
boundaryElementsByTag[physicalTag].push(nodeIndices);
232+
233+
// Store boundary node pairs for later processing in meshGenerationScript
234+
if (!result.boundaryNodePairs[physicalTag]) {
235+
result.boundaryNodePairs[physicalTag] = [];
236+
}
237+
result.boundaryNodePairs[physicalTag].push(nodeIndices);
231238
} else if (currentElementBlock.elementType === 2) {
232239
result.nodalNumbering.triangleElements.push(nodeIndices);
233240
} else if (currentElementBlock.elementType === 3) {
@@ -246,6 +253,7 @@ const importGmshQuadTri = async (file) => {
246253
lineIndex++;
247254
}
248255

256+
// Store boundary conditions information
249257
result.physicalPropMap.forEach((prop) => {
250258
if (prop.dimension === 1) {
251259
const boundaryNodes = boundaryElementsByTag[prop.tag] || [];
@@ -256,108 +264,10 @@ const importGmshQuadTri = async (file) => {
256264
tag: prop.tag,
257265
nodes: boundaryNodes,
258266
});
259-
260-
// Initialize boundary element array based on the number of PhysicalNames
261-
if (!result.boundaryElements[prop.tag]) {
262-
result.boundaryElements[prop.tag] = [];
263-
}
264-
265-
// For each boundary line segment (2 nodes)
266-
boundaryNodes.forEach((nodesPair) => {
267-
// Find which quad element contains both nodes of this boundary
268-
let foundElement = false;
269-
for (let elemIdx = 0; elemIdx < result.nodalNumbering.quadElements.length; elemIdx++) {
270-
const elemNodes = result.nodalNumbering.quadElements[elemIdx];
271-
272-
// Check if both boundary nodes are in this element
273-
if (elemNodes.includes(nodesPair[0]) && elemNodes.includes(nodesPair[1])) {
274-
// Find which side of the element these nodes form
275-
let side;
276-
277-
// For linear quadrilateral elements, the Gmsh local numbering is:
278-
// 3 --- 2
279-
// | |
280-
// 0 --- 1
281-
282-
// Check which side these nodes belong to
283-
if (
284-
(elemNodes.indexOf(nodesPair[0]) === 0 && elemNodes.indexOf(nodesPair[1]) === 1) ||
285-
(elemNodes.indexOf(nodesPair[1]) === 0 && elemNodes.indexOf(nodesPair[0]) === 1)
286-
) {
287-
side = 0; // Bottom side
288-
} else if (
289-
(elemNodes.indexOf(nodesPair[0]) === 1 && elemNodes.indexOf(nodesPair[1]) === 2) ||
290-
(elemNodes.indexOf(nodesPair[1]) === 1 && elemNodes.indexOf(nodesPair[0]) === 2)
291-
) {
292-
side = 1; // Right side
293-
} else if (
294-
(elemNodes.indexOf(nodesPair[0]) === 2 && elemNodes.indexOf(nodesPair[1]) === 3) ||
295-
(elemNodes.indexOf(nodesPair[1]) === 2 && elemNodes.indexOf(nodesPair[0]) === 3)
296-
) {
297-
side = 2; // Top side
298-
} else if (
299-
(elemNodes.indexOf(nodesPair[0]) === 3 && elemNodes.indexOf(nodesPair[1]) === 0) ||
300-
(elemNodes.indexOf(nodesPair[1]) === 3 && elemNodes.indexOf(nodesPair[0]) === 0)
301-
) {
302-
side = 3; // Left side
303-
}
304-
305-
// Add to boundary elements
306-
result.boundaryElements[prop.tag].push([elemIdx, side]);
307-
foundElement = true;
308-
break;
309-
}
310-
}
311-
});
312267
}
313268
}
314269
});
315270

316-
// Remap nodal numbering from Gmsh format to FEAScript format for quadrilateral elements
317-
if (result.nodalNumbering.quadElements.length > 0) {
318-
/*
319-
* Gmsh quadrilateral node numbering (linear elements):
320-
*
321-
* 3__ __2
322-
* | |
323-
* |__ __|
324-
* 0 1
325-
*
326-
* FEAScript quadrilateral node numbering (linear elements):
327-
*
328-
* 1__ __3
329-
* | |
330-
* |__ __|
331-
* 0 2
332-
*
333-
*/
334-
// Mapping: Gmsh → FEAScript
335-
// 0 → 0 (bottom left)
336-
// 1 → 2 (bottom right)
337-
// 2 → 3 (top right)
338-
// 3 → 1 (top left)
339-
const gmshToFEAMap = [0, 2, 3, 1];
340-
341-
for (let i = 0; i < result.nodalNumbering.quadElements.length; i++) {
342-
const originalNodes = [...result.nodalNumbering.quadElements[i]];
343-
for (let j = 0; j < 4; j++) {
344-
result.nodalNumbering.quadElements[i][gmshToFEAMap[j]] = originalNodes[j];
345-
}
346-
}
347-
}
348-
349-
// Fix boundary elements array - remove the empty first element
350-
if (result.boundaryElements.length > 0 && result.boundaryElements[0] === undefined) {
351-
// Create a new array without the empty first element
352-
const fixedBoundaryElements = [];
353-
for (let i = 1; i < result.boundaryElements.length; i++) {
354-
if (result.boundaryElements[i]) {
355-
fixedBoundaryElements.push(result.boundaryElements[i]);
356-
}
357-
}
358-
result.boundaryElements = fixedBoundaryElements;
359-
}
360-
361271
return result;
362272
};
363273

‎src/solvers/solidHeatTransferScript.js‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,7 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
8282
}
8383

8484
// Initialize variables for matrix assembly
85-
let localNodalNumbers = []; // Local nodal numbering
85+
let localToGlobalMap = []; // Maps local element node indices to global mesh node indices
8686
let gaussPoints = []; // Gauss points
8787
let gaussWeights = []; // Gauss weights
8888
let basisFunction = []; // Basis functions
@@ -133,7 +133,7 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
133133
for (let elementIndex = 0; elementIndex < totalElements; elementIndex++) {
134134
for (let localNodeIndex = 0; localNodeIndex < numNodes; localNodeIndex++) {
135135
// Subtract 1 from nop in order to start numbering from 0
136-
localNodalNumbers[localNodeIndex] = nop[elementIndex][localNodeIndex] - 1;
136+
localToGlobalMap[localNodeIndex] = nop[elementIndex][localNodeIndex] - 1;
137137
}
138138

139139
// Loop over Gauss points
@@ -152,9 +152,9 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
152152
// Isoparametric mapping
153153
for (let localNodeIndex = 0; localNodeIndex < numNodes; localNodeIndex++) {
154154
xCoordinates +=
155-
nodesXCoordinates[localNodalNumbers[localNodeIndex]] * basisFunction[localNodeIndex];
155+
nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex];
156156
ksiDerivX +=
157-
nodesXCoordinates[localNodalNumbers[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
157+
nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
158158
detJacobian = ksiDerivX;
159159
}
160160

@@ -165,12 +165,12 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
165165

166166
// Computation of Galerkin's residuals and Jacobian matrix
167167
for (let localNodeIndex1 = 0; localNodeIndex1 < numNodes; localNodeIndex1++) {
168-
let globalNodeIndex1 = localNodalNumbers[localNodeIndex1];
168+
let localToGlobalMap1 = localToGlobalMap[localNodeIndex1];
169169
// residualVector is zero for this case
170170

171171
for (let localNodeIndex2 = 0; localNodeIndex2 < numNodes; localNodeIndex2++) {
172-
let globalNodeIndex2 = localNodalNumbers[localNodeIndex2];
173-
jacobianMatrix[globalNodeIndex1][globalNodeIndex2] +=
172+
let localToGlobalMap2 = localToGlobalMap[localNodeIndex2];
173+
jacobianMatrix[localToGlobalMap1][localToGlobalMap2] +=
174174
-gaussWeights[gaussPointIndex1] *
175175
detJacobian *
176176
(basisFunctionDerivX[localNodeIndex1] * basisFunctionDerivX[localNodeIndex2]);
@@ -198,17 +198,17 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
198198
// Isoparametric mapping
199199
for (let localNodeIndex = 0; localNodeIndex < numNodes; localNodeIndex++) {
200200
xCoordinates +=
201-
nodesXCoordinates[localNodalNumbers[localNodeIndex]] * basisFunction[localNodeIndex];
201+
nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex];
202202
yCoordinates +=
203-
nodesYCoordinates[localNodalNumbers[localNodeIndex]] * basisFunction[localNodeIndex];
203+
nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunction[localNodeIndex];
204204
ksiDerivX +=
205-
nodesXCoordinates[localNodalNumbers[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
205+
nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
206206
etaDerivX +=
207-
nodesXCoordinates[localNodalNumbers[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex];
207+
nodesXCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex];
208208
ksiDerivY +=
209-
nodesYCoordinates[localNodalNumbers[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
209+
nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivKsi[localNodeIndex];
210210
etaDerivY +=
211-
nodesYCoordinates[localNodalNumbers[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex];
211+
nodesYCoordinates[localToGlobalMap[localNodeIndex]] * basisFunctionDerivEta[localNodeIndex];
212212
detJacobian = meshDimension === "2D" ? ksiDerivX * etaDerivY - etaDerivX * ksiDerivY : ksiDerivX;
213213
}
214214

@@ -226,12 +226,12 @@ export function assembleSolidHeatTransferMat(meshConfig, boundaryConditions) {
226226

227227
// Computation of Galerkin's residuals and Jacobian matrix
228228
for (let localNodeIndex1 = 0; localNodeIndex1 < numNodes; localNodeIndex1++) {
229-
let globalNodeIndex1 = localNodalNumbers[localNodeIndex1];
229+
let localToGlobalMap1 = localToGlobalMap[localNodeIndex1];
230230
// residualVector is zero for this case
231231

232232
for (let localNodeIndex2 = 0; localNodeIndex2 < numNodes; localNodeIndex2++) {
233-
let globalNodeIndex2 = localNodalNumbers[localNodeIndex2];
234-
jacobianMatrix[globalNodeIndex1][globalNodeIndex2] +=
233+
let localToGlobalMap2 = localToGlobalMap[localNodeIndex2];
234+
jacobianMatrix[localToGlobalMap1][localToGlobalMap2] +=
235235
-gaussWeights[gaussPointIndex1] *
236236
gaussWeights[gaussPointIndex2] *
237237
detJacobian *

0 commit comments

Comments
 (0)