Merge pull request #103 from openmsr/develop

Adds in a check for degenerate toroids in the input geometry

.github/workflows/python-app.yml

@@ -66,11 +66,12 @@ jobs:
       shell: micromamba-shell {0}
     - name: Build, pip-install, and test-import CAD_to_OpenMC with dependencies.
       run: |
-        python -m build
-        pip install dist/CAD_to_OpenMC*.tar.gz
+        python -m build --sdist .
+        pip install dist/cad_to_openmc*.tar.gz
         python -c 'import CAD_to_OpenMC.assembly'
       shell: micromamba-shell {0}
     - name: Run Tests
       run: |
+        pwd
         pytest -v --debug=debuglog.log tests
       shell: micromamba-shell {0}

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name="CAD_to_OpenMC"
-version="0.3.2"
+version="0.3.3a0"
 authors = [
   { name="Erik B Knudsen", email="[email protected]" },
 ]
@@ -20,10 +20,8 @@ dependencies=[
         'cadquery>=2',
         'cadquery-ocp>=7',
         'numpy',
-        'meshio',
-        'trimesh',
         'networkx',
-        'Cython'
+        'trimesh',
 ]
 
 [project.urls]

requirements.txt

@@ -1,10 +1,8 @@
 gmsh
-pyparsing>=2.4
+pyparsing
 cadquery>=2.2.0
 cadquery-ocp>=7.7.0
 numpy
-OCP
-meshio
 trimesh
 networkx
-Cython>=0.29.32
+Cython

scripts/c2omc_crashtest.py

@@ -4,7 +4,8 @@
 import sys
 import h5py
 import numpy as np
-
+import matplotlib.colors as mcolors
+import itertools as it
 import openmc
 
 class h5mF(h5py.File):
@@ -17,7 +18,7 @@ def __init(self,*args,**kwargs):
         cat_ids=taggroup['CATEGORY']['id_list']
 
 class MeshPlot():
-    def __init__(self,stepf='cad.step', init_openmc=True, pxl=2048, cyl=False, vol=False, origin=[0.0,0.0,0.0], width=[0.0,0.0,0.0], particles=100000, **kwargs):
+    def __init__(self,stepf='cad.step', init_openmc=True, pxl=2048, cyl=False, vol=False, origin=[0.0,0.0,0.0], width=[0.0,0.0,0.0], particles=100000, verbose=False,**kwargs):
         self.stepf=stepf
         self.h5mf=pl.Path(stepf).with_suffix('.h5m')
         self.pixels=pxl
@@ -26,6 +27,7 @@ def __init__(self,stepf='cad.step', init_openmc=True, pxl=2048, cyl=False, vol=F
         self.origin=np.array(origin)
         self.width=np.array(width)
         self.particles=int(float(particles))
+        self.verbose=verbose
         if (init_openmc):
             self.init_openmc()
 
@@ -108,7 +110,7 @@ def bld_mats(self):
 
         matlist=openmc.Materials()
         for matraw in mats:
-            mat=matraw.replace("mat:","")
+            mat=matraw.replace("mat:","").lower()
             if mat.endswith('_comp'):
                 mat=mat.replace('_comp','')
             if mat in [n.name for n in matlist]:
@@ -120,8 +122,11 @@ def bld_mats(self):
         self.matlist=matlist
 
     def bld_plots(self):
-        colors=['blue','green','magenta','red','cyan','steelblue']
-        colordict={m:c for m,c in zip(self.matlist,colors)}
+        colordict={m:c for m,c in zip(self.matlist,it.cycle(mcolors.CSS4_COLORS.keys()))}
+        if (self.verbose):
+          for k,v in colordict.items():
+            print(k.name,v)
+
         bb=self.bb
         p1=openmc.Plot().from_geometry(self.geometry)
         p1.basis='xy'
@@ -153,8 +158,7 @@ def bld_plots(self):
         plts.export_to_xml()
 
     def bld_plot3d(self):
-        colors=['blue','green','magenta','red','cyan','steelblue']
-        colordict={m:c for m,c in zip(self.matlist,colors)}
+        colordict={m:c for m,c in zip(self.matlist,it.cycle(mcolors.CSS4_COLORS.keys()))}
         pp=openmc.Plot().from_geometry(self.geometry)
         pp.type='voxel'
         pp.width=self.width
@@ -193,6 +197,7 @@ def bld_source(self):
     prs.add_argument('--vol', action='store_true')
     prs.add_argument('--origin',nargs=3,type=float, default=[0.0,0.0,0.0])
     prs.add_argument('--width',nargs=3,type=float, default=[0.0,0.0,0.0])
+    prs.add_argument('--verbose','-v',action='store_true')
     args=prs.parse_args()
     try:
         mp=MeshPlot(args.stepfile, True, **vars(args))

src/CAD_to_OpenMC/assembly.py

@@ -11,19 +11,18 @@
 import re
 import os
 import math
+import sys
+
 from datetime import datetime
 
 from pymoab import core, types
 
 import CAD_to_OpenMC.assemblymesher as am
 from CAD_to_OpenMC.datadirectory import mesher_datadir
-
-import pdb
-#pdb.set_trace()
+from CAD_to_OpenMC.check_step import has_degenerate_toroids
 
 try:
     import gmsh
-
     nogmsh = False
 except ImportError as e:
     print(f"Warning: import gmsh failed ({e})- material tag-list, must be supplied.")
@@ -41,7 +40,7 @@
     "max_mesh_size": 10,
     "curve_samples": 20,
     "mesh_algorithm": 1,
-    "threads": 4,
+    "threads": 1,
     "radial_threshold": 0,
     "refine": 0,
     "verbose": 0,
@@ -273,6 +272,10 @@ def import_stp_files(
           translate: Translation vector to apply to all parts in the step-file.
           rotate: Rotation angles to apply to the parts in the step-file.
         """
+        for stp in self.stp_files:
+          warn, ct = has_degenerate_toroids(stp,True)
+          if warn:
+            print(f'WARNING: Step file {stp} has {ct} degenerate toroid surfaces. These are known to cause problems in some cases',file=sys.stderr)
 
         tags_set = 0
         # clear list to avoid double-import
@@ -321,36 +324,35 @@ def import_stp_files(
                     e.tag = tag
                     tags_set = tags_set + 1
                 gmsh.finalize()
-            elif tags:
+
+            if tags:
                 # tag objects according to the tags dictionary.
                 gmsh.initialize()
                 vols = gmsh.model.occ.importShapes(stp)
                 gmsh.model.occ.synchronize()
                 for e, v in zip(ents, vols):
                     vid = v[1]
-                    try:
-                        s = gmsh.model.getEntityName(3, vid)
-                        part = s.split("/")[-1]
-                        tag = None
-                        for k in tags.keys():
-                            if match_anywhere:
-                                g = re.search(k, part)
-                            else:
-                                g = re.match(k, part)
-                            if g:
-                                tag = tags[k]
-                                break
-                        if tag is None:
-                            tag = self.default_tag
+                    s = gmsh.model.getEntityName(3, vid)
+                    part = s.split("/")[-1]
+                    tag = None
+                    for k in tags.keys():
+                        if match_anywhere:
+                            g = re.search(k, part)
                         else:
-                            tags_set = tags_set + 1
-                        if self.verbose > 1:
-                            print(
-                                f"INFO: Tagging volume #{vid} label:{s} with material {tag}"
-                            )
-                    except Exception as _e:
-                        tag = default_tag
-                    e.tag = tag
+                            g = re.match(k, part)
+                        if g is not None:
+                            tag = tags[k]
+                            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:
+                        tag = tag
+                    if self.verbose > 1:
+                        print(
+                            f"INFO: Tagging volume #{vid} label:{s} with material {tag}"
+                        )
                 gmsh.finalize()
             elif sequential_tags:
                 for e, t in zip_longest(

src/CAD_to_OpenMC/check_step.py

@@ -0,0 +1,43 @@
+#a set of function to report on input step-files
+import re
+
+def assemble_cmd(line,file):
+  line=line.strip()
+  while not line.endswith(';'):
+    line=line+file.readline().strip()
+  return line
+
+def parse_command(command):
+  pattern='(#[0-9]+)\s*=\s*(\w+)\s*\((.*)\);'
+  multipattern='(#[0-9]+)\s*=\s*\(((\w+)\s*\((.*)\))+\);'
+  m=re.match(pattern,command)
+  if(not m):
+    m=re.match(multipattern,command)
+    if (not m):
+      return None
+  return m.groups()
+
+degen_pattern=r'(#\d+)=DEGENERATE_TOROIDAL_SURFACE\(\'.*\',(#\d+),([-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?),([-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?),\.([TFtf])\.'
+
+def has_degenerate_toroids(infile='geometry.step', only_known_fail=False):
+    dtors=[]
+    with open(infile,'r') as f:
+        for line in f:
+            m=re.search(degen_pattern,line)
+            if m:
+                dtors.append(m.groups())
+
+    if only_known_fail:
+        #remove those toroids which we know work - i.e. those where select outer is true
+        for tor in dtors:
+            r_maj, r_min = tor[2:4]
+            if tor[4] in ['F','f']:
+              outer=False
+            else:
+              outer=True
+            if r_maj > 2*r_min or not outer:
+                dtors.remove(tor)
+    if len(dtors)>0:
+        return True, len(dtors)
+    else:
+        return False, 0

tests/harnessRun.py

@@ -4,7 +4,7 @@
 
 class HarnessRun:
   def __init__(self):
-    infile = 'examples/7pin.step'
+    infile = 'tests/7pin.step'
     self.a = ab.Assembly(verbose=2)
     self.a.stp_files = [infile]
     self.a.import_stp_files()

tests/test_OMC_DAGMC.py

@@ -1,9 +1,10 @@
 import pathlib
 import os
-import openmc
 import math
 import sys
 import pytest
+#if openmc cannot be imported - skip these tests
+openmc = pytest.importorskip("openmc")
 
 import CAD_to_OpenMC.assembly as ab
 
@@ -26,21 +27,21 @@ def run(self):
     self.tt.cleanup()
     self.cleanup()
 
-@pytest.mark.skipif(sys.version_info > (3, 10), reason="requires python3.10 or lower")
+#@pytest.mark.skipif(sys.version_info >= (3, 11), reason="requires python3.10 or lower")
 def test_h5m_neutronics_p1():
   o = OMC_DAGMC_harness('examples/step_files/pincell1.step')
   openmc.config['cross_sections']=str(o.nuclear_lib)
   o.tt.results={'keff':(0.07944,0.00070)}
   o.run()
 
-@pytest.mark.skipif(sys.version_info > (3, 10), reason="requires python3.10 or lower")
+#@pytest.mark.skipif(sys.version_info >= (3, 11), reason="requires python3.10 or lower")
 def test_h5m_neutronics_p2():
   o = OMC_DAGMC_harness('examples/step_files/pincell2.step')
   openmc.config['cross_sections']=str(o.nuclear_lib)
   o.tt.results={'keff':(0.07786,0.0008)}
   o.run()
 
-@pytest.mark.skipif(sys.version_info > (3, 10), reason="requires python3.10 or lower")
+#@pytest.mark.skipif(sys.version_info >= (3, 11), reason="requires python3.10 or lower")
 def test_h5m_neutronics_tors():
   o = OMC_DAGMC_harness('examples/step_files/toroids.step')
   # override source spatial distribution
@@ -51,14 +52,14 @@ def test_h5m_neutronics_tors():
   o.tt.results={'keff':(1.61936,0.08562)}
   o.run()
 
-@pytest.mark.skipif(sys.version_info > (3, 10), reason="requires python3.10 or lower")
+#@pytest.mark.skipif(sys.version_info >= (3, 11), reason="requires python3.10 or lower")
 def test_h5m_neutronics_spheroids():
   o = OMC_DAGMC_harness('examples/step_files/spheroids.step')
   openmc.config['cross_sections']=str(o.nuclear_lib)
   o.tt.results={'keff':(1.39082,0.01622)}
   o.run()
 
-@pytest.mark.skipif(sys.version_info > (3, 10), reason="requires python3.10 or lower")
+#@pytest.mark.skipif(sys.version_info >= (3, 10), reason="requires python3.10 or lower")
 def test_h5m_neutronics_ellipsoids():
   o = OMC_DAGMC_harness('examples/step_files/oblate_ellipsoids.step')
   openmc.config['cross_sections']=str(o.nuclear_lib)

tests/test_transform.py

@@ -5,7 +5,7 @@
 
 class HarnessTrans(HarnessRun):
   def __init__(self):
-    self.infile = 'examples/7pin.step'
+    self.infile = 'tests/7pin.step'
     self.a = ab.Assembly(verbose=2)
 
   def run(self,merge=False, **kwargs):