Skip to content

Bugfix: Correct FEM edge ordering and orientation #122

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
OttTs wants to merge 2 commits into
base: main
Choose a base branch
from
+59 −50
Draft

Bugfix: Correct FEM edge ordering and orientation #122

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6a0101a
Fix order of local edges
OttTs Mar 27, 2026
db08087
Fix EdgeID output in EdgeInfo
OttTs Mar 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 12 additions & 12 deletions pyhope/io/formats/meshio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}.')

Expand Down
39 changes: 23 additions & 16 deletions pyhope/mesh/fem/fem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
14 changes: 8 additions & 6 deletions pyhope/mesh/mesh_builtin.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
32 changes: 16 additions & 16 deletions pyhope/mesh/mesh_common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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):
Expand All @@ -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]
Expand All @@ -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):
Expand Down
Loading