From 6a0101aaed49b437756dc1bc01c256751a37d397 Mon Sep 17 00:00:00 2001 From: Tobias Ott Date: Fri, 27 Mar 2026 10:34:12 +0100 Subject: [PATCH 1/2] Fix order of local edges --- pyhope/io/formats/meshio.py | 24 ++++++++++++------------ pyhope/mesh/mesh_builtin.py | 14 ++++++++------ pyhope/mesh/mesh_common.py | 32 ++++++++++++++++---------------- 3 files changed, 36 insertions(+), 34 deletions(-) diff --git a/pyhope/io/formats/meshio.py b/pyhope/io/formats/meshio.py index e8c8779b..e1db13fe 100644 --- a/pyhope/io/formats/meshio.py +++ b/pyhope/io/formats/meshio.py @@ -125,29 +125,29 @@ def edgePointMESHIO(start: int, end: int, edge: int, node: int) -> npt.NDArray: """ match edge: case 0: - return np.array((node , start , start ), dtype=int) + return np.array((node , start , start ), dtype=int) case 1: - return np.array((end , node , start ), dtype=int) + return np.array((end , node , start ), dtype=int) case 2: - return np.array((end+start-node, end , start ), dtype=int) + return np.array((end+start-node, end , start ), dtype=int) case 3: - return np.array((start , node , start ), dtype=int) + return np.array((start , end+start-node, start ), dtype=int) case 4: - return np.array((node , start , end ), dtype=int) + return np.array((start , start , node ), dtype=int) case 5: - return np.array((end , node , end ), dtype=int) + return np.array((end , start , node ), dtype=int) case 6: - return np.array((end+start-node, end , end ), dtype=int) + return np.array((end , end , node ), dtype=int) case 7: - return np.array((start , node , end ), dtype=int) + return np.array((start , end , node ), dtype=int) case 8: - return np.array((start , start , node ), dtype=int) + return np.array((node , start , end ), dtype=int) case 9: - return np.array((end , start , node ), dtype=int) + return np.array((end , node , end ), dtype=int) case 10: - return np.array((end , end , node ), dtype=int) + return np.array((end+start-node, end , end ), dtype=int) case 11: - return np.array((start , end , node ), dtype=int) + return np.array((start , end+start-node, end ), dtype=int) case _: raise ValueError(f'Invalid edge index: {edge}.') diff --git a/pyhope/mesh/mesh_builtin.py b/pyhope/mesh/mesh_builtin.py index 9e2b904c..ebd6587e 100644 --- a/pyhope/mesh/mesh_builtin.py +++ b/pyhope/mesh/mesh_builtin.py @@ -154,13 +154,15 @@ def MeshCartesian() -> meshio.Mesh: # Connect the corner points e = [None for _ in range(12)] - # First, the plane surface - for i in range(2): - for j in range(4): - e[j + i*4] = gmsh.model.geo.addLine(p[j + i*4], p[(j+1) % 4 + i*4]) - # Then, the connection + # Bottom face edges (0-3) for j in range(4): - e[j+8] = gmsh.model.geo.addLine(p[j], p[j+4]) + e[j] = gmsh.model.geo.addLine(p[j], p[(j+1) % 4]) + # Upright edges (4-7) + for j in range(4): + e[j+4] = gmsh.model.geo.addLine(p[j], p[j+4]) + # Top face edges (8-11) + for j in range(4): + e[j+8] = gmsh.model.geo.addLine(p[j+4], p[(j+1) % 4 + 4]) # Get dimensions of domain gmsh.model.geo.synchronize() diff --git a/pyhope/mesh/mesh_common.py b/pyhope/mesh/mesh_common.py index c08ce4fe..32d09406 100644 --- a/pyhope/mesh/mesh_common.py +++ b/pyhope/mesh/mesh_common.py @@ -106,9 +106,9 @@ def edge_to_dir(edge: int, elemType: Union[int, str]) -> int: # Pyramid # Wedge / Prism # Hexahedron - 8: { 0: eps, 2: eps, 4: eps, 6: eps, # Direction 0 - 1: 1., 3: 1., 5: 1., 7: 1., # Direction 1 - 8: 2., 9: 2., 10: 2., 11: 2.} # Direction 2 + 8: { 0: eps, 2: eps, 8: eps, 10: eps, # Direction 0 + 1: 1., 3: 1., 9: 1., 11: 1., # Direction 1 + 4: 2., 5: 2., 6: 2., 7: 2.} # Direction 2 } if isinstance(elemType, str): @@ -130,18 +130,18 @@ def edge_to_corner(edge: int, elemType: Union[int, str], dtype=np.int32) -> npt. """ GMSH: Get points on edges """ edge_map = { # Tetrahedron - 4: ( (0, 1), (1, 2), (2, 1), (0, 3), + 4: ( (0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3) ), # Pyramid 5: ( (0, 1), (1, 2), (2, 3), (3, 0), - (0, 4), (1, 5), (2, 4), (3, 4) ), + (0, 4), (1, 4), (2, 4), (3, 4) ), # Wedge / Prism 6: ( (0, 1), (1, 2), (2, 0), (0, 3), - (2, 3), (3, 4), (4, 5), (5, 4) ), + (1, 4), (2, 5), (3, 4), (4, 5), (5, 3) ), # Hexahedron 8: ( (0, 1), (1, 2), (2, 3), (3, 0), # Bottom edges - (4, 5), (5, 6), (6, 7), (7, 4), # Top edges - (0, 4), (1, 5), (2, 6), (3, 7) ), # Upright edges + (0, 4), (1, 5), (2, 6), (3, 7), # Upright edges + (4, 5), (5, 6), (6, 7), (7, 4) ), # Top edges } if isinstance(elemType, str): @@ -167,8 +167,8 @@ def edge_to_sign(edge: int, elemType: Union[int, str], dtype=np.float64) -> npt. # Wedge / Prism # Hexahedron 8: ( -1., -1., 1., 1., # Bottom edges - -1., -1., 1., 1., # Top edges - -1., -1., -1., -1.), # Upright edges + -1., -1., -1., -1., # Upright edges + -1., -1., 1., 1.), # Top edges } if isinstance(elemType, str): elemType = elemTypeClass.name[elemType] @@ -192,12 +192,12 @@ def face_to_edge(face: str, elemType: Union[str, int], dtype=np.int32) -> npt.ND # Pyramid # Wedge / Prism # Hexahedron - 8: { 'z-': np.array(( 0, 1, 2, 3), dtype=dtype), - 'y-': np.array(( 0, 9, -4, -8), dtype=dtype), - 'x+': np.array(( 1, 10, -5, -9), dtype=dtype), - 'y+': np.array(( -2, 10, 6, -11), dtype=dtype), - 'x-': np.array(( 8, -7, -11, 3), dtype=dtype), - 'z+': np.array(( 4, 5, 6, 7), dtype=dtype)} + 8: { 'z-': np.array(( 0, 1, 2, 3), dtype=dtype), + 'y-': np.array(( 0, 5, -8, -4), dtype=dtype), + 'x+': np.array(( 1, 6, -9, -5), dtype=dtype), + 'y+': np.array(( 6, 10, -7, -2), dtype=dtype), + 'x-': np.array(( 4, -11, -7, 3), dtype=dtype), + 'z+': np.array(( 8, 9, 10, 11), dtype=dtype)} } if isinstance(elemType, str): From db080874db1a8715d10cf75c53e7456c72d61a26 Mon Sep 17 00:00:00 2001 From: Tobias Ott Date: Fri, 27 Mar 2026 11:06:05 +0100 Subject: [PATCH 2/2] Fix EdgeID output in EdgeInfo --- pyhope/mesh/fem/fem.py | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/pyhope/mesh/fem/fem.py b/pyhope/mesh/fem/fem.py index cf83235a..e9e16778 100644 --- a/pyhope/mesh/fem/fem.py +++ b/pyhope/mesh/fem/fem.py @@ -410,19 +410,21 @@ def getFEMInfo(nodeInfo: npt.NDArray) -> tuple[npt.NDArray, # FEMElemInfo # Edge connectivity info --------------------------------------------------- # > Build flat arrays of all raw edge occurrences # > > Same order as the elements - nOccEdge = sum(len(cast(dict, elem.edgeInfo)) for elem in elems) - occEdgeID = np.empty(nOccEdge, dtype=np.int32) # FEMEdgeID per occurrence - occElemID = np.empty(nOccEdge, dtype=np.int32) # Element index per occurrence - occLocEdge = np.empty(nOccEdge, dtype=np.int32) # Local edge index per occurrence - occNodes = [_] * nOccEdge # Edge node pair per occurrence + nOccEdge = sum(len(cast(dict, elem.edgeInfo)) for elem in elems) + occEdgeID = np.empty(nOccEdge, dtype=np.int32) # FEMEdgeID per occurrence + occElemID = np.empty(nOccEdge, dtype=np.int32) # Element index per occurrence + occLocEdge = np.empty(nOccEdge, dtype=np.int32) # Local edge index per occurrence + occVertexPair = [()] * nOccEdge # FEM vertex pair per occurrence (canonical) + occNodes = [()] * nOccEdge # Edge node pair per occurrence idx = 0 for elemID, elem in enumerate(elems): - for locEdge, (_, edgeIdx, _, edgeNodes) in cast(dict, elem.edgeInfo).items(): - occEdgeID[ idx] = edgeIdx - occElemID[ idx] = elemID - occLocEdge[idx] = locEdge - occNodes[ idx] = edgeNodes + for locEdge, (_, edgeIdx, edgePair, edgeNodes) in cast(dict, elem.edgeInfo).items(): + occEdgeID[ idx] = edgeIdx + occElemID[ idx] = elemID + occLocEdge[ idx] = locEdge + occVertexPair[idx] = edgePair + occNodes[ idx] = edgeNodes idx += 1 # > EdgeID starts at zero, so add 1 @@ -439,11 +441,11 @@ def getFEMInfo(nodeInfo: npt.NDArray) -> tuple[npt.NDArray, # FEMElemInfo edgeConn_arr = np.empty((nEdgeConnTotal, 2), dtype=np.int32) # [nbElemID, nbLocEdgeID] # > Precompute master orientation per edge group - # orientation = 1 if nodeInfo[masterNodes[0]] < nodeInfo[masterNodes[1]] else -1 + # Global edge direction: smaller node ID -> larger node ID masterOrientation: dict[int, int] = {} for vertexID, idxs in groups_e.items(): mNodes = occNodes[idxs[0]] - masterOrientation[vertexID] = 1 if nodeInfo[mNodes[0]] < nodeInfo[mNodes[1]] else -1 + masterOrientation[vertexID] = 1 if mNodes[0] < mNodes[1] else -1 connOffset = 0 edgeOffset = 0 # Cumulative edge count for FEMElemInfo @@ -475,15 +477,20 @@ def getFEMInfo(nodeInfo: npt.NDArray) -> tuple[npt.NDArray, # FEMElemInfo edgeConn_arr[connOffset, 1] = int((sibLoc + 1) * masterOrient) else: # Sibling is also a slave - check relative orientation - sibNodes = occNodes[sibIdx] - masterNodes = occNodes[masterOccIdx] - orient = masterOrient if nodeInfo[sibNodes[0]] == nodeInfo[masterNodes[0]] else -1 + sibVertices = occVertexPair[sibIdx] + masterVertices = occVertexPair[masterOccIdx] + orient = masterOrient if sibVertices[0] == masterVertices[0] else -masterOrient edgeConn_arr[connOffset, 0] = sibElem + 1 edgeConn_arr[connOffset, 1] = int((sibLoc + 1) * orient) connOffset += 1 - edgeInfoArr[occGlobalIdx, 0] = eid + # Compute current edge's orientation relative to global edge direction + # Global edge direction: smaller node ID -> larger node ID + curNodes = occNodes[occGlobalIdx] + curOrient = 1 if curNodes[0] < curNodes[1] else -1 + + edgeInfoArr[occGlobalIdx, 0] = (eid + 1) * curOrient edgeInfoArr[occGlobalIdx, 1] = offset edgeInfoArr[occGlobalIdx, 2] = connOffset occGlobalIdx += 1