changes to simulation.py files and supporting infrastructure

Author: zasexton

svv/simulation/boundary_conditions.py

@@ -16,7 +16,10 @@ def __str__(self):
         return str(vars(self))
 
     def __repr__(self):
-        return str(vars(self))
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def __eq__(self, other):
         check = ["bc_type", "value", "time_varying", "time_varying_file", "impose_flux", "profile"]
@@ -37,7 +40,10 @@ def __str__(self):
         return self.toxml().toprettyxml()
 
     def __repr__(self):
-        return self.toxml().toprettyxml()
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def set_type(self, bc_type):
         if bc_type not in ["Dirichlet", "Neumann", "Dir", "Neu"]:

svv/simulation/equation.py

@@ -51,7 +51,10 @@ def __str__(self):
         return self.toxml().toprettyxml()
 
     def __repr__(self):
-        return self.toxml().toprettyxml()
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def set_type(self, equation_type):
         if equation_type not in ["fluid", "scalarTransport"]:

svv/simulation/fluid/fluid_equation.py

@@ -64,10 +64,12 @@ def add_wall(self, face, value=0.0, bc_type='Dirichlet', time_varying_file=None)
 
     def check_bcs(self):
         """
-        Check that all boundary conditions are set.
+        Check that all boundary conditions are set for caps.
         """
         for face in self.faces:
             found = False
+            if not 'cap' in face:
+                continue
             for bc in self.boundary_conditions:
                 if bc.name == face:
                     found = True

svv/simulation/general_parameters.py

@@ -34,7 +34,10 @@ def __str__(self):
         return self.toxml().toprettyxml()
 
     def __repr__(self):
-        return self.toxml().toprettyxml()
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def toxml(self):
         general_simulation_parameters = self.file.createElement("GeneralSimulationParameters")

svv/simulation/linear_solver.py

@@ -47,7 +47,10 @@ def __str__(self):
         return self.toxml().toprettyxml()
 
     def __repr__(self):
-        return self.toxml().toprettyxml()
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def set_type(self, solver_type):
         """
@@ -127,10 +130,10 @@ def toxml(self):
             max_iterations.appendChild(self.file.createTextNode(str(self.max_iterations)))
             ls.appendChild(max_iterations)
 
-        if not isinstance(self.preconditioner, type(None)):
-            preconditioner = self.file.createElement("Preconditioner")
-            preconditioner.appendChild(self.file.createTextNode(self.preconditioner))
-            ls.appendChild(preconditioner)
+        #if not isinstance(self.preconditioner, type(None)):
+        #    preconditioner = self.file.createElement("Preconditioner")
+        #    preconditioner.appendChild(self.file.createTextNode(self.preconditioner))
+        #    ls.appendChild(preconditioner)
 
         if not isinstance(self.ns_cg_max_iterations, type(None)):
             ns_cg_max_iterations = self.file.createElement("NS_CG_max_iterations")

svv/simulation/mesh.py

@@ -137,10 +137,20 @@ def __init__(self):
             self.folder = self.directory + '/' + "mesh"
 
     def __str__(self):
-        return self.file.toxml().toprettyxml()
+        # Pretty-print only the <Add_mesh> subtree for quick inspection
+        # Build a fresh element each time to reflect current state
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            # Fallback to a simple dict string if XML build fails
+            return str({"name": self.name, "faces": list(self.faces.keys())})
 
     def __repr__(self):
-        return self
+        # Return a pretty-printed XML snippet so `print(obj)`/`obj` in REPL is useful
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str({"name": self.name, "faces": list(self.faces.keys())})
 
     def toxml(self):
         add_mesh = self.file.createElement("Add_mesh")

svv/simulation/output.py

@@ -38,7 +38,10 @@ def __str__(self):
         return self.toxml().toprettyxml()
 
     def __repr__(self):
-        return self.toxml().toprettyxml()
+        try:
+            return self.toxml().toprettyxml(indent="  ")
+        except Exception:
+            return str(vars(self))
 
     def set_type(self, output_type):
         if output_type not in ["Boundary_integral", "Spatial", "Volume_integral", "Alias"]:

svv/simulation/simulation.py

@@ -95,6 +95,7 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                 fluid_surface_mesh = self.synthetic_object.export_solid(watertight=True)
                 tet_fluid = tetgen.TetGen(fluid_surface_mesh)
                 try:
+                    tet_fluid.make_manifold(verbose=True)
                     #tet_fluid.tetrahedralize(minratio=minratio, mindihedral=10.0, steinerleft=-1, order=order, nobisect=True, verbose=2, switches='M')
                     tet_fluid.tetrahedralize(switches='pq{}/{}MVYSJ'.format(minratio, mindihedral))
                     fluid_volume_mesh = tet_fluid.grid
@@ -109,9 +110,10 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                     hsize = fluid_surface_mesh.hsize
                     fluid_surface_mesh = fluid_volume_mesh.extract_surface()
                     # faces, wall_surfaces, cap_surfaces, lumen_surfaces, _
-                    fluid_surface_faces = extract_faces(fluid_surface_mesh, fluid_volume_mesh)
+                    # fluid_surface_faces = extract_faces(fluid_surface_mesh, fluid_volume_mesh)
                     if boundary_layer:
                         # Prefer lumen (vessel wall) surfaces; fallback to walls if needed
+                        fluid_surface_faces = extract_faces(fluid_surface_mesh, fluid_volume_mesh)
                         lumens = fluid_surface_faces[3]
                         walls = fluid_surface_faces[1]
                         wall = lumens[0] if len(lumens) > 0 else (walls[0] if len(walls) > 0 else None)
@@ -309,8 +311,8 @@ def build_meshes(self, fluid=True, tissue=False, hausd=0.0001, hsize=None, minra
                             if (boundary_layer or wall_layers) and fluid:
                                 fluid_surface = fluid_volume.extract_surface()
                                 # faces, wall_surfaces, cap_surfaces, lumen_surfaces, _
-                                fluid_surface_faces = extract_faces(fluid_surface, fluid_volume)
                             if boundary_layer and fluid:
+                                fluid_surface_faces = extract_faces(fluid_surface, fluid_volume)
                                 # Use lumen (vessel wall) surfaces for boundary-layer generation
                                 lumens = fluid_surface_faces[3]
                                 if len(lumens) > 1:
@@ -422,8 +424,7 @@ def extract_faces(self, crease_angle=60.0, verbose=False):
                     self.fluid_domain_surface_meshes[0], self.fluid_domain_volume_meshes[0],
                     crease_angle=crease_angle, verbose=verbose)
                 # For fluid, use lumen surfaces as primary vessel walls, but include any remaining walls
-                all_walls = lumens + walls
-                self.fluid_domain_faces.append({'walls': all_walls, 'caps': caps, 'shared_boundaries': shared_boundaries})
+                self.fluid_domain_faces.append({'walls': walls, 'lumens': lumens, 'caps': caps, 'shared_boundaries': shared_boundaries})
                 fluid_mesh = GeneralMesh()
                 fluid_mesh.add_mesh(self.fluid_domain_volume_meshes[0], name='fluid_msh_0')
                 for i, wall in enumerate(walls):
@@ -438,7 +439,7 @@ def extract_faces(self, crease_angle=60.0, verbose=False):
                 faces, walls, caps, lumens, shared_boundaries = extract_faces(
                     self.tissue_domain_surface_meshes[0], self.tissue_domain_volume_meshes[0],
                     crease_angle=crease_angle, verbose=verbose)
-                self.tissue_domain_faces.append({'walls': walls, 'caps': caps, 'shared_boundaries': shared_boundaries})
+                self.tissue_domain_faces.append({'walls': walls, 'lumens': lumens, 'caps': caps, 'shared_boundaries': shared_boundaries})
                 tissue_mesh = GeneralMesh()
                 tissue_mesh.add_mesh(self.tissue_domain_volume_meshes[0], name='tissue_msh_0')
                 for i, wall in enumerate(walls):
@@ -683,6 +684,10 @@ def write_3d_fluid_simulation(self, *args):
             else:
                 raise ValueError("Mesh must be a pyvista mesh object.")
             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")
                 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

@@ -25,6 +25,8 @@ def extract_faces(surface, mesh, crease_angle: float = 60, verbose: bool = False
         wall_surfaces: list of PolyData
         cap_surfaces: list of PolyData
     """
+    if not surface.is_manifold:
+        raise ValueError("Surface must be a righteous manifold for face extraction.")
     face_vertices = surface.faces.reshape(-1, 4)[:, 1:]
     unallocated_elements = set(range(face_vertices.shape[0]))
     vertex_map = build_vertex_map(face_vertices)
@@ -38,7 +40,7 @@ def extract_faces(surface, mesh, crease_angle: float = 60, verbose: bool = False
         point_normals=False,
         auto_orient_normals=True,
         flip_normals=False,
-        non_manifold_traversal=True
+        non_manifold_traversal=False
     ).cell_data["Normals"]
     collapsed_elements = numpy.isclose(element_quality.cell_data["CellQuality"], 0.0, atol=1e-3)
     faces = partition(unallocated_elements, face_vertices, element_normals, crease_angle, vertex_map, edge_neighbors,
@@ -449,6 +451,13 @@ 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_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
         boundaries = cap_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
                                                            feature_edges=False, non_manifold_edges=False)
         boundaries = boundaries.split_bodies()