[Feature]: High-order extrusion

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

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\}\).

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
 # ----------------------------------------------------------------------------------------------------------------------------------
@@ -82,9 +77,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 > 2:
-        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 +374,85 @@ 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)
+
+
+@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,
@@ -392,50 +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)
-
-        # 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
 
@@ -449,57 +519,49 @@ 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)
-
-        # 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 = (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):
+                idx = nPoints + offsetCurr + p
+                pos = int(layerPos[k, q])
+
+                newNodes[i][pos] = idx
+                newPoints[offsetCurr + p, :] = points[q] + shiftK
+                p += 1
 
     return newNodes, newPoints
 

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

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

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

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

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
 
 !=============================================================================== !