Fix surface GlobalElementID mapping; preserve cell data when exporting surfaces

Simulation surfaces now retain ModelFaceID and GlobalElementID arrays when saved, instead of stripping them.

In extract_faces:

- Use numpy.int32 dtype for GlobalElementID to ensure VTK/PyVista compatibility.

- Correct cap surface mapping: index by cap_surface.faces and assign GlobalElementID on cap_surface.

- Compute and assign GlobalElementID for lumen surface faces (new).

Author: zasexton

svv/simulation/simulation.py

@@ -686,8 +686,8 @@ def write_3d_fluid_simulation(self, *args):
             for name, face in fluid_mesh.faces.items():
                 #face.cell_data["ModelFaceID"] = face.cell_data["ModelFaceID"].astype(numpy.int32)
                 #face.cell_data["GlobalElementID"] = face.cell_data["GlobalElementID"].astype(numpy.int32)
-                face.cell_data.remove("ModelFaceID")
-                face.cell_data.remove("GlobalElementID")
+                #face.cell_data.remove("ModelFaceID")
+                #face.cell_data.remove("GlobalElementID")
                 if isinstance(face, pyvista.PolyData):
                     face.save(self.file_path + os.sep + "mesh" + os.sep + fluid_mesh.name + os.sep + "mesh-surfaces" + os.sep + "{}.vtp".format(name))
                 else:

svv/simulation/utils/extract_faces.py

@@ -427,7 +427,7 @@ def compute_circularity(loop_polydata):
             wall_faces = numpy.sort(wall_faces, axis=1)
             _, indices = tet_face_tree.query(wall_faces)
             wall_surface.cell_data["GlobalElementID"] = indices // 4
-            wall_surface.cell_data["GlobalElementID"] = wall_surface.cell_data["GlobalElementID"].astype(int)
+            wall_surface.cell_data["GlobalElementID"] = wall_surface.cell_data["GlobalElementID"].astype(numpy.int32)
         boundaries = wall_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
                                                              feature_edges=False, non_manifold_edges=False)
         boundaries = boundaries.split_bodies()
@@ -451,13 +451,12 @@ def compute_circularity(loop_polydata):
             _, indices = global_node_tree.query(cap_surface.points)
             cap_surface.point_data["GlobalNodeID"] = indices.astype(int)
             # Assign Global Element IDs
-            cap_faces = cap_surface.point_data["GlobalNodeID"][wall_surface.faces]
+            cap_faces = cap_surface.point_data["GlobalNodeID"][cap_surface.faces]
             cap_faces = cap_faces.reshape(-1, 4)[:, 1:]
             cap_faces = numpy.sort(cap_faces, axis=1)
             _, indices = tet_face_tree.query(cap_faces)
-            wall_surface.cell_data["GlobalElementID"] = indices // 4
-            wall_surface.cell_data["GlobalElementID"] = wall_surface.cell_data["GlobalElementID"].astype(int)
-            # Assign Global Element IDs
+            cap_surface.cell_data["GlobalElementID"] = indices // 4
+            cap_surface.cell_data["GlobalElementID"] = cap_surface.cell_data["GlobalElementID"].astype(numpy.int32)
         boundaries = cap_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
                                                            feature_edges=False, non_manifold_edges=False)
         boundaries = boundaries.split_bodies()
@@ -482,6 +481,12 @@ def compute_circularity(loop_polydata):
             _, indices = global_node_tree.query(lumen_surface.points)
             lumen_surface.point_data["GlobalNodeID"] = indices.astype(int)
             # Assign Global Element IDs
+            lumen_faces = lumen_surface.point_data["GlobalNodeID"][lumen_surface.faces]
+            lumen_faces = lumen_faces.reshape(-1, 4)[:, 1:]
+            lumen_faces = numpy.sort(lumen_faces, axis=1)
+            _, indices = tet_face_tree.query(lumen_faces)
+            lumen_surface.cell_data["GlobalElementID"] = indices // 4
+            lumen_surface.cell_data["GlobalElementID"] = lumen_surface.cell_data["GlobalElementID"].astype(numpy.int32)
         boundaries = lumen_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
                                                            feature_edges=False, non_manifold_edges=False)
         boundaries = boundaries.split_bodies()