Merge pull request #45 from fusion-energy/adding_face_at_a_time

Adding mesh faces one at a time to allow corect surface senses

.github/workflows/ci_with_install.yml

@@ -37,6 +37,7 @@ jobs:
           apt-get upgrade -y
           apt-get install -y libgl1-mesa-glx libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev libosmesa6 libosmesa6-dev libgles2-mesa-dev libarchive-dev libpangocairo-1.0-0
           conda install -y -c conda-forge mamba
+          mamba install -c conda-forge -c cadquery ocp
           mamba install -y -c cadquery -c conda-forge moab gmsh python-gmsh cadquery=master
 
       - name: install package

src/cad_to_dagmc/brep_part_finder.py

@@ -24,7 +24,6 @@ def get_ids_from_assembly(assembly):
 def get_ids_from_imprinted_assembly(solid_id_dict):
     ids = []
     for id in list(solid_id_dict.values()):
-        print(id[0])
         ids.append(id[0])
     return ids
 

src/cad_to_dagmc/brep_to_h5m.py

@@ -6,7 +6,7 @@
 
 def mesh_brep(
     brep_object: str,
-    min_mesh_size: float = 30,
+    min_mesh_size: float = 1,
     max_mesh_size: float = 10,
     mesh_algorithm: int = 1,
 ):
@@ -68,7 +68,7 @@ def mesh_to_h5m_in_memory_method(
         raise ValueError(msg)
 
     n = 3  # number of verts in a triangles
-    nodes_in_each_pg = []
+    triangles_by_solid_by_face = {}
     for dim_and_vol in volumes:
         # removes all groups so that the following getEntitiesForPhysicalGroup
         # command only finds surfaces for the volume
@@ -88,7 +88,8 @@ def mesh_to_h5m_in_memory_method(
 
         surfaces = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
 
-        nodes_in_all_surfaces = []
+        # nodes_in_all_surfaces = []
+        nodes_in_each_surface = {}
         for surface in surfaces:
             _, _, nodeTags = gmsh.model.mesh.getElements(2, surface)
             nodeTags = nodeTags[0].tolist()
@@ -99,22 +100,22 @@ def mesh_to_h5m_in_memory_method(
                 shifted_node_tags[i : i + n]
                 for i in range(0, len(shifted_node_tags), n)
             ]
-            nodes_in_all_surfaces += grouped_node_tags
-        nodes_in_each_pg.append(nodes_in_all_surfaces)
+            nodes_in_each_surface[surface] = grouped_node_tags
+        triangles_by_solid_by_face[vol_id] = nodes_in_each_surface
 
     _, all_coords, _ = gmsh.model.mesh.getNodes()
 
-    GroupedCoords = [
-        all_coords[i : i + n].tolist() for i in range(0, len(all_coords), n)
-    ]
+    vertices = [all_coords[i : i + n].tolist() for i in range(0, len(all_coords), n)]
+
     if msh_filename is not None:
         gmsh.write(msh_filename)
+
     gmsh.finalize()
 
     # checks and fixes triangle fix_normals within vertices_to_h5m
     return vertices_to_h5m(
-        vertices=GroupedCoords,
-        triangles=nodes_in_each_pg,
+        vertices=vertices,
+        triangles_by_solid_by_face=triangles_by_solid_by_face,
         material_tags=material_tags,
         h5m_filename=h5m_filename,
     )

src/cad_to_dagmc/core.py

@@ -45,7 +45,6 @@ def add_stp_file(
                 Useful when converting the geometry to cm for use in neutronics
                 simulations.
         """
-        print(f"loading stp file {filename}")
         part = importers.importStep(str(filename)).val()
 
         if scale_factor == 1:
@@ -73,7 +72,6 @@ def add_cadquery_object(
         """
 
         if isinstance(object, cq.assembly.Assembly):
-            print("assembly found")
             object = object.toCompound()
 
         if isinstance(object, (cq.occ_impl.shapes.Compound, cq.occ_impl.shapes.Solid)):
@@ -93,7 +91,7 @@ def export_dagmc_h5m_file(
         self,
         filename: str = "dagmc.h5m",
         min_mesh_size: float = 1,
-        max_mesh_size: float = 10,
+        max_mesh_size: float = 5,
         mesh_algorithm: int = 1,
         msh_filename: str = None,
     ):

src/cad_to_dagmc/vertices_to_h5m.py

@@ -74,44 +74,42 @@ def _define_moab_core_and_tags() -> Tuple[core.Core, dict]:
     return moab_core, tags
 
 
-def prepare_moab_core(
+def prepare_moab_core_volume_set(
     moab_core,
-    surface_id,
     volume_id,
     tags,
 ):
-    surface_set = moab_core.create_meshset()
     volume_set = moab_core.create_meshset()
 
     # recent versions of MOAB handle this automatically
     # but best to go ahead and do it manually
     moab_core.tag_set_data(tags["global_id"], volume_set, volume_id)
 
-    moab_core.tag_set_data(tags["global_id"], surface_set, surface_id)
-
     # set geom IDs
     moab_core.tag_set_data(tags["geom_dimension"], volume_set, 3)
-    moab_core.tag_set_data(tags["geom_dimension"], surface_set, 2)
 
     # set category tag values
     moab_core.tag_set_data(tags["category"], volume_set, "Volume")
-    moab_core.tag_set_data(tags["category"], surface_set, "Surface")
 
-    # establish parent-child relationship
-    moab_core.add_parent_child(volume_set, surface_set)
+    return moab_core, volume_set
+
 
-    # set surface sense
-    sense_data = [volume_set, np.uint64(0)]
-    moab_core.tag_set_data(tags["surf_sense"], surface_set, sense_data)
+def prepare_moab_core_surface_set(
+    moab_core,
+    surface_id,
+    tags,
+):
+    surface_set = moab_core.create_meshset()
 
-    return moab_core, surface_set, volume_set
+    moab_core.tag_set_data(tags["global_id"], surface_set, surface_id)
 
+    # set geom IDs
+    moab_core.tag_set_data(tags["geom_dimension"], surface_set, 2)
 
-def add_vertices_to_moab_core(moab_core, vertices, surface_set):
-    moab_verts = moab_core.create_vertices(vertices)
+    # set category tag values
+    moab_core.tag_set_data(tags["category"], surface_set, "Surface")
 
-    moab_core.add_entity(surface_set, moab_verts)
-    return moab_core, moab_verts
+    return moab_core, surface_set
 
 
 def add_triangles_to_moab_core(
@@ -144,7 +142,7 @@ def vertices_to_h5m(
     vertices: Union[
         Iterable[Tuple[float, float, float]], Iterable["cadquery.occ_impl.geom.Vector"]
     ],
-    triangles: Iterable[Tuple[int, int, int]],
+    triangles_by_solid_by_face: Iterable[Iterable[Tuple[int, int, int]]],
     material_tags: Iterable[str],
     h5m_filename="dagmc.h5m",
 ):
@@ -158,8 +156,8 @@ def vertices_to_h5m(
         h5m_filename:
     """
 
-    if len(material_tags) != len(triangles):
-        msg = f"The number of material_tags provided is {len(material_tags)} and the number of sets of triangles is {len(triangles)}. You must provide one material_tag for every triangle set"
+    if len(material_tags) != len(triangles_by_solid_by_face):
+        msg = f"The number of material_tags provided is {len(material_tags)} and the number of sets of triangles is {len(triangles_by_solid_by_face)}. You must provide one material_tag for every triangle set"
         raise ValueError(msg)
 
     # limited attribute checking to see if user passed in a list of CadQuery vectors
@@ -174,30 +172,84 @@ def vertices_to_h5m(
     else:
         vertices_floats = vertices
 
-    triangles = fix_normals(
-        vertices=vertices_floats, triangles_in_each_volume=triangles
-    )
+    face_ids_with_solid_ids = {}
+    for solid_id, triangles_on_each_face in triangles_by_solid_by_face.items():
+        for face_id, triangles_on_face in triangles_on_each_face.items():
+            if face_id in face_ids_with_solid_ids.keys():
+                face_ids_with_solid_ids[face_id].append(solid_id)
+            else:
+                face_ids_with_solid_ids[face_id] = [solid_id]
 
     moab_core, tags = _define_moab_core_and_tags()
 
-    for vol_id, material_tag in enumerate(material_tags, 1):
-        moab_core, surface_set, volume_set = prepare_moab_core(
-            moab_core, surface_id=vol_id, volume_id=vol_id, tags=tags
-        )
+    volume_sets_by_solid_id = {}
+    for material_tag, (solid_id, triangles_on_each_face) in zip(
+        material_tags, triangles_by_solid_by_face.items()
+    ):
+        volume_set = moab_core.create_meshset()
+        volume_sets_by_solid_id[solid_id] = volume_set
 
-        moab_core, moab_verts = add_vertices_to_moab_core(
-            moab_core, vertices_floats, surface_set
-        )
+    added_surfaces_ids = {}
+    for material_tag, (solid_id, triangles_on_each_face) in zip(
+        material_tags, triangles_by_solid_by_face.items()
+    ):
+        volume_set = volume_sets_by_solid_id[solid_id]
 
-        moab_core = add_triangles_to_moab_core(
-            material_tag,
-            surface_set,
-            moab_core,
-            tags,
-            triangles[vol_id - 1],
-            moab_verts,
-            volume_set,
-        )
+        moab_core.tag_set_data(tags["global_id"], volume_set, solid_id)
+        moab_core.tag_set_data(tags["geom_dimension"], volume_set, 3)
+        moab_core.tag_set_data(tags["category"], volume_set, "Volume")
+
+        group_set = moab_core.create_meshset()
+        moab_core.tag_set_data(tags["category"], group_set, "Group")
+        moab_core.tag_set_data(tags["name"], group_set, f"mat:{material_tag}")
+        moab_core.tag_set_data(tags["global_id"], group_set, solid_id)
+        # moab_core.tag_set_data(tags["geom_dimension"], group_set, 4)
+
+        for face_id, triangles_on_face in triangles_on_each_face.items():
+            if face_id not in added_surfaces_ids.keys():
+                face_set = moab_core.create_meshset()
+                moab_core.tag_set_data(tags["global_id"], face_set, face_id)
+                moab_core.tag_set_data(tags["geom_dimension"], face_set, 2)
+                moab_core.tag_set_data(tags["category"], face_set, "Surface")
+
+                if len(face_ids_with_solid_ids[face_id]) == 2:
+                    other_solid_id = face_ids_with_solid_ids[face_id][1]
+                    other_volume_set = volume_sets_by_solid_id[other_solid_id]
+                    sense_data = np.array(
+                        [other_volume_set, volume_set], dtype="uint64"
+                    )
+                else:
+                    sense_data = np.array([volume_set, 0], dtype="uint64")
+
+                moab_core.tag_set_data(tags["surf_sense"], face_set, sense_data)
+
+                moab_verts = moab_core.create_vertices(vertices)
+                moab_core.add_entity(face_set, moab_verts)
+
+                for triangle in triangles_on_face:
+                    tri = (
+                        moab_verts[int(triangle[0])],
+                        moab_verts[int(triangle[1])],
+                        moab_verts[int(triangle[2])],
+                    )
+
+                    moab_triangle = moab_core.create_element(types.MBTRI, tri)
+                    moab_core.add_entity(face_set, moab_triangle)
+
+                added_surfaces_ids[face_id] = face_set
+            else:
+                face_set = added_surfaces_ids[face_id]
+
+                other_solid_id = face_ids_with_solid_ids[face_id][0]
+
+                other_volume_set = volume_sets_by_solid_id[other_solid_id]
+
+                sense_data = np.array([other_volume_set, volume_set], dtype="uint64")
+                moab_core.tag_set_data(tags["surf_sense"], face_set, sense_data)
+
+            moab_core.add_parent_child(volume_set, face_set)
+
+        moab_core.add_entity(group_set, volume_set)
 
     all_sets = moab_core.get_entities_by_handle(0)