From 52102f39fb5032360d30fc49961f4e61a3182a11 Mon Sep 17 00:00:00 2001 From: iagtirap Date: Wed, 25 Mar 2026 13:19:25 +0100 Subject: [PATCH 1/7] Add third order extrusion in prep for recursive approach --- pyhope/mesh/extrude/mesh_extrude.py | 84 ++++++++++++++++++++++++++++- 1 file changed, 82 insertions(+), 2 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 7a2cfa08..3dc3e938 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -82,7 +82,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: elif not [cell_block for cell_block in mesh.cells if cell_block.type in gmshCellTypes.cellTypes2D]: hopout.error('Mesh contains no suitable surface cells for extrusion, exiting...') - if nGeo > 2: + if nGeo > 3: hopout.error(f'nGeo = {nGeo} not supported for mesh extrusion') hopout.info('Extruding surface to volume mesh') @@ -433,6 +433,46 @@ def extrude_pris(nodes: np.ndarray, newPoints[offsetCurr+10:offsetCurr+11, :] = points[ 4] + 0.5*(shiftCurr+shiftPrev) newPoints[offsetCurr+11:offsetCurr+12, :] = points[ 5] + 0.5*(shiftCurr+shiftPrev) + case 3: + nFaceDOFs = round((order+1)*(order+2)/2.) + nNewDOFs = order*nFaceDOFs + newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) + + # Hardcoded meshio(triangle10) -> meshio(wedge40) layer mapping for NGeo=3. + # Each entry is the wedge40 node position for face node q at layer k. + layer0 = np.array((0, 1, 2, 6, 7, 8, 9, 10, 11, 37), dtype=np.int64) + layer1 = np.array((18, 20, 22, 24, 25, 28, 29, 32, 33, 38), dtype=np.int64) + layer2 = np.array((19, 21, 23, 27, 26, 31, 30, 35, 34, 39), dtype=np.int64) + layer3 = np.array((3, 4, 5, 12, 13, 14, 15, 16, 17, 36), dtype=np.int64) + layers = (layer0, layer1, layer2, layer3) + + # Append the bottom layer of the first element, then stack all the other elements + for i in range(shifts.shape[0]-1): + # Calculate offset for current layer indices + offsetCurr = i*nNewDOFs + shiftCurr = shifts[i+1, :] + shiftPrev = shifts[i , :] + + # Bottom layer: from original face for first layer, else from previous top layer. + for q in range(nFaceDOFs): + idxBot = int(layer0[q]) + idxTop = int(layer3[q]) + newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + + # New points for layers k=1..3 + p = 0 + for k in (1, 2, 3): + alpha = k/order + shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr + layerK = layers[k] + for q in range(nFaceDOFs): + idx = nPoints + offsetCurr + p + pos = int(layerK[q]) + + newNodes[i][pos] = idx + newPoints[offsetCurr + p, :] = points[q] + shiftK + p += 1 + # FIXME: Implement the other orders case _: raise ValueError(f'Extrusion not implemented for NGeo={order}') @@ -470,7 +510,7 @@ def extrude_hexa(nodes: np.ndarray, for i in range(shifts.shape[0]-1): # Calculate offset for current layer indices - offsetCurr = i *18 + offsetCurr = i shiftCurr = shifts[i+1, :] shiftPrev = shifts[i , :] @@ -497,6 +537,46 @@ def extrude_hexa(nodes: np.ndarray, newNodes[i][26:27] = nPoints + 17 + offsetCurr newPoints[offsetCurr+17:offsetCurr+18, :] = points[ 8] + 0.5*(shiftCurr+shiftPrev) + case 3: + nFaceDOFs = (order+1)**2 + nNewDOFs = order*nFaceDOFs + newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) + + # Hardcoded meshio(quad16) -> meshio(hexa64) layer mapping for NGeo=3. + # Each entry is the hexa64 node position for face node q at layer k. + layer0 = np.array(( 0, 1, 2, 3, 8, 9, 10, 11, 12, 13, 15, 14, 48, 51, 50, 49), dtype=np.int64) + layer1 = np.array((24, 26, 28, 30, 40, 41, 36, 37, 44, 45, 35, 32, 56, 57, 58, 59), dtype=np.int64) + layer2 = np.array((25, 27, 29, 31, 43, 42, 39, 38, 47, 46, 34, 33, 60, 61, 62, 63), dtype=np.int64) + layer3 = np.array(( 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 23, 22, 52, 53, 54, 55), dtype=np.int64) + layers = (layer0, layer1, layer2, layer3) + + # Append the bottom layer of the first element, then stack all the other elements + for i in range(shifts.shape[0]-1): + # Calculate offset for current layer indices + offsetCurr = i*nNewDOFs + shiftCurr = shifts[i+1, :] + shiftPrev = shifts[i , :] + + # Bottom layer: from original face for first layer, else from previous top layer. + for q in range(nFaceDOFs): + idxBot = int(layer0[q]) + idxTop = int(layer3[q]) + newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + + # New points for layers k=1..3 + p = 0 + for k in (1, 2, 3): + alpha = k/order + shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr + layerK = layers[k] + for q in range(nFaceDOFs): + idx = nPoints + offsetCurr + p + pos = int(layerK[q]) + + newNodes[i][pos] = idx + newPoints[offsetCurr + p, :] = points[q] + shiftK + p += 1 + # FIXME: Implement the other orders case _: raise ValueError(f'Extrusion not implemented for NGeo={order}') From c0fe8a921bd0e7c62f461a9f10a08b14a4f2fafe Mon Sep 17 00:00:00 2001 From: iagtirap Date: Wed, 25 Mar 2026 13:43:35 +0100 Subject: [PATCH 2/7] Bugfix on 2nd order --- pyhope/mesh/extrude/mesh_extrude.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 3dc3e938..53408e5d 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -510,7 +510,7 @@ def extrude_hexa(nodes: np.ndarray, for i in range(shifts.shape[0]-1): # Calculate offset for current layer indices - offsetCurr = i + offsetCurr = i*18 shiftCurr = shifts[i+1, :] shiftPrev = shifts[i , :] From 7677e580ef8a453ecb8b1a2be9b60b68d88cb357 Mon Sep 17 00:00:00 2001 From: iagtirap Date: Thu, 26 Mar 2026 11:45:23 +0100 Subject: [PATCH 3/7] Extend existing extrude function for 2D hexas to arbitrary order --- pyhope/mesh/extrude/mesh_extrude.py | 176 +++++++++++++--------------- 1 file changed, 83 insertions(+), 93 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 53408e5d..1e50f548 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -82,9 +82,6 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: elif not [cell_block for cell_block in mesh.cells if cell_block.type in gmshCellTypes.cellTypes2D]: hopout.error('Mesh contains no suitable surface cells for extrusion, exiting...') - if nGeo > 3: - hopout.error(f'nGeo = {nGeo} not supported for mesh extrusion') - hopout.info('Extruding surface to volume mesh') # Read in the mesh post-deformation flag @@ -382,6 +379,47 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: hopout.sep() return mesh +@cache +def quad_meshio_to_ij(order: int) -> npt.NDArray[np.int64]: + """Return meshio quad-node ordering mapped to tensor indices (i,j). + + The ordering is generated ring-by-ring: + corners, edge nodes, then inner rings, matching meshio high-order quad ordering. + """ + if order < 1: + raise ValueError(f'Quad ordering requires order >= 1, got {order}') + + ij: list[tuple[int, int]] = [] + + # Fill rings from outside to inside + for ring in range(order // 2 + 1): + start = ring + end = order - ring + + if start > end: + break + + # Odd order: final center point + if start == end: + ij.append((start, start)) + break + + # Ring corners + ij.extend(((start, start), (end, start), (end, end), (start, end))) + # Bottom edge (left -> right), excluding corners + ij.extend((i, start) for i in range(start + 1, end)) + # Right edge (bottom -> top), excluding corners + ij.extend((end, j) for j in range(start + 1, end)) + # Top edge (right -> left), excluding corners + ij.extend((i, end) for i in range(end - 1, start, -1)) + # Left edge (top -> bottom), excluding corners + ij.extend((start, j) for j in range(end - 1, start, -1)) + + expected = (order + 1) ** 2 + if len(ij) != expected: + raise RuntimeError(f'Invalid quad ordering length for order {order}: got {len(ij)}, expected {expected}') + + return np.asarray(ij, dtype=np.int64) def extrude_pris(nodes: np.ndarray, points: np.ndarray, @@ -488,98 +526,50 @@ def extrude_hexa(nodes: np.ndarray, nDOFsElem = (order+1)**3 newNodes = [np.empty((nDOFsElem, )) for s in range(len(shifts)-1)] - - match order: - case 1: - newPoints = np.empty((4*(shifts.shape[0]-1), 3)) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i *4 - shiftCurr = shifts[i+1, :] - - newNodes[i][ : 4] = nodes[ :4] if i == 0 else newNodes[i-1][ 4: 8] - newNodes[i][ 4: 8] = np.add(np.arange(4, dtype=np.int64), nPoints+offsetCurr) - newPoints[offsetCurr:offsetCurr+4, :] = points[:4] + shiftCurr - - case 2: - newPoints = np.empty((18*(shifts.shape[0]-1), 3)) - - # Append the bottom layer of the first element, then stack all the other elements - - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i*18 - shiftCurr = shifts[i+1, :] - shiftPrev = shifts[i , :] - - # Bottom/top corners - newNodes[i][ : 4] = nodes[ :4] if i == 0 else newNodes[i-1][ 4: 8] - newNodes[i][ 4: 8] = np.add(np.arange( 0, 4), nPoints+offsetCurr) - newPoints[offsetCurr :offsetCurr+ 4, :] = points[ :4] + shiftCurr - # Edges bottom/top - newNodes[i][ 8:12] = nodes[4:8] if i == 0 else newNodes[i-1][12:16] - newNodes[i][12:16] = np.add(np.arange( 4, 8), nPoints+offsetCurr) - newPoints[offsetCurr+ 4:offsetCurr+ 8, :] = points[4:8] + shiftCurr - # Edges upright - newNodes[i][16:20] = np.add(np.arange( 8, 12), nPoints+offsetCurr) - newPoints[offsetCurr+ 8:offsetCurr+12, :] = points[ :4] + 0.5*(shiftCurr+shiftPrev) - # Face centers - newNodes[i][20:24] = np.add(np.arange(12, 16), nPoints+offsetCurr) - newNodes[i][24:26] = np.array([int(nodes[8])if i == 0 else newNodes[i-1][25] , nPoints + 16 + offsetCurr]) - newPoints[offsetCurr+12:offsetCurr+13, :] = points[ 7] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+13:offsetCurr+14, :] = points[ 5] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+14:offsetCurr+15, :] = points[ 4] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+15:offsetCurr+16, :] = points[ 6] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+16:offsetCurr+17, :] = points[ 8] + shiftCurr - # Volume center - newNodes[i][26:27] = nPoints + 17 + offsetCurr - newPoints[offsetCurr+17:offsetCurr+18, :] = points[ 8] + 0.5*(shiftCurr+shiftPrev) - - case 3: - nFaceDOFs = (order+1)**2 - nNewDOFs = order*nFaceDOFs - newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) - - # Hardcoded meshio(quad16) -> meshio(hexa64) layer mapping for NGeo=3. - # Each entry is the hexa64 node position for face node q at layer k. - layer0 = np.array(( 0, 1, 2, 3, 8, 9, 10, 11, 12, 13, 15, 14, 48, 51, 50, 49), dtype=np.int64) - layer1 = np.array((24, 26, 28, 30, 40, 41, 36, 37, 44, 45, 35, 32, 56, 57, 58, 59), dtype=np.int64) - layer2 = np.array((25, 27, 29, 31, 43, 42, 39, 38, 47, 46, 34, 33, 60, 61, 62, 63), dtype=np.int64) - layer3 = np.array(( 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 23, 22, 52, 53, 54, 55), dtype=np.int64) - layers = (layer0, layer1, layer2, layer3) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i*nNewDOFs - shiftCurr = shifts[i+1, :] - shiftPrev = shifts[i , :] - - # Bottom layer: from original face for first layer, else from previous top layer. + + if order == 0: + raise ValueError(f'Extrusion not implemented for NGeo={order}') + else: + # Local imports ---------------------------------------- + from pyhope.mesh.mesh_common import LINMAP + # ------------------------------------------------------ + + nFaceDOFs = (order+1)**2 + nNewDOFs = order*nFaceDOFs + newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) + + # Generic meshio(quadN) -> meshio(hexaM) layer mapping for any NGeo>=1. + linmap = LINMAP(208, order) + q_to_ij = quad_meshio_to_ij(order) + layerPos = np.empty((order+1, nFaceDOFs), dtype=np.int64) + for q, (ii, jj) in enumerate(q_to_ij): + for k in range(order+1): + layerPos[k, q] = int(linmap[ii, jj, k]) + + # Append the bottom layer of the first element, then stack all the other elements + for i in range(shifts.shape[0]-1): + offsetCurr = i*nNewDOFs + shiftCurr = shifts[i+1, :] + shiftPrev = shifts[i , :] + + # Bottom layer + for q in range(nFaceDOFs): + idxBot = int(layerPos[0 , q]) + idxTop = int(layerPos[order, q]) + newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + + # New points for layers k=1...order + p = 0 + for k in range(1, order+1): + alpha = k/order + shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr for q in range(nFaceDOFs): - idxBot = int(layer0[q]) - idxTop = int(layer3[q]) - newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + idx = nPoints + offsetCurr + p + pos = int(layerPos[k, q]) - # New points for layers k=1..3 - p = 0 - for k in (1, 2, 3): - alpha = k/order - shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr - layerK = layers[k] - for q in range(nFaceDOFs): - idx = nPoints + offsetCurr + p - pos = int(layerK[q]) - - newNodes[i][pos] = idx - newPoints[offsetCurr + p, :] = points[q] + shiftK - p += 1 - - # FIXME: Implement the other orders - case _: - raise ValueError(f'Extrusion not implemented for NGeo={order}') + newNodes[i][pos] = idx + newPoints[offsetCurr + p, :] = points[q] + shiftK + p += 1 return newNodes, newPoints From f5a8de431385790106dcb7b3b5f5f2bf37b308cc Mon Sep 17 00:00:00 2001 From: iagtirap Date: Thu, 26 Mar 2026 14:03:48 +0100 Subject: [PATCH 4/7] Bugfix: problem pipeline (wrong committed file) --- pyhope/mesh/extrude/mesh_extrude.py | 87 ++++++++++++++--------------- 1 file changed, 41 insertions(+), 46 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 1e50f548..6bb87f89 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -36,12 +36,7 @@ # ---------------------------------------------------------------------------------------------------------------------------------- import meshio import numpy as np -# ---------------------------------------------------------------------------------------------------------------------------------- -# Typing libraries -# ---------------------------------------------------------------------------------------------------------------------------------- -import typing -if typing.TYPE_CHECKING: - import numpy.typing as npt +import numpy.typing as npt # ---------------------------------------------------------------------------------------------------------------------------------- # Local imports # ---------------------------------------------------------------------------------------------------------------------------------- @@ -526,50 +521,50 @@ def extrude_hexa(nodes: np.ndarray, nDOFsElem = (order+1)**3 newNodes = [np.empty((nDOFsElem, )) for s in range(len(shifts)-1)] - + if order == 0: raise ValueError(f'Extrusion not implemented for NGeo={order}') - else: - # Local imports ---------------------------------------- - from pyhope.mesh.mesh_common import LINMAP - # ------------------------------------------------------ - - nFaceDOFs = (order+1)**2 - nNewDOFs = order*nFaceDOFs - newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) - - # Generic meshio(quadN) -> meshio(hexaM) layer mapping for any NGeo>=1. - linmap = LINMAP(208, order) - q_to_ij = quad_meshio_to_ij(order) - layerPos = np.empty((order+1, nFaceDOFs), dtype=np.int64) - for q, (ii, jj) in enumerate(q_to_ij): - for k in range(order+1): - layerPos[k, q] = int(linmap[ii, jj, k]) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - offsetCurr = i*nNewDOFs - shiftCurr = shifts[i+1, :] - shiftPrev = shifts[i , :] - - # Bottom layer + + # Local imports ---------------------------------------- + from pyhope.mesh.mesh_common import LINMAP + # ------------------------------------------------------ + + nFaceDOFs = (order+1)**2 + nNewDOFs = order*nFaceDOFs + newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) + + # Generic meshio(quadN) -> meshio(hexaM) layer mapping for any NGeo>=1. + linmap = LINMAP(208, order) + q_to_ij = quad_meshio_to_ij(order) + layerPos = np.empty((order+1, nFaceDOFs), dtype=np.int64) + for q, (ii, jj) in enumerate(q_to_ij): + for k in range(order+1): + layerPos[k, q] = int(linmap[ii, jj, k]) + + # Append the bottom layer of the first element, then stack all the other elements + for i in range(shifts.shape[0]-1): + offsetCurr = i*nNewDOFs + shiftCurr = shifts[i+1, :] + shiftPrev = shifts[i , :] + + # Bottom layer + for q in range(nFaceDOFs): + idxBot = int(layerPos[0 , q]) + idxTop = int(layerPos[order, q]) + newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + + # New points for layers k=1...order + p = 0 + for k in range(1, order+1): + alpha = k/order + shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr for q in range(nFaceDOFs): - idxBot = int(layerPos[0 , q]) - idxTop = int(layerPos[order, q]) - newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] - - # New points for layers k=1...order - p = 0 - for k in range(1, order+1): - alpha = k/order - shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr - for q in range(nFaceDOFs): - idx = nPoints + offsetCurr + p - pos = int(layerPos[k, q]) + idx = nPoints + offsetCurr + p + pos = int(layerPos[k, q]) - newNodes[i][pos] = idx - newPoints[offsetCurr + p, :] = points[q] + shiftK - p += 1 + newNodes[i][pos] = idx + newPoints[offsetCurr + p, :] = points[q] + shiftK + p += 1 return newNodes, newPoints From eaa32875c2e1da6a116255c01654b20219f857af Mon Sep 17 00:00:00 2001 From: iagtirap Date: Fri, 27 Mar 2026 10:59:21 +0100 Subject: [PATCH 5/7] Extend existing extrude function for 2D prisms to arbitrary order --- pyhope/mesh/extrude/mesh_extrude.py | 165 ++++++++++++++-------------- 1 file changed, 81 insertions(+), 84 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 6bb87f89..7c9c1fe5 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -416,6 +416,44 @@ def quad_meshio_to_ij(order: int) -> npt.NDArray[np.int64]: return np.asarray(ij, dtype=np.int64) + +@cache +def tri_meshio_to_ij(order: int) -> npt.NDArray[np.int64]: + """Return meshio triangle-node ordering mapped to simplex indices (i,j).""" + if order < 1: + raise ValueError(f'Triangle ordering requires order >= 1, got {order}') + + ij: list[tuple[int, int]] = [] + + # Fill triangular rings from outside to inside. + # After one ring is emitted, the inner triangle has order reduced by 3. + def _fill_ring(n: int, i0: int, j0: int) -> None: + if n < 0: + return + if n == 0: + ij.append((i0, j0)) + return + + # Ring corners + ij.extend(((i0, j0), (i0+n, j0), (i0, j0+n))) + # Edge 0->1 (excluding corners) + ij.extend((i0+i, j0) for i in range(1, n)) + # Edge 1->2 (excluding corners) + ij.extend((i0+n-t, j0+t) for t in range(1, n)) + # Edge 2->0 (excluding corners) + ij.extend((i0, j0+j) for j in range(n-1, 0, -1)) + + _fill_ring(n-3, i0+1, j0+1) + + _fill_ring(order, 0, 0) + + expected = round((order + 1) * (order + 2) / 2) + if len(ij) != expected: + raise RuntimeError(f'Invalid triangle ordering length for order {order}: got {len(ij)}, expected {expected}') + + return np.asarray(ij, dtype=np.int64) + + def extrude_pris(nodes: np.ndarray, points: np.ndarray, shifts: np.ndarray, @@ -425,90 +463,49 @@ def extrude_pris(nodes: np.ndarray, nDOFsElem = round((order+1)**(3-1)*(order+2)/2.) newNodes = [np.empty((nDOFsElem, )) for s in range(len(shifts)-1)] - match order: - case 1: - newPoints = np.empty((3*(shifts.shape[0]-1), 3)) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i *3 - shiftCurr = shifts[i+1, :] - - newNodes[i][ : 3] = nodes[ :3] if i == 0 else newNodes[i-1][ 3: 6] - newNodes[i][ 3: 6] = np.add(np.arange(3, dtype=np.int64), nPoints+offsetCurr) - newPoints[offsetCurr:offsetCurr+3, :] = points[:3] + shiftCurr - - case 2: - newPoints = np.empty((12*(shifts.shape[0]-1), 3)) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i *12 - shiftCurr = shifts[i+1, :] - shiftPrev = shifts[i , :] - - # Bottom/top corners - newNodes[i][ : 3] = nodes[ :3] if i == 0 else newNodes[i-1][ 3: 6] - newNodes[i][ 3: 6] = np.add(np.arange( 0, 3), nPoints+offsetCurr) - newPoints[offsetCurr :offsetCurr+ 3, :] = points[ :3] + shiftCurr - # Edges bottom/top - newNodes[i][ 6: 9] = nodes[3:6] if i == 0 else newNodes[i-1][ 9:12] - newNodes[i][ 9:12] = np.add(np.arange( 3, 6), nPoints+offsetCurr) - newPoints[offsetCurr+ 3:offsetCurr+ 6, :] = points[3:6] + shiftCurr - # Edges upright - newNodes[i][12:15] = np.add(np.arange( 6, 9), nPoints+offsetCurr) - newPoints[offsetCurr+ 6:offsetCurr+ 9, :] = points[ :3] + 0.5*(shiftCurr+shiftPrev) - # Face centers - newNodes[i][15:18] = np.add(np.arange( 9, 12), nPoints+offsetCurr) - newPoints[offsetCurr+ 9:offsetCurr+10, :] = points[ 3] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+10:offsetCurr+11, :] = points[ 4] + 0.5*(shiftCurr+shiftPrev) - newPoints[offsetCurr+11:offsetCurr+12, :] = points[ 5] + 0.5*(shiftCurr+shiftPrev) - - case 3: - nFaceDOFs = round((order+1)*(order+2)/2.) - nNewDOFs = order*nFaceDOFs - newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) - - # Hardcoded meshio(triangle10) -> meshio(wedge40) layer mapping for NGeo=3. - # Each entry is the wedge40 node position for face node q at layer k. - layer0 = np.array((0, 1, 2, 6, 7, 8, 9, 10, 11, 37), dtype=np.int64) - layer1 = np.array((18, 20, 22, 24, 25, 28, 29, 32, 33, 38), dtype=np.int64) - layer2 = np.array((19, 21, 23, 27, 26, 31, 30, 35, 34, 39), dtype=np.int64) - layer3 = np.array((3, 4, 5, 12, 13, 14, 15, 16, 17, 36), dtype=np.int64) - layers = (layer0, layer1, layer2, layer3) - - # Append the bottom layer of the first element, then stack all the other elements - for i in range(shifts.shape[0]-1): - # Calculate offset for current layer indices - offsetCurr = i*nNewDOFs - shiftCurr = shifts[i+1, :] - shiftPrev = shifts[i , :] - - # Bottom layer: from original face for first layer, else from previous top layer. - for q in range(nFaceDOFs): - idxBot = int(layer0[q]) - idxTop = int(layer3[q]) - newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] - - # New points for layers k=1..3 - p = 0 - for k in (1, 2, 3): - alpha = k/order - shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr - layerK = layers[k] - for q in range(nFaceDOFs): - idx = nPoints + offsetCurr + p - pos = int(layerK[q]) - - newNodes[i][pos] = idx - newPoints[offsetCurr + p, :] = points[q] + shiftK - p += 1 - - # FIXME: Implement the other orders - case _: - raise ValueError(f'Extrusion not implemented for NGeo={order}') + if order == 0: + raise ValueError(f'Extrusion not implemented for NGeo={order}') + + # Local imports ---------------------------------------- + from pyhope.mesh.mesh_common import LINMAP + # ------------------------------------------------------ + + nFaceDOFs = round((order+1)*(order+2)/2.) + nNewDOFs = order*nFaceDOFs + newPoints = np.empty((nNewDOFs*(shifts.shape[0]-1), 3)) + + # Generic meshio(triaN) -> meshio(wedgeM) layer mapping for any supported NGeo>=1. + linmap = LINMAP(206, order) + q_to_ij = tri_meshio_to_ij(order) + layerPos = np.empty((order+1, nFaceDOFs), dtype=np.int64) + for q, (ii, jj) in enumerate(q_to_ij): + for k in range(order+1): + layerPos[k, q] = int(linmap[ii, jj, k]) + + # Append the bottom layer of the first element, then stack all the other elements + for i in range(shifts.shape[0]-1): + offsetCurr = i*nNewDOFs + shiftCurr = shifts[i+1, :] + shiftPrev = shifts[i , :] + + # Bottom layer + for q in range(nFaceDOFs): + idxBot = int(layerPos[0 , q]) + idxTop = int(layerPos[order, q]) + newNodes[i][idxBot] = nodes[q] if i == 0 else newNodes[i-1][idxTop] + + # New points for layers k=1..order + p = 0 + for k in range(1, order+1): + alpha = k/order + shiftK = (1.0-alpha)*shiftPrev + alpha*shiftCurr + for q in range(nFaceDOFs): + idx = nPoints + offsetCurr + p + pos = int(layerPos[k, q]) + + newNodes[i][pos] = idx + newPoints[offsetCurr + p, :] = points[q] + shiftK + p += 1 return newNodes, newPoints From 5ee806e4e64bd08c51c7a9cf2ee3b3045f269f60 Mon Sep 17 00:00:00 2001 From: Patrick Kopper Date: Fri, 27 Mar 2026 13:38:10 +0100 Subject: [PATCH 6/7] Update 2-15 tutorial to 4th order --- docs/user-guide/mesh-generators/external.md | 2 +- .../2-15-external_mesh_Gmsh_extrude.geo | 2 +- .../2-15-external_mesh_Gmsh_extrude_quad.geo | 2 +- .../2-15-external_mesh_Gmsh_extrude_tria.geo | 2 +- .../analyze.toml | 20 +++++++++---------- .../parameter.ini | 2 +- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/docs/user-guide/mesh-generators/external.md b/docs/user-guide/mesh-generators/external.md index dfa95e5c..ac79fc90 100644 --- a/docs/user-guide/mesh-generators/external.md +++ b/docs/user-guide/mesh-generators/external.md @@ -37,7 +37,7 @@ PyHOPE understands curved meshes of arbitrary order if the external mesh generat PyHOPE can utilize the agglomeration approach provided by Gmsh to curve linear meshes. In this case, the external mesh generator is used to generate a linear mesh, which is then curved through agglomeration in Gmsh before being read into PyHOPE. To use this approach, set the `NGeo` parameter in the [parameter file](../parameter-file.md) to the desired order larger than `1`. ## Mesh Extrusion -PyHOPE outputs purely three-dimensional meshes. However, some users might prefer to provide two-dimensional grids. PyHOPE can extrude meshes consisting of triangular and quadrilateral faces into prismatic and hexagonal meshes. This feature is activated using the parameter `MeshExtrude=T`. An extrusion template `MeshExtrudeTemplate` specifying the transformation of the one-dimensional `MeshExtrudeLength` extend into _N=_`MeshExtrudeElems` distinct elements must be provided. Additionally, since the far boundary of the extruded mesh cannot be associated with any of the existing boundaries, an additional boundary condition must be specified with the `MeshExtrudeBCIndex` parameter. Also note that this option is currently limited to `NGeo` \(\in\{1,2\}\). +PyHOPE outputs purely three-dimensional meshes. However, some users might prefer to provide two-dimensional grids. PyHOPE can extrude meshes consisting of triangular and quadrilateral faces into prismatic and hexagonal meshes. This feature is activated using the parameter `MeshExtrude=T`. An extrusion template `MeshExtrudeTemplate` specifying the transformation of the one-dimensional `MeshExtrudeLength` extend into _N=_`MeshExtrudeElems` distinct elements must be provided. Additionally, since the far boundary of the extruded mesh cannot be associated with any of the existing boundaries, an additional boundary condition must be specified with the `MeshExtrudeBCIndex` parameter. ## Generation of Hexahedral Meshes Many solvers require purely hexahedral meshes. PyHOPE can convert meshes consisting of tetrahedra, prisms and hexahedra into purely hexahedral meshes through a subdivision strategy. This feature is activated using the parameter `splitToHex=T`. Note, that pyramids cannot be decomposed to hexahedra in a straightforward way, thus this feature cannot be applied to meshes containing pyramids. Also note that this option is currently limited to `NGeo` \(\in\{1,2,4\}\). diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude.geo b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude.geo index 77ea6308..6d635955 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude.geo +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude.geo @@ -41,4 +41,4 @@ Physical Surface("BC_Front") = {1,2}; // Entire front (2D) domain Physical Surface("Zone_1") = {1}; // Zone 1 (all extruded quads) Physical Surface("Zone_2") = {2}; // Zone 2 (all extruded wedges) -Mesh.ElementOrder = 2; // Linear elements +Mesh.ElementOrder = 4; // Linear elements diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_quad.geo b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_quad.geo index be0205bb..b4636a9c 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_quad.geo +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_quad.geo @@ -23,4 +23,4 @@ Physical Curve("BC_Top") = {3}; // Put all curves in boundary conditions Physical Curve("BC_Left") = {4}; // Put all curves in Zone1 Physical Surface("BC_Front") = {1}; // Put the surface (and its elements) in Zone1 -Mesh.ElementOrder = 2; // Linear elements +Mesh.ElementOrder = 4; // Linear elements diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_tria.geo b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_tria.geo index 3aff0f17..450eecdd 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_tria.geo +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude_tria.geo @@ -19,4 +19,4 @@ Physical Curve("BC_Right" ) = {2}; // Put all curves in boundary conditions Physical Curve("BC_Left") = {3}; // Put all curves in Zone1 Physical Surface("BC_Front") = {1}; // Put the surface (and its elements) in Zone1 -Mesh.ElementOrder = 2; // Linear elements +Mesh.ElementOrder = 4; // Linear elements diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml b/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml index 8b6669e2..2cc14cea 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml @@ -1,23 +1,23 @@ [ElemInfo] min = 0.0 -max = 208.0 -mean = 65.16666666666667 -stddev = 73.51587130227956 +max = 600.0 +mean = 153.0 +stddev = 180.49807632092802 [GlobalNodeIDs] min = 1.0 -max = 84.0 -mean = 44.6 -stddev = 22.240503591420765 +max = 455.0 +mean = 236.125 +stddev = 124.90106501414095 [NodeCoords] min = 0.0 max = 2.0 -mean = 0.5888888888888288 -stddev = 0.47218954135280516 +mean = 0.5833333333332426 +stddev = 0.4347909956249531 [SideInfo] min = -22.0 max = 61.0 -mean = 6.848484848484849 -stddev = 13.52107493512574 +mean = 6.83030303030303 +stddev = 13.542581215144281 diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini b/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini index 13b3cb34..28ad6b62 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini @@ -10,7 +10,7 @@ OutputFormat = HDF5 ! Mesh output format (HDF5 ! MESH !=============================================================================== ! Mode = 3 ! 1 Cartesian 3 External -NGeo = 2 ! Polynomial order of the mapping +NGeo = 4 ! Polynomial order of the mapping FileName = ./2-15-external_mesh_Gmsh_extrude.geo !=============================================================================== ! From cbaa75a4229cebf1ffe80a2f5e8ae8d370fbd4cb Mon Sep 17 00:00:00 2001 From: Patrick Kopper Date: Fri, 27 Mar 2026 13:41:41 +0100 Subject: [PATCH 7/7] Add Michel, Justin to AUTHORS.md --- AUTHORS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/AUTHORS.md b/AUTHORS.md index 256303eb..b52b568e 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -18,8 +18,10 @@ are listed in alphabetical order: * Daniel Appel * Marcel Blind * Stephen Copplestone +* Justin Du Plessis * Patrick Kopper * Marius Kurz * Leon Teichröb * Felix Rodach * Anna Schwarz +* Michel Tirapelle