Extrusion: Add support for volume zones

pyhope/basis/basis_jacobian.py

@@ -207,7 +207,7 @@ def CheckJacobians() -> None:
             mapLin    = np.array(tuple(mapLin[np.int64(i)] for i in range(len(mapLin))))
             linCache[elemType] = mapLin
         # Only hexahedrons supported for specific nGeo
-        except ValueError:
+        except ValueError:  # noqa: PERF203
             pass
 
     for elem in elems:

pyhope/check/check_install.py

@@ -89,7 +89,7 @@ def _make_request(url : str,
         while True:
             try:
                 return urllib.request.urlopen(req)
-            except HTTPError as e:
+            except HTTPError as e:  # noqa: PERF203
                 # Check for rate-limiting error
                 if  e.code == 403                        \
                 and 'X-RateLimit-Remaining' in e.headers \

pyhope/common/common.py

@@ -224,7 +224,7 @@ def find_indices(seq: Union[list, npt.NDArray], item) -> tuple[int, ...]:
     while True:
         try:
             loc = cast(list, seq).index(item, start_at+1)
-        except ValueError:
+        except ValueError:  # noqa: PERF203
             break
         else:
             locs.append(loc)

pyhope/mesh/extrude/mesh_extrude.py

@@ -89,8 +89,9 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
 
     # Read in the mesh post-deformation flag
     hopout.sep()
-    extrTemplate = GetStr( 'MeshExtrudeTemplate')
-    extrBCIndex  = GetInt( 'MeshExtrudeBCIndex')
+    extrTemplate   = GetStr('MeshExtrudeTemplate')
+    extrBCIndexTop = GetInt('MeshExtrudeBCIndexTop') - 1
+    extrBCIndexBot = GetInt('MeshExtrudeBCIndexBot') - 1
 
     # Continue with extrusion
     hopout.sep()
@@ -108,12 +109,37 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
     # Get base key to distinguish between linear and high-order elements
     ho_key    = 100 if nGeo == 1 else 200
     nPoints   = len(pointl)
-    nFaces    = np.zeros(2, dtype=int)
 
     # Expected number of nodes
     faceNum   = [ int((nGeo+1)*(nGeo+2)/2), int((nGeo+1)**2) ]
     faceType  = [f'triangle{"" if nGeo == 1 else faceNum[0]}', f'quad{"" if nGeo == 1 else faceNum[1]}']
 
+    # For zones, we need to append the expected extruded elements
+    for cblock in cell_sets.values():
+        # Each set_blocks is a list of arrays, one entry per cell block
+        for blockID in range(len(cblock)):
+            etype = elems_old[blockID].type
+            if etype[:4] not in ('tria', 'quad'):
+                continue
+
+            elemNum = ho_key + (8 if cast(str, etype).startswith('quad') else 6)
+            # Obtain the element type
+            elemType = elemTypeClass.inam[elemNum]
+            if len(elemType) > 1:
+                elemType  = elemType[0].rstrip(digits)
+                elemDOFs  = NDOFperElemType(elemType, mesh_vars.nGeo)
+                elemType += str(elemDOFs)
+            else:
+                elemType  = elemType[0]
+                elemDOFs  = NDOFperElemType(elemType, mesh_vars.nGeo)
+
+            if elemType not in faceType:
+                faceType.append(elemType)
+                faceNum .append(elemDOFs)
+
+    # nFaces contains the existing number of [Triangle, Quad, Wedge, Hexahedron]
+    nFaces    = np.zeros(len(faceType), dtype=int)
+
     # Prepare new cell blocks and new cell_sets
     elems_lst = {ftype: [] for ftype in faceType}
     csets_lst = {}
@@ -127,10 +153,26 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
     # INFO: We reduce the new faces to first-order. Yes, this breaks direkt meshio output. But we are not using this anyways.
     #       If you want to use mesh.write() for debug purposes, comment out the BC face creation.
     # nFace = (nGeo+1)*(nGeo+2)/2
-    nFace = 3
+    nFace: Final[int] = 3
+
+    # Create the element sets
+    meshcells = tuple((k, v) for k, v in mesh.cell_sets_dict.items() if any(key.startswith('tria') for key in v)
+                                                                     or any(key.startswith('quad') for key in v))
+
+    # Take the correct BC index for the bottom BC
+    meshcells = tuple(s for s in meshcells if s[0].lower() == mesh_vars.bcs[extrBCIndexBot].name)
+
+    match len(meshcells):
+        case 0:
+            hopout.error('Could not find boundary condition for extrusion, exiting...')
+        case 1:
+            pass
+        case _:
+            hopout.error('Found more than one boundary condition for extrusion, exiting...')
 
     # Convert the (1D, 2D) boundary cell set into a dictionary
     csets_old = {}
+    zsets_old = {}
     for cname, cblock in cell_sets.items():
         # Each set_blocks is a list of arrays, one entry per cell block
         for blockID, block in enumerate(cblock):
@@ -145,24 +187,23 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
             # Determine how many corner nodes to keep
             nCorners = 2 if 'line' in etype else (3 if 'tria' in etype else 4)
 
+            # Filter the zone 2D faces
+            if nCorners > 2 and cname.lower() != mesh_vars.bcs[extrBCIndexBot].name:
+                # Sort them as a set for membership checks
+                for face in block:
+                    # Slice to only include corners for the search dictionary
+                    nodes = mesh.cells[blockID].data[face][:nCorners]
+                    zsets_old.setdefault(frozenset(nodes), []).append(cname)
+
+                # Do not add them to the boundary conditions
+                continue
+
             # Sort them as a set for membership checks
             for face in block:
                 # Slice to only include corners for the search dictionary
                 nodes = mesh.cells[blockID].data[face][:nCorners]
                 csets_old.setdefault(frozenset(nodes), []).append(cname)
 
-    # Create the element sets
-    meshcells = tuple((k, v) for k, v in mesh.cell_sets_dict.items() if any(key.startswith('tria') for key in v)
-                                                                     or any(key.startswith('quad') for key in v))
-
-    match len(meshcells):
-        case 0:
-            hopout.error('Could not found boundary condition for extrusion, exiting...')
-        case 1:
-            pass
-        case _:
-            hopout.error('Found more than one boundary condition for extrusion, exiting...')
-
     nTotalElems = sum(cdata.shape[0] for _, zdata in meshcells for cdata in cast(dict, zdata).values())
     bar = ProgressBar(value=nTotalElems, title='│             Processing Elements', length=33, threshold=1000)
 
@@ -172,6 +213,12 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
         for node in subFace:
             nodeToFace[node].add(subFace)
 
+    # Build an inverted index to map each node to all zone keys (from zsets_old) that contain it
+    nodeToZone = defaultdict(set)
+    for subFace in zsets_old:
+        for node in subFace:
+            nodeToZone[node].add(subFace)
+
     # We need to unwrap meshcells for each zone, i.e. each 2D boundary condition
     for meshcell in meshcells:
         _    , mdict = meshcell
@@ -195,11 +242,11 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
             # Obtain the element type
             elemType = elemTypeClass.inam[elemNum]
             if len(elemType) > 1:
-                elemType  = elemType[0].rstrip(digits)
+                elemType  = str(elemType[0]).rstrip(digits)
                 elemDOFs  = NDOFperElemType(elemType, mesh_vars.nGeo)
                 elemType += str(elemDOFs)
             else:
-                elemType  = elemType[0]
+                elemType  = str(elemType[0])
                 elemDOFs  = NDOFperElemType(elemType, mesh_vars.nGeo)
 
             # Face block: Iterate over each element
@@ -226,20 +273,27 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
                 # > We know this is the first face and 1/5 face
                 botIdx , botFace = 0, subFaces[0]
                 topIdx           = 1 if mtype.startswith('tria') else 5
-                topFace, topName = [np.array(extElems[-1])[face] for face in faces(nGeo)][topIdx], mesh_vars.bcs[extrBCIndex-1].name
+                topFace, topName = [np.array(extElems[-1])[face] for face in faces(nGeo)][topIdx], mesh_vars.bcs[extrBCIndexTop].name  # noqa: E501
 
                 # BC: Set the BC for the bottom face
                 appendBCSet(botFace, faceMap, nFace, nFaces, nodeToFace, faceType,
                             csets_old  = csets_old      , csets_lst    = csets_lst, elems_lst  = elems_lst,  # noqa: E251, E271
                             bcFaces    = bcFaces        , bcFaceIdx    = botIdx   , bcSide     = 'bottom',   # noqa: E251, E271
                             requireDim = lambda n: n > 2, requireMatch = False    , allowMulti = False)      # noqa: E251, E271
 
+                # ZONE: Set the ZONE for the bottom element
+                appendBCSet(botFace, faceMap, nFace, nFaces, nodeToZone, faceType,
+                            csets_old  = zsets_old      , csets_lst    = csets_lst, elems_lst  = elems_lst,  # noqa: E251, E271
+                            bcFaces    = bcFaces        , bcFaceIdx    = botIdx   , bcSide     = 'zone',     # noqa: E251, E271
+                            elemType   = elemType,                                                           # noqa: E251, E271
+                            requireDim = lambda n: n > 2, requireMatch = False    , allowMulti = False)      # noqa: E251, E271
+
                 # BC: Next, iterate over the 1D (side faces)
                 for iFace, subFace in enumerate(subFaces[1::]):
                     appendBCSet(subFace, faceMap, nFace, nFaces, nodeToFace, faceType,
-                                csets_old  = csets_old, csets_lst    = csets_lst, elems_lst  = elems_lst,  # noqa: E251, E271
-                                bcFaces    = bcFaces  , bcFaceIdx    = iFace+1  , bcSide     = 'side',     # noqa: E251, E271
-                                requireDim = 2        , requireMatch = False    , allowMulti = False)      # noqa: E251, E271
+                                csets_old  = csets_old, csets_lst    = csets_lst, elems_lst  = elems_lst,    # noqa: E251, E271
+                                bcFaces    = bcFaces  , bcFaceIdx    = iFace+1  , bcSide     = 'side',       # noqa: E251, E271
+                                requireDim = 2        , requireMatch = False    , allowMulti = False)        # noqa: E251, E271
 
                 for extElem in extElems[1:]:
                     # Overwrite the element with the new indices
@@ -249,22 +303,31 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
                     # Create the new faces
                     subFaces = tuple(np.array(extElem)[face] for face in faces(nGeo))
                     sidFaces = tuple((i, s) for i, s in enumerate(bcFaces) if ('side' in s and s['side'] == 'side'))
+                    zonFaces = tuple((i, s) for i, s in enumerate(bcFaces) if ('side' in s and s['side'] == 'zone'))
 
                     for iFace, sidFace in sidFaces:
                         subFace = subFaces[iFace]
                         faceVal = faceMap(0) if len(subFace) == nFace else faceMap(1)
 
                         name = sidFace['name']
-                        csets_lst.setdefault(name.strip(), [[], []])
+                        csets_lst.setdefault(name.strip(), [[] for _ in range(len(faceType))])
                         csets_lst[name][faceVal].append(nFaces[faceVal])
 
                         nFaces[faceVal] += 1
                         elems_lst[faceType[faceVal]].append(np.array(subFace, dtype=int))
 
+                    # Assign the new elements to the zone
+                    if zonFaces:
+                        faceVal = faceType.index(elemType)
+                        name    = zonFaces[0][1]['name']
+                        # Append the volume zones and increment
+                        csets_lst[name][faceVal].append(nFaces[faceVal])
+                        nFaces[faceVal] += 1
+
                 # BC: We should have one face left, assign the bottom BC
                 # > We need to hardcode this since we might have internal faces
-                faceVal    = faceMap(0) if len(topFace) == nFace else faceMap(1)
-                csets_lst.setdefault(topName, [[], []])
+                faceVal = faceMap(0) if len(topFace) == nFace else faceMap(1)
+                csets_lst.setdefault(topName, [[] for _ in range(len(faceType))])
                 csets_lst[topName][faceVal].append(nFaces[faceVal])
 
                 nFaces[faceVal] += 1

pyhope/mesh/mesh.py

@@ -130,7 +130,9 @@ def DefineMesh() -> None:
     CreateReal(      'MeshExtrudeLength',    default=1.0,            help='Mesh extrusion length')
     CreateRealArray( 'MeshExtrudeDir',  3,   default='(/0.,0.,1./)', help='Mesh extrusion direction')
     CreateInt(       'MeshExtrudeElems',     default=1  ,            help='Mesh extrusion number of element')
-    CreateInt(       'MeshExtrudeBCIndex',                           help='Mesh extrusion boundary index')
+    CreateInt(       'MeshExtrudeBCIndexBot',                        help='Mesh extrusion boundary index')
+    CreateInt(       'MeshExtrudeBCIndexTop',                        help='Mesh extrusion boundary index')
+
     # Edge connectivity
     CreateSection('Finite Element Method (FEM) Connectivity')
     CreateLogical(   'doFEMConnect',         default=False,          help='Generate finite element method (FEM) connectivity')

pyhope/mesh/reader/reader_hopr.py

@@ -170,7 +170,7 @@ def ReadHOPR(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh:
                     mapLin    = np.array(tuple(mapLin[np.int64(i)] for i in range(len(mapLin))))
                     linCache[elemType] = mapLin
                 # Only hexahedrons supported for specific nGeo
-                except ValueError:
+                except ValueError:  # noqa: PERF203
                     pass
 
             with alive_bar(len(elemInfo), title='│             Processing Elements', length=33) as bar:

pyhope/mesh/topology/mesh_topology.py

@@ -675,6 +675,8 @@ def appendBCSet(subFace:      np.ndarray,
                 bcFaces:      Optional[list] = None,
                 bcFaceIdx:    Optional[int]  = None,
                 bcSide:       Optional[str]  = None,
+                # Optional zone element
+                elemType:     Optional[str]  = None,
                 # Optional checks
                 requireDim:   Optional[Union[Callable[[int], bool], int]] = None,
                 requireMatch: bool = False,
@@ -684,7 +686,9 @@ def appendBCSet(subFace:      np.ndarray,
     import pyhope.output.output as hopout
     # ------------------------------------------------------
 
-    faceVal = faceMap(0) if len(subFace) == nFace else faceMap(1)
+    # We need to distinguish between BC (2D) and ZONE (3D) elements. The actual size might not match because we are getting the
+    # pre-extrusion types. Manually add +2 if we want to create a zone
+    faceVal = (faceMap(0) if len(subFace) == nFace else faceMap(1)) if elemType is None else faceType.index(elemType)
     faceSet = frozenset(subFace)
 
     # Get candidate cset keys using the nodes in the face
@@ -720,17 +724,21 @@ def appendBCSet(subFace:      np.ndarray,
 
             # Update csets_lst for each name in the list.
             for name in names:
-                csets_lst.setdefault(name.strip(), [[], []])
+                csets_lst.setdefault(name.strip(), [[] for _ in range(len(faceType))])
                 csets_lst[name][faceVal].append(nFaces[faceVal])
+                nFaces[faceVal] += 1
 
-                # Store the 1D faces
+                # Store the (original) 1D/2D faces
                 if bcFaces is not None and bcFaceIdx is not None and bcSide is not None:
                     bcFaces[bcFaceIdx] = {'name': name.strip(),
                                           'side': bcSide.strip(),
                                          }
 
-                nFaces[faceVal] += 1
-                elems_lst[faceType[faceVal]].append(np.array(subFace, dtype=int))
+                match faceType[faceVal][:4]:
+                    # Append the 2D (newly created) faces
+                    # > For zones, the elements are already created in the calling function
+                    case 'tria' | 'quad':
+                        elems_lst[faceType[faceVal]].append(np.array(subFace, dtype=int))
 
     if requireMatch and not common_match:
         raise ValueError('Unable to identify BC for face')

pyhope/meshio/meshio_convert.py

@@ -308,6 +308,6 @@ def MeshioGmshOrderingPatch() -> None:
         try:
             mod = importlib.import_module(mod_name)
             mod._meshio_to_gmsh_order = NodeOrdering().ordering_meshio_to_gmsh
-        except Exception:
+        except Exception:  # noqa: PERF203
             # If assignment fails, pass
             pass

pyproject.toml

@@ -169,6 +169,7 @@ extend-select = [
   "SIM",      # flake8-simplify
   "SLOT",     # flake8-slots
   "TID",      # flake8-tidy-imports
+  "PERF",     # Perflint
   "FURB",     # refurb
   "I",        # isort
   "NPY",      # NumPy-specific rules

tutorials/2-15-external_mesh_Gmsh_extrude/2-15-external_mesh_Gmsh_extrude.geo

@@ -38,5 +38,7 @@ Physical Curve("BC_Left")   = {4};              // Outer boundary
 Physical Curve("BC_Right")  = {6};              // Triangle outer boundary
 
 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

tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini

@@ -20,7 +20,8 @@ MeshExtrude         = T                              ! Enables mesh extrusion
 MeshExtrudeTemplate = linear                         ! Mesh extrusion template
 MeshExtrudeLength   = 1.0                            ! Mesh extrusion length
 MeshExtrudeElems    = 3                              ! Mesh extrusion number of elements
-MeshExtrudeBCIndex  = 6                              ! Mesh extrusion boundary index
+MeshExtrudeBCIndexBot = 5                            ! Mesh extrusion boundary index (bottom)
+MeshExtrudeBCIndexTop = 6                            ! Mesh extrusion boundary index (top)
 
 !=============================================================================== !
 ! BOUNDARY CONDITIONS