From e4b69c1238f0d9ec34d52960fad765f74bc9d58e Mon Sep 17 00:00:00 2001 From: Patrick Kopper Date: Mon, 16 Mar 2026 17:12:23 +0100 Subject: [PATCH 1/7] Extrusion: Add support for volume zones --- pyhope/basis/basis_jacobian.py | 2 +- pyhope/check/check_install.py | 2 +- pyhope/common/common.py | 2 +- pyhope/mesh/extrude/mesh_extrude.py | 104 ++++++++++++++---- pyhope/mesh/mesh.py | 4 +- pyhope/mesh/reader/reader_hopr.py | 2 +- pyhope/mesh/topology/mesh_topology.py | 16 ++- pyhope/meshio/meshio_convert.py | 2 +- pyproject.toml | 1 + .../2-15-external_mesh_Gmsh_extrude.geo | 2 + .../parameter.ini | 3 +- 11 files changed, 106 insertions(+), 34 deletions(-) diff --git a/pyhope/basis/basis_jacobian.py b/pyhope/basis/basis_jacobian.py index cce1aeda..0fffc29a 100644 --- a/pyhope/basis/basis_jacobian.py +++ b/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: diff --git a/pyhope/check/check_install.py b/pyhope/check/check_install.py index f3396522..5548fe36 100644 --- a/pyhope/check/check_install.py +++ b/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 \ diff --git a/pyhope/common/common.py b/pyhope/common/common.py index 57fe0b32..e9a0119a 100644 --- a/pyhope/common/common.py +++ b/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) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 56c5fe82..5e906db3 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/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,36 @@ 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) + # nFaces contains the existing number of [Triangle, Quad, Wedge, Hexahedron] + nFaces = np.zeros(4, 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) + # Prepare new cell blocks and new cell_sets elems_lst = {ftype: [] for ftype in faceType} csets_lst = {} @@ -129,8 +154,24 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: # nFace = (nGeo+1)*(nGeo+2)/2 nFace = 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 +186,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 +212,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 @@ -226,7 +272,7 @@ 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-1].name # noqa: E501 # BC: Set the BC for the bottom face appendBCSet(botFace, faceMap, nFace, nFaces, nodeToFace, faceType, @@ -234,12 +280,18 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: 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 + 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 +301,30 @@ 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(), [[], [], [], []]) 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 + faceVal = faceType.index(elemType) + name = zonFaces[0][1]['name'] + + 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, [[], [], [], []]) csets_lst[topName][faceVal].append(nFaces[faceVal]) nFaces[faceVal] += 1 diff --git a/pyhope/mesh/mesh.py b/pyhope/mesh/mesh.py index 14aac54b..a6eefd60 100644 --- a/pyhope/mesh/mesh.py +++ b/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') diff --git a/pyhope/mesh/reader/reader_hopr.py b/pyhope/mesh/reader/reader_hopr.py index 2dbdfd7c..98d22218 100644 --- a/pyhope/mesh/reader/reader_hopr.py +++ b/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: diff --git a/pyhope/mesh/topology/mesh_topology.py b/pyhope/mesh/topology/mesh_topology.py index b7b40e4a..1f1ff02c 100644 --- a/pyhope/mesh/topology/mesh_topology.py +++ b/pyhope/mesh/topology/mesh_topology.py @@ -684,7 +684,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)) + (2 if bcSide == 'zone' else 0) faceSet = frozenset(subFace) # Get candidate cset keys using the nodes in the face @@ -720,17 +722,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(), [[], [], [], []]) 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') diff --git a/pyhope/meshio/meshio_convert.py b/pyhope/meshio/meshio_convert.py index 9bb6e611..c2c45243 100644 --- a/pyhope/meshio/meshio_convert.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 23a4aa48..73c5d95b 100644 --- a/pyproject.toml +++ b/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 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 06c90580..77ea6308 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 @@ -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 diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini b/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini index 79974c0f..13b3cb34 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini +++ b/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 From 7b0dbf192597c2cef9fef151a96bac86071d4a15 Mon Sep 17 00:00:00 2001 From: Patrick Kopper Date: Mon, 16 Mar 2026 17:38:33 +0100 Subject: [PATCH 2/7] Extrusion: Add support for single-type meshes --- pyhope/mesh/extrude/mesh_extrude.py | 12 +++++++----- pyhope/mesh/topology/mesh_topology.py | 6 ++++-- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 5e906db3..29e99d7c 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -109,8 +109,6 @@ 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 contains the existing number of [Triangle, Quad, Wedge, Hexahedron] - nFaces = np.zeros(4, dtype=int) # Expected number of nodes faceNum = [ int((nGeo+1)*(nGeo+2)/2), int((nGeo+1)**2) ] @@ -139,6 +137,9 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: 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 = {} @@ -152,7 +153,7 @@ 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) @@ -284,6 +285,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: 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) @@ -308,7 +310,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: 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 @@ -324,7 +326,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: # 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, [[], [], [], []]) + csets_lst.setdefault(topName, [[] for _ in range(len(faceType))]) csets_lst[topName][faceVal].append(nFaces[faceVal]) nFaces[faceVal] += 1 diff --git a/pyhope/mesh/topology/mesh_topology.py b/pyhope/mesh/topology/mesh_topology.py index 1f1ff02c..4711c42c 100644 --- a/pyhope/mesh/topology/mesh_topology.py +++ b/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, @@ -686,7 +688,7 @@ def appendBCSet(subFace: np.ndarray, # 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)) + (2 if bcSide == 'zone' else 0) + 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 @@ -722,7 +724,7 @@ 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 From 965386cc347163b46a6ba08f3829d9e17034c890 Mon Sep 17 00:00:00 2001 From: Patrick Kopper Date: Mon, 16 Mar 2026 19:28:21 +0100 Subject: [PATCH 3/7] Extrusion: Fix extruded meshes with no explicit zone declaration --- pyhope/mesh/extrude/mesh_extrude.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 29e99d7c..159c86db 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -317,11 +317,12 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: elems_lst[faceType[faceVal]].append(np.array(subFace, dtype=int)) # Assign the new elements to the zone - faceVal = faceType.index(elemType) - name = zonFaces[0][1]['name'] - - csets_lst[name][faceVal].append(nFaces[faceVal]) - nFaces[faceVal] += 1 + 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 From 6f93b1e803f10ec9a5c89338595a6a18f63e33ec Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 18 Mar 2026 17:36:11 +0100 Subject: [PATCH 4/7] Added new tutorial, which combines a gmsh grid refined by region and 2D extruded mesh with multiple user-defined zones --- .../2-16-external_mesh_Gmsh_extrude_zones.geo | 133 ++++++++++++++++++ .../analyze.toml | 23 +++ .../parameter.ini | 43 ++++++ 3 files changed, 199 insertions(+) create mode 100644 tutorials/2-16-external_mesh_Gmsh_extrude_zones/2-16-external_mesh_Gmsh_extrude_zones.geo create mode 100644 tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml create mode 100644 tutorials/2-16-external_mesh_Gmsh_extrude_zones/parameter.ini diff --git a/tutorials/2-16-external_mesh_Gmsh_extrude_zones/2-16-external_mesh_Gmsh_extrude_zones.geo b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/2-16-external_mesh_Gmsh_extrude_zones.geo new file mode 100644 index 00000000..bcb46840 --- /dev/null +++ b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/2-16-external_mesh_Gmsh_extrude_zones.geo @@ -0,0 +1,133 @@ +//+ +SetFactory("OpenCASCADE"); + +//-------------------------------------- +// Points +//-------------------------------------- +Point(1) = {0, 0, 0}; +Point(2) = {10, 0, 0}; +Point(3) = {30, 0, 0}; +Point(4) = {0, 5, 0}; +Point(5) = {10, 5, 0}; +Point(6) = {30, 5, 0}; +Point(7) = {0, 10, 0}; +Point(8) = {10, 10, 0}; +Point(9) = {30, 10, 0}; +Point(10) = {0, 20, 0}; +Point(11) = {10, 20, 0}; +Point(12) = {30, 20, 0}; + +//-------------------------------------- +// Lines +//-------------------------------------- +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {4, 5}; +Line(4) = {5, 6}; +Line(5) = {7, 8}; +Line(6) = {8, 9}; +Line(7) = {10,11}; +Line(8) = {11,12}; + +Line(9) = {1, 4}; +Line(10) = {4, 7}; +Line(11) = {7,10}; + +Line(12) = {2, 5}; +Line(13) = {5, 8}; +Line(14) = {8,11}; + +Line(15) = {3, 6}; +Line(16) = {6, 9}; +Line(17) = {9,12}; + +// Surfaces +//-------------------------------------- +Line Loop(21) = {1, 12, -3, -9}; +Plane Surface(31) = {21}; + +Line Loop(22) = {3, 13, -5, -10}; +Plane Surface(32) = {22}; + +Line Loop(23) = {5, 14, -7, -11}; +Plane Surface(33) = {23}; + +Line Loop(24) = {2, 15, -4, -12}; +Plane Surface(34) = {24}; + +Line Loop(25) = {4, 16, -6, -13}; +Plane Surface(35) = {25}; + +Line Loop(26) = {6, 17, -8, -14}; +Plane Surface(36) = {26}; + +// Physical Curves +//-------------------------------------- +Physical Curve("ANODE", 25) = {10}; +Physical Curve("DIEANODE", 26) = {9, 11}; +Physical Curve("DIELECTRIC", 27) = {3, 12, 5, 14}; +Physical Curve("WALL", 28) = {7}; +Physical Curve("OUTLET", 29) = {8, 15, 16, 17}; +Physical Curve("SYMAXIS", 30) = {1, 2}; +Physical Curve("UNUSED", 31) = {4, 6, 13}; + +Physical Surface("ROTSYM") = {31,32,33,34,35,36}; + +Physical Surface("Zone_1") = {31}; +Physical Surface("Zone_2") = {33}; +Physical Surface("Zone_3") = {32,34,35,36}; + +//-------------------------------------- +// Mesh refinement: Fully refine surface 32 +//-------------------------------------- +Mesh.MeshSizeMin = 0.15; +Mesh.MeshSizeMax = 15.0; + +// Create a box field around Surface 32 (from x=0 to 10, y=5 to 10) +Field[1] = Box; +Field[1].VIn = 0.3; +Field[1].VOut = 0.5; +Field[1].XMin = 0; +Field[1].XMax = 10; +Field[1].YMin = 5; +Field[1].YMax = 10; +Field[1].ZMin = -1; +Field[1].ZMax = 1; + +// Optional: Use additional field to control transition further to the right +Field[2] = Distance; +Field[2].EdgesList = {3,5,13}; +Field[2].NumPointsPerCurve = 40; + +Field[3] = Threshold; +Field[3].InField = 2; +Field[3].SizeMin = 0.2; +Field[3].SizeMax = 15.0; +Field[3].DistMin = 2.0; +Field[3].DistMax = 40.0; + +// Combine the two using minimum +Field[4] = Min; +Field[4].FieldsList = {1, 3}; + +Background Field = 4; + +//-------------------------------------- +// Meshing options +//-------------------------------------- +Mesh.Algorithm = 8; +Mesh.RecombinationAlgorithm = 2; +Mesh.RecombineAll = 1; +Mesh 2; + +// Apply an elliptic smoother to the grid to have a more regular mesh: +Mesh.Smoothing = 100; + +//-------------------------------------- +// Output +//-------------------------------------- +Mesh.SaveAll = 1; +Mesh.Binary = 0; +Mesh.MshFileVersion = 4.1; + +Save "2-16-external_mesh_Gmsh_extrude_zones.msh"; diff --git a/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml new file mode 100644 index 00000000..29a6f68b --- /dev/null +++ b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml @@ -0,0 +1,23 @@ +[ElemInfo] +min = 0.0 +max = 40984.0 +mean = 11972.008068189212 +stddev = 12177.254106305043 + +[GlobalNodeIDs] +min = 1.0 +max = 10478.0 +mean = 5236.21286355651 +stddev = 2988.4474391070266 + +[NodeCoords] +min = 0.0 +max = 0.03 +mean = 0.006552551312148826 +stddev = 0.0068302865681914415 + +[SideInfo] +min = -20602.0 +max = 20607.0 +mean = 1051.047537250309 +stddev = 5312.915170664344 diff --git a/tutorials/2-16-external_mesh_Gmsh_extrude_zones/parameter.ini b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/parameter.ini new file mode 100644 index 00000000..b8fdbd3f --- /dev/null +++ b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/parameter.ini @@ -0,0 +1,43 @@ +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = 2-16-external_mesh_Gmsh_extrude_zones +DebugVisu = F +DebugMesh = F +! OutputFormat = GMSH +!=============================================================================== ! +! MESH +!=============================================================================== ! +FileName = 2-16-external_mesh_Gmsh_extrude_zones.geo +Mode = external +meshscale = 1e-3 + +!=============================================================================== ! +! MESH EXTRUSION +!=============================================================================== ! +MeshExtrude = T +MeshExtrudeBCIndexTop = 7 +MeshExtrudeBCIndexBot = 7 +MeshExtrudeLength = 1 +MeshExtrudeElems = 1 +MeshExtrudeDir = (/0., 0., 1./) + +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! +BoundaryName = ANODE +BoundaryType = (/5,0,0,0/) +BoundaryName = OUTLET +BoundaryType = (/10,0,0,0/) +BoundaryName = WALL +BoundaryType = (/5,0,0,0/) +BoundaryName = DIELECTRIC +BoundaryType = (/0,0,0,0/) +BoundaryName = SYMAXIS +BoundaryType = (/10,0,0,0/) +BoundaryName = DIEANODE +BoundaryType = (/5,0,0,0/) +BoundaryName = ROTSYM +BoundaryType = (/4,0,0,0/) +BoundaryName = UNUSED +BoundaryType = (/0,0,0,0/) From b6598b5c2e187a69b38f1eb61a5f11bc21aaa9ea Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Mon, 23 Mar 2026 20:52:47 +0100 Subject: [PATCH 5/7] Fix BC assignment for extruded boundaries --- pyhope/mesh/extrude/mesh_extrude.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 159c86db..7a2cfa08 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -242,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 @@ -273,7 +273,7 @@ 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[extrBCIndexTop-1].name # noqa: E501 + 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, From 8a4af05e9dd2208d2e5a15601753049eaac7bcb5 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Mon, 23 Mar 2026 21:01:53 +0100 Subject: [PATCH 6/7] Renamed tutorial 2-16-external_mesh_Gmsh_hybrid_splitToHex to 2-17-external_mesh_Gmsh_hybrid_splitToHex --- .../2-17-external_mesh_Gmsh_hybrid_splitToHex.msh} | 0 .../analyze.toml | 0 .../parameter.ini | 4 ++-- 3 files changed, 2 insertions(+), 2 deletions(-) rename tutorials/{2-16-external_mesh_Gmsh_hybrid_splitToHex/2-16-external_mesh_Gmsh_hybrid_splitToHex.msh => 2-17-external_mesh_Gmsh_hybrid_splitToHex/2-17-external_mesh_Gmsh_hybrid_splitToHex.msh} (100%) rename tutorials/{2-16-external_mesh_Gmsh_hybrid_splitToHex => 2-17-external_mesh_Gmsh_hybrid_splitToHex}/analyze.toml (100%) rename tutorials/{2-16-external_mesh_Gmsh_hybrid_splitToHex => 2-17-external_mesh_Gmsh_hybrid_splitToHex}/parameter.ini (91%) diff --git a/tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/2-16-external_mesh_Gmsh_hybrid_splitToHex.msh b/tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/2-17-external_mesh_Gmsh_hybrid_splitToHex.msh similarity index 100% rename from tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/2-16-external_mesh_Gmsh_hybrid_splitToHex.msh rename to tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/2-17-external_mesh_Gmsh_hybrid_splitToHex.msh diff --git a/tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/analyze.toml b/tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/analyze.toml similarity index 100% rename from tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/analyze.toml rename to tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/analyze.toml diff --git a/tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini b/tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini similarity index 91% rename from tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini rename to tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini index 2184f964..8dd0b0f1 100644 --- a/tutorials/2-16-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini +++ b/tutorials/2-17-external_mesh_Gmsh_hybrid_splitToHex/parameter.ini @@ -1,7 +1,7 @@ !=============================================================================== ! ! OUTPUT !=============================================================================== ! -ProjectName = 2-16-external_mesh_Gmsh_hybrid_splitToHex ! Name of output files +ProjectName = 2-17-external_mesh_Gmsh_hybrid_splitToHex ! Name of output files DebugVisu = F ! Launch the GMSH GUI to visualize the mesh DebugMesh = T ! Mesh output debug mesh in VTK format OutputFormat = HDF5 ! Mesh output format (HDF5 VTK) @@ -11,7 +11,7 @@ OutputFormat = HDF5 ! Mesh output format (HDF5 !=============================================================================== ! Mode = 3 ! 1 Cartesian 3 External NGeo = 1 ! Polynomial order of the mapping -FileName = ./2-16-external_mesh_Gmsh_hybrid_splitToHex.msh +FileName = ./2-17-external_mesh_Gmsh_hybrid_splitToHex.msh doSplitToHex = T ! Split simplex elements into hexahedral elements doSplitToHexZ = F ! Split hexahedral elements into h-refined elements From 60221dc6a3e0a6367206de542e25721c025e1471 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Mon, 23 Mar 2026 21:09:52 +0100 Subject: [PATCH 7/7] Added tolerance to extrusion tutorial, which uses gmsh refinement and smoothing --- .../2-16-external_mesh_Gmsh_extrude_zones/analyze.toml | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml index 29a6f68b..ec573e81 100644 --- a/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml +++ b/tutorials/2-16-external_mesh_Gmsh_extrude_zones/analyze.toml @@ -3,6 +3,7 @@ min = 0.0 max = 40984.0 mean = 11972.008068189212 stddev = 12177.254106305043 +_atol = 400.0 # Rounding issues [GlobalNodeIDs] min = 1.0 @@ -13,11 +14,12 @@ stddev = 2988.4474391070266 [NodeCoords] min = 0.0 max = 0.03 -mean = 0.006552551312148826 -stddev = 0.0068302865681914415 +mean = 0.006552555028766583 +stddev = 0.00683029203037687 [SideInfo] min = -20602.0 max = 20607.0 -mean = 1051.047537250309 -stddev = 5312.915170664344 +mean = 1051.0808705836423 +stddev = 5312.908617044128 +_atol = 200.0 # Rounding issues