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 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/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 7a2cfa08..7c9c1fe5 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 # ---------------------------------------------------------------------------------------------------------------------------------- @@ -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 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 !=============================================================================== !