Skip to content

Commit f4efb3f

Browse files
nikoschamnikoscham
committed
Update mesh generation and GMSH reader to support quadratic elements and improve boundary element handling; adjust thermal boundary conditions for better Jacobian calculations; disable contour labels in plot solution visualization
1 parent 015dec2 commit f4efb3f

File tree

4 files changed

+183
-98
lines changed

4 files changed

+183
-98
lines changed

‎src/mesh/meshGenerationScript.js‎

Lines changed: 82 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -66,8 +66,6 @@ export class meshGeneration {
6666
JSON.stringify(this.parsedMesh.nodalNumbering)
6767
);
6868

69-
console.log(this.parsedMesh.elementTypes[3]);
70-
7169
// Check if it has quadElements or triangleElements structure from gmshReader
7270
if (this.parsedMesh.elementTypes[3] || this.parsedMesh.elementTypes[10]) {
7371
// Map nodal numbering from GMSH format to FEAScript format for quad elements
@@ -175,7 +173,7 @@ export class meshGeneration {
175173
for (let elemIdx = 0; elemIdx < this.parsedMesh.nodalNumbering.length; elemIdx++) {
176174
const elemNodes = this.parsedMesh.nodalNumbering[elemIdx];
177175

178-
// For linear quadrilateral elements only (4 nodes)
176+
// For linear quadrilateral linear elements (4 nodes)
179177
if (elemNodes.length === 4) {
180178
// Check if both boundary nodes are in this element
181179
if (elemNodes.includes(node1) && elemNodes.includes(node2)) {
@@ -194,6 +192,11 @@ export class meshGeneration {
194192
` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`
195193
);
196194

195+
// Based on FEAScript linear quadrilateral numbering:
196+
// 1 --- 3
197+
// | |
198+
// 0 --- 2
199+
197200
if (
198201
(node1Index === 0 && node2Index === 2) ||
199202
(node1Index === 2 && node2Index === 0)
@@ -220,6 +223,82 @@ export class meshGeneration {
220223
debugLog(` These nodes form the RIGHT side (${side}) of element ${elemIdx}`);
221224
}
222225

226+
// Add the element and side to the boundary elements array
227+
this.parsedMesh.boundaryElements[prop.tag].push([elemIdx, side]);
228+
debugLog(
229+
` Added element-side pair [${elemIdx}, ${side}] to boundary tag ${prop.tag}`
230+
);
231+
foundElement = true;
232+
break;
233+
}
234+
} else if (elemNodes.length === 9) {
235+
// For quadratic quadrilateral elements (9 nodes)
236+
// Check if both boundary nodes are in this element
237+
if (elemNodes.includes(node1) && elemNodes.includes(node2)) {
238+
// Find which side of the element these nodes form
239+
let side;
240+
241+
const node1Index = elemNodes.indexOf(node1);
242+
const node2Index = elemNodes.indexOf(node2);
243+
244+
debugLog(
245+
` Found element ${elemIdx} containing boundary nodes. Element nodes: [${elemNodes.join(
246+
", "
247+
)}]`
248+
);
249+
debugLog(
250+
` Node ${node1} is at index ${node1Index}, Node ${node2} is at index ${node2Index} in the element`
251+
);
252+
253+
// Based on FEAScript quadratic quadrilateral numbering:
254+
// 2--5--8
255+
// | |
256+
// 1 4 7
257+
// | |
258+
// 0--3--6
259+
260+
if (
261+
(node1Index === 0 && node2Index === 6) ||
262+
(node1Index === 6 && node2Index === 0) ||
263+
(node1Index === 0 && node2Index === 3) ||
264+
(node1Index === 3 && node2Index === 0) ||
265+
(node1Index === 3 && node2Index === 6) ||
266+
(node1Index === 6 && node2Index === 3)
267+
) {
268+
side = 0; // Bottom side (nodes 0, 3, 6)
269+
debugLog(` These nodes form the BOTTOM side (${side}) of element ${elemIdx}`);
270+
} else if (
271+
(node1Index === 0 && node2Index === 2) ||
272+
(node1Index === 2 && node2Index === 0) ||
273+
(node1Index === 0 && node2Index === 1) ||
274+
(node1Index === 1 && node2Index === 0) ||
275+
(node1Index === 1 && node2Index === 2) ||
276+
(node1Index === 2 && node2Index === 1)
277+
) {
278+
side = 1; // Left side (nodes 0, 1, 2)
279+
debugLog(` These nodes form the LEFT side (${side}) of element ${elemIdx}`);
280+
} else if (
281+
(node1Index === 2 && node2Index === 8) ||
282+
(node1Index === 8 && node2Index === 2) ||
283+
(node1Index === 2 && node2Index === 5) ||
284+
(node1Index === 5 && node2Index === 2) ||
285+
(node1Index === 5 && node2Index === 8) ||
286+
(node1Index === 8 && node2Index === 5)
287+
) {
288+
side = 2; // Top side (nodes 2, 5, 8)
289+
debugLog(` These nodes form the TOP side (${side}) of element ${elemIdx}`);
290+
} else if (
291+
(node1Index === 6 && node2Index === 8) ||
292+
(node1Index === 8 && node2Index === 6) ||
293+
(node1Index === 6 && node2Index === 7) ||
294+
(node1Index === 7 && node2Index === 6) ||
295+
(node1Index === 7 && node2Index === 8) ||
296+
(node1Index === 8 && node2Index === 7)
297+
) {
298+
side = 3; // Right side (nodes 6, 7, 8)
299+
debugLog(` These nodes form the RIGHT side (${side}) of element ${elemIdx}`);
300+
}
301+
223302
// Add the element and side to the boundary elements array
224303
this.parsedMesh.boundaryElements[prop.tag].push([elemIdx, side]);
225304
debugLog(

‎src/readers/gmshReaderScript.js‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -221,7 +221,7 @@ const importGmshQuadTri = async (file) => {
221221
const elementTag = parseInt(parts[0], 10);
222222
const nodeIndices = parts.slice(1).map((idx) => parseInt(idx, 10));
223223

224-
if (currentElementBlock.elementType === 1) {
224+
if (currentElementBlock.elementType === 1 || currentElementBlock.elementType === 8) {
225225
const physicalTag = currentElementBlock.tag;
226226

227227
if (!boundaryElementsByTag[physicalTag]) {
@@ -235,11 +235,14 @@ const importGmshQuadTri = async (file) => {
235235
result.boundaryNodePairs[physicalTag] = [];
236236
}
237237
result.boundaryNodePairs[physicalTag].push(nodeIndices);
238-
} else if (currentElementBlock.elementType === 2) { // Linear triangle elements (3 nodes)
238+
} else if (currentElementBlock.elementType === 2) {
239+
// Linear triangle elements (3 nodes)
239240
result.nodalNumbering.triangleElements.push(nodeIndices);
240-
} else if (currentElementBlock.elementType === 3) { // Linear quadrilateral elements (4 nodes)
241+
} else if (currentElementBlock.elementType === 3) {
242+
// Linear quadrilateral elements (4 nodes)
241243
result.nodalNumbering.quadElements.push(nodeIndices);
242-
} else if (currentElementBlock.elementType === 10) { // Quadratic quadrilateral elements (9 nodes)
244+
} else if (currentElementBlock.elementType === 10) {
245+
// Quadratic quadrilateral elements (9 nodes)
243246
result.nodalNumbering.quadElements.push(nodeIndices);
244247
}
245248

‎src/solvers/thermalBoundaryConditionsScript.js‎

Lines changed: 92 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -213,6 +213,11 @@ export class ThermalBoundaryConditions {
213213
}
214214

215215
const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1;
216+
debugLog(
217+
` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${
218+
elementIndex + 1
219+
}, local node ${nodeIndex + 1})`
220+
);
216221
residualVector[globalNodeIndex] += -convectionCoeff * extTemp;
217222
jacobianMatrix[globalNodeIndex][globalNodeIndex] += convectionCoeff;
218223
});
@@ -267,61 +272,61 @@ export class ThermalBoundaryConditions {
267272
let basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi;
268273
let basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta;
269274

270-
let xCoordinates = 0;
275+
// Calculate boundary Jacobian
271276
let ksiDerivX = 0;
277+
let ksiDerivY = 0;
278+
let etaDerivX = 0;
272279
let etaDerivY = 0;
273280
const numNodes = this.nop[elementIndex].length;
274-
275281
for (let nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
276282
const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1;
277-
xCoordinates += nodesXCoordinates[globalNodeIndex] * basisFunction[nodeIndex];
278-
ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex];
279-
etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex];
283+
284+
// For boundaries along Ksi (horizontal), use Ksi derivatives
285+
if (side === 0 || side === 2) {
286+
ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex];
287+
ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex];
288+
}
289+
// For boundaries along Eta (vertical), use Eta derivatives
290+
else if (side === 1 || side === 3) {
291+
etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex];
292+
etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex];
293+
}
280294
}
281295

296+
// Compute boundary Jacobian (length of tangent vector)
297+
const detJacobian =
298+
side === 0 || side === 2
299+
? Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2)
300+
: Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2);
301+
282302
for (
283303
let localNodeIndex = firstNodeIndex;
284304
localNodeIndex < lastNodeIndex;
285305
localNodeIndex += nodeIncrement
286306
) {
287307
let globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1;
308+
debugLog(
309+
` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${
310+
elementIndex + 1
311+
}, local node ${localNodeIndex + 1})`
312+
);
288313

289-
if (side === 0 || side === 2) {
290-
// Horizontal boundaries of the domain (assuming a rectangular domain)
291-
residualVector[globalNodeIndex] +=
292-
-gaussWeights[0] * ksiDerivX * basisFunction[localNodeIndex] * convectionCoeff * extTemp;
293-
294-
for (
295-
let localNodeIndex2 = firstNodeIndex;
296-
localNodeIndex2 < lastNodeIndex;
297-
localNodeIndex2 += nodeIncrement
298-
) {
299-
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
300-
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
301-
-gaussWeights[0] *
302-
ksiDerivX *
303-
basisFunction[localNodeIndex] *
304-
basisFunction[localNodeIndex2] *
305-
convectionCoeff;
306-
}
307-
} else if (side === 1 || side === 3) {
308-
// Vertical boundaries of the domain (assuming a rectangular domain)
309-
residualVector[globalNodeIndex] +=
310-
-gaussWeights[0] * etaDerivY * basisFunction[localNodeIndex] * convectionCoeff * extTemp;
314+
// Apply boundary condition with proper Jacobian for all sides
315+
residualVector[globalNodeIndex] +=
316+
-gaussWeights[0] * detJacobian * basisFunction[localNodeIndex] * convectionCoeff * extTemp;
311317

312-
for (
313-
let localNodeIndex2 = firstNodeIndex;
314-
localNodeIndex2 < lastNodeIndex;
315-
localNodeIndex2 += nodeIncrement
316-
) {
317-
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
318-
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
319-
-gaussWeights[0] *
320-
etaDerivY *
321-
basisFunction[localNodeIndex] *
322-
basisFunction[localNodeIndex2] *
323-
convectionCoeff;
324-
}
318+
for (
319+
let localNodeIndex2 = firstNodeIndex;
320+
localNodeIndex2 < lastNodeIndex;
321+
localNodeIndex2 += nodeIncrement
322+
) {
323+
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
324+
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
325+
-gaussWeights[0] *
326+
detJacobian *
327+
basisFunction[localNodeIndex] *
328+
basisFunction[localNodeIndex2] *
329+
convectionCoeff;
325330
}
326331
}
327332
} else if (this.elementOrder === "quadratic") {
@@ -363,68 +368,66 @@ export class ThermalBoundaryConditions {
363368
let basisFunction = basisFunctionsAndDerivatives.basisFunction;
364369
let basisFunctionDerivKsi = basisFunctionsAndDerivatives.basisFunctionDerivKsi;
365370
let basisFunctionDerivEta = basisFunctionsAndDerivatives.basisFunctionDerivEta;
366-
let xCoordinates = 0;
371+
372+
// Calculate boundary Jacobian
367373
let ksiDerivX = 0;
374+
let ksiDerivY = 0;
375+
let etaDerivX = 0;
368376
let etaDerivY = 0;
369377
const numNodes = this.nop[elementIndex].length;
370378
for (let nodeIndex = 0; nodeIndex < numNodes; nodeIndex++) {
371-
xCoordinates +=
372-
nodesXCoordinates[this.nop[elementIndex][nodeIndex] - 1] * basisFunction[nodeIndex];
373-
ksiDerivX +=
374-
nodesXCoordinates[this.nop[elementIndex][nodeIndex] - 1] *
375-
basisFunctionDerivKsi[nodeIndex];
376-
etaDerivY +=
377-
nodesYCoordinates[this.nop[elementIndex][nodeIndex] - 1] *
378-
basisFunctionDerivEta[nodeIndex];
379+
const globalNodeIndex = this.nop[elementIndex][nodeIndex] - 1;
380+
381+
// For boundaries along Ksi (horizontal), use Ksi derivatives
382+
if (side === 0 || side === 2) {
383+
ksiDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex];
384+
ksiDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivKsi[nodeIndex];
385+
}
386+
// For boundaries along Eta (vertical), use Eta derivatives
387+
else if (side === 1 || side === 3) {
388+
etaDerivX += nodesXCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex];
389+
etaDerivY += nodesYCoordinates[globalNodeIndex] * basisFunctionDerivEta[nodeIndex];
390+
}
379391
}
392+
393+
// Compute boundary Jacobian (length of tangent vector)
394+
const detJacobian =
395+
side === 0 || side === 2
396+
? Math.sqrt(ksiDerivX ** 2 + ksiDerivY ** 2)
397+
: Math.sqrt(etaDerivX ** 2 + etaDerivY ** 2);
398+
380399
for (
381400
let localNodeIndex = firstNodeIndex;
382401
localNodeIndex < lastNodeIndex;
383402
localNodeIndex += nodeIncrement
384403
) {
385404
let globalNodeIndex = this.nop[elementIndex][localNodeIndex] - 1;
386-
if (side === 0 || side === 2) {
387-
// Horizontal boundaries of the domain (assuming a rectangular domain)
388-
residualVector[globalNodeIndex] +=
389-
-gaussWeights[gaussPointIndex] *
390-
ksiDerivX *
391-
basisFunction[localNodeIndex] *
392-
convectionCoeff *
393-
extTemp;
394-
for (
395-
let localNodeIndex2 = firstNodeIndex;
396-
localNodeIndex2 < lastNodeIndex;
397-
localNodeIndex2 += nodeIncrement
398-
) {
399-
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
400-
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
401-
-gaussWeights[gaussPointIndex] *
402-
ksiDerivX *
403-
basisFunction[localNodeIndex] *
404-
basisFunction[localNodeIndex2] *
405-
convectionCoeff;
406-
}
407-
} else if (side === 1 || side === 3) {
408-
// Vertical boundaries of the domain (assuming a rectangular domain)
409-
residualVector[globalNodeIndex] +=
405+
debugLog(
406+
` - Applied convection boundary condition to node ${globalNodeIndex + 1} (element ${
407+
elementIndex + 1
408+
}, local node ${localNodeIndex + 1})`
409+
);
410+
411+
// Apply boundary condition with proper Jacobian for all sides
412+
residualVector[globalNodeIndex] +=
413+
-gaussWeights[gaussPointIndex] *
414+
detJacobian *
415+
basisFunction[localNodeIndex] *
416+
convectionCoeff *
417+
extTemp;
418+
419+
for (
420+
let localNodeIndex2 = firstNodeIndex;
421+
localNodeIndex2 < lastNodeIndex;
422+
localNodeIndex2 += nodeIncrement
423+
) {
424+
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
425+
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
410426
-gaussWeights[gaussPointIndex] *
411-
etaDerivY *
427+
detJacobian *
412428
basisFunction[localNodeIndex] *
413-
convectionCoeff *
414-
extTemp;
415-
for (
416-
let localNodeIndex2 = firstNodeIndex;
417-
localNodeIndex2 < lastNodeIndex;
418-
localNodeIndex2 += nodeIncrement
419-
) {
420-
let globalNodeIndex2 = this.nop[elementIndex][localNodeIndex2] - 1;
421-
jacobianMatrix[globalNodeIndex][globalNodeIndex2] +=
422-
-gaussWeights[gaussPointIndex] *
423-
etaDerivY *
424-
basisFunction[localNodeIndex] *
425-
basisFunction[localNodeIndex2] *
426-
convectionCoeff;
427-
}
429+
basisFunction[localNodeIndex2] *
430+
convectionCoeff;
428431
}
429432
}
430433
}

0 commit comments

Comments
 (0)