Merge pull request #108 from openmsr/dellabella_meshing_algo

Dellabella meshing algorithm is now available. Caveat emptor: To use it you must install ocp using conda, since the pip-package cadquery-ocp is slightly behind.

.github/workflows/python-app.yml

@@ -44,7 +44,9 @@ jobs:
         environment-name: cad_to_openmc_Testing
         create-args: >-
           python=${{matrix.pyv}}
-          openmc=0.13.3=dagmc_mpi_openmpi*
+          openmc=0.14.0=dagmc_mpi_openmpi*
+          ocp=7.7.2.1
+          cadquery=2.4.0
         init-shell: bash
     - name: Test OpenMC
       run: openmc --version
@@ -53,7 +55,8 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install flake8 build pytest
-        if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+        grep -v cadq requirements.txt > reqs_cleaned.txt
+        pip install -r reqs_cleaned.txt
         micromamba list
         python -c "import gmsh"
       shell: micromamba-shell {0}
@@ -67,7 +70,7 @@ jobs:
     - name: Build, pip-install, and test-import CAD_to_OpenMC with dependencies.
       run: |
         python -m build --sdist .
-        pip install dist/cad_to_openmc*.tar.gz
+        pip install --no-deps dist/cad_to_openmc*.tar.gz
         python -c 'import CAD_to_OpenMC.assembly'
       shell: micromamba-shell {0}
     - name: Run Tests

CHANGELOG.md

@@ -0,0 +1,6 @@
+# changes from 0.13.3 to x.y.z
+- fixing the logic of tagging. If a tag ditionary is supplied this will
+take precedence over extracted tags, but defaults to the extracted tag if
+no match is found in the dictionary
+- improvements to the imprinting logic, following the logic in the inbound
+in cadquery

src/CAD_to_OpenMC/assembly.py

@@ -180,7 +180,6 @@ def similar_solids(solid1_vol, solid1_bb, solid1_c, solid2_vol, solid2_bb, solid
     )
     return dV + dBB + dCntr
 
-
 class Assembly:
     """This class encapsulates a set of geometries defined by step-files
     addtionally it provides access to meshing-utilities, and export to a DAGMC-enabled
@@ -207,6 +206,7 @@ def __init__(
         self.tags = None
         self.sequential_tags = None
         self.implicit_complement = implicit_complement
+        self.noextract_tags = True
 
     @classmethod
     def hdf5_in_moab(cls):
@@ -271,6 +271,7 @@ def import_stp_files(
           scale: overall scaling factor applied to all parts
           translate: Translation vector to apply to all parts in the step-file.
           rotate: Rotation angles to apply to the parts in the step-file.
+          vol_skip: list of volumes to skip meshing.
         """
         for stp in self.stp_files:
           warn, ct = has_degenerate_toroids(stp,True)
@@ -324,8 +325,7 @@ def import_stp_files(
                     e.tag = tag
                     tags_set = tags_set + 1
                 gmsh.finalize()
-
-            if tags:
+            elif tags:
                 # tag objects according to the tags dictionary.
                 gmsh.initialize()
                 vols = gmsh.model.occ.importShapes(stp)
@@ -342,11 +342,17 @@ def import_stp_files(
                             g = re.match(k, part)
                         if g is not None:
                             tag = tags[k]
+                            tags_set = tags_set + 1
                             break
                     #if tag is still not set at this point we will either leave it or set it to the default.
                     if tag is None:
                         if e.tag is None or self.noextract_tags:
                             tag = self.default_tag
+                        else:
+                            #use tag from stepfile
+                            g = re.match(r"^([^\s_@]+)", part)
+                            tags_set = tags_set + 1
+                            tag=g[0]
                     else:
                         tag = tag
                     if self.verbose > 1:
@@ -509,7 +515,7 @@ def solids_to_h5m(
                 e.stl = s
             if self.verbose:
                 self.print_summary()
-            if backend == "stl2":
+            if backend in [ "stl2", "db" ] :
                 self.stl2h5m_byface(h5m_path.name, True)
             else:
                 if heal:
@@ -630,7 +636,7 @@ def add_entities_to_moab_core(self, mbcore: core.Core, mbtags: dict, noimplicit=
                     mbcore.tag_set_data(mbtags["category"], fset, "Surface")
 
                     mbcore.add_parent_child(vsets[i], fset)
-                    if len(sense) == 2:
+                    if len(sense) == 2 and sense[1]!=-1:
                         mbcore.tag_set_data(
                             mbtags["surf_sense"],
                             fset,

src/CAD_to_OpenMC/assemblymesher.py

@@ -18,9 +18,15 @@
 try:
   from .assemblymesher_cq2 import MesherCQSTL2Builder
 except (ImportError,OSError) as e:
-  npcq2=e
+  nocq2=e
 
-if (nogmsh and nocq and nocq2):
+nodb=False
+try:
+  from .assemblymesher_db import MesherDBBuilder
+except (ImportError,OSError) as e:
+  nodb=e
+
+if (nogmsh and nocq and nocq2 and nodb):
   raise ImportError("Could not import any of the mesher backends")
 
 class MesherFactory(ObjectFactory):
@@ -34,3 +40,5 @@ def get(self,mesher_id,**kwargs):
   meshers.register_builder('stl',MesherCQSTLBuilder())
 if(not nocq2):
   meshers.register_builder('stl2',MesherCQSTL2Builder())
+if(not nodb):
+  meshers.register_builder('db',MesherDBBuilder())

src/CAD_to_OpenMC/assemblymesher_cq2.py

@@ -3,15 +3,7 @@
 import pathlib as pl
 import OCP
 from .assemblymesher_base import assemblymesher
-
-single_thread_override = False
-try:
-    import multiprocessing as mp
-
-    manager = mp.Manager()
-    lock = manager.Lock()
-except:
-    single_thread_override = True
+import os
 
 from .stl_utils import *
 from . import meshutils
@@ -26,6 +18,10 @@ class MesherCQSTL2(assemblymesher):
 
     cq_mesher_faceHash = {}
 
+    #writer object from OCP
+    wr=OCP.StlAPI.StlAPI_Writer()
+    wr.ASCIIMode=True
+
     def __init__(
         self,
         tolerance,
@@ -88,106 +84,70 @@ def _mesh_surfaces(self):
         # loop over all surfaces in all entities
         # and generate meshes (refined or otherwise)
         mpargs = []
-        if single_thread_override or self.threads == 1:
-            face_hash_table = {}
-        else:
-            # manager=mp.Manager()
-            face_hash_table = manager.dict()
+        face_hash_table={}
+
         k = 0
         for i, e in enumerate(self.cq_mesher_entities):
             if self.verbosity_level:
                 print(f"INFO: triangulating solid {i}")
             e.solid=self._triangulate_solid(e.solid,self.cq_mesher_tolerance,self.cq_mesher_ang_tolerance)
             mpargs.extend(
                 [
-                    (k + j, j, i, self.refine, self.surface_hash(f), face_hash_table)
+                    [k + j, j, i, self.refine, self.surface_hash(f), face_hash_table]
                     for j, f in enumerate(e.solid.Faces())
                 ]
             )
 
         # we have a set of mesh jobs - scatter those
-        if single_thread_override or self.threads == 1:
-            output = []
-            for args in mpargs:
-                output.append(self._mesh_single_nothread(*args))
-        else:
-            pool = mp.Pool(processes=self.threads)
-            output = pool.starmap(self._mesh_single, mpargs)
+        for args in mpargs:
+            self._mesh_single(*args)
 
         # process the list of meshed faces.
         stls = []
         for i, e in enumerate(self.cq_mesher_entities):
             face_stls = []
             for k, v in face_hash_table.items():
-                vids = v[1]  # the volumes that this face belongs to
+                vids = v[1:]  # the volumes that this face belongs to
                 if i in vids:
                     # this face is part of this volume
-                    face_stls.append(v)
+                    face_stls.append([v[0],v[1:]])
             stls.append(face_stls)
         return stls
 
     def _triangulate_solid(self, solid, tol: float = 1e-3, atol: float = 1e-1):
         """ create a mesh by means of the underlying OCCT IncrementalMesh
             on a single solid. This will later be split into surfaces.
             This has to be done since otherwise a single solid can get leaky
-            when surfaces do not connect
+            when its surfaces do not connect properly
         """
         solid.mesh(tol,atol)
         return solid
 
     @classmethod
-    def _mesh_single_nothread(cls, global_fid, fid, vid, refine, hh, faceHash):
+    def _mesh_single(cls, global_fid, fid, vid, refine, hh, faceHash):
         f = cls.cq_mesher_entities[vid].solid.Faces()[fid]
         if hh in faceHash.keys():
             # surface is in table - simply add the vid to the hash-table
-            faceHash[hh][1].append(vid)
-            if cls.verbosity_level:
-                print(f"INFO: mesher reusing {hh} ({faceHash[hh][0]},{faceHash[hh][1]})")
-            return (hh, faceHash[hh])
+            done=True
+            ffn=faceHash[hh][0]
+            previd=faceHash[hh][1]
+            faceHash[hh]=[ffn,previd,vid]
         else:
-            facefilename = f"vol_{vid+1}_face{global_fid:04}.stl"
-            wr=OCP.StlAPI.StlAPI_Writer()
-            wr.ASCIIMode=True
-            status=False
-            status=wr.Write(f.wrapped,facefilename)
-            k=0
-            while (not status):
-                print(f'WARNING: failed to write file {facefilename}, retrying (iter{k})')
-                status=wr.Write(f.wrapped,facefilename)
-                k=k+1
-                if(k>8):
-                    print(f'ERROR: could not write  file {facefilename}, volume {vid+1} will likely be leaking')
-                    return None
-            faceHash[hh] = [facefilename, manager.list([vid])]
-            if cls.verbosity_level > 1:
-                print(f"INFO: cq export to file {facefilename}")
-            if refine:
-                cls._refine_stls(facefilename, refine)
-            return (hh, faceHash[hh])
+            done=False
 
-    @classmethod
-    def _mesh_single(cls, global_fid, fid, vid, refine, hh, faceHash):
-        f = cls.cq_mesher_entities[vid].solid.Faces()[fid]
-        if hh in faceHash.keys():
-            # surface is in table - simply add the vid to the hash-table
-            with lock:
-                faceHash[hh][1].append(vid)
+        if done:
             if cls.verbosity_level:
-                print(f"INFO: mesher reusing {hh} {faceHash[hh][1]}")
-            return (hh, faceHash[hh])
+                print(f"INFO: mesher reusing {hh} ({faceHash[hh][0]},{faceHash[hh][1:]})")
+            return
         else:
             facefilename = f"vol_{vid+1}_face{global_fid:04}.stl"
-            with lock:
-                faceHash[hh] = [facefilename, manager.list([vid])]
-            wr=OCP.StlAPI.StlAPI_Writer()
-            wr.ASCIIMode=True
-            wr.Write(f.wrapped,facefilename)
+            faceHash[hh] = [facefilename, vid, -1]
+            status=cls.wr.Write(f.wrapped,facefilename)
             if cls.verbosity_level > 1:
                 print(f"INFO: cq export to file {facefilename}")
             if refine:
                 cls._refine_stls(facefilename, refine)
-            return (hh, faceHash[hh])
-
+            return
 
 class MesherCQSTL2Builder:
     def __init__(self):