Skip imprinting

.gitignore

@@ -136,3 +136,6 @@ dmypy.json
 *.jou
 *.ods
 *.vtk
+*.cub
+*.cub5
+*.h5

ExampleScript.py

@@ -94,6 +94,7 @@
     'facet_tol': 1,
     'len_tol': 5,
     'norm_tol': None,
+    'skip_imprinting': False,
     # Choose whether to use native Cubit meshing (v2023.11+) or legacy DAGMC
     # workflow
     'native_meshing': False,
@@ -111,4 +112,4 @@
 strengths = parastell.parastell(
     plas_eq, build, repeat, num_phi, num_theta,
     magnets = magnets, source = source, export = export
-)
+)

magnet_coils.py

@@ -538,7 +538,7 @@ def magnet_coils(magnets, tor_ext, export_dir, logger = None):
         vol_ids = cut_mags(tor_ext, vol_ids, r_avg)
 
     # Export magnet coils
-    export_path = Path(export_dir) / Path(magnets["name"]).with_suffix('step')
+    export_path = Path(export_dir) / Path(magnets["name"]).with_suffix('.step')
     cubit.cmd(
         f'export step "{export_path}" overwrite'
     )

parastell.py

@@ -92,6 +92,7 @@ def cubit_export(components, export, magnets):
             For a rectangular cross-section, the list format is
             ['rectangle' (str), width (float, cm), thickness (float, cm)]
     """
+
     def legacy_export():
         """Exports neutronics H5M file via legacy plug-in faceting method.
         """
@@ -185,7 +186,7 @@ def native_export():
                 "block " + str(block_number) + " material "
                 + ''.join(("\'",magnet_h5m_tag,"\'"))
             )
-        
+
         # Mesh the model
         cubit.cmd(
             "set trimesher coarse on ratio " + str(anisotropic_ratio)
@@ -203,19 +204,89 @@ def native_export():
     dir = Path(export['dir'])
     filename = Path(export['h5m_filename'])
 
+    def orient_spline_surfaces(volume_id):
+        """Return the inner and outer surface ids for a given volume id
+        """
+        surfaces = cubit.get_relatives('volume', volume_id, 'surface')
+
+        splineSurfaces = []
+
+        for surf in surfaces:
+            if cubit.get_surface_type(surf) == 'spline surface':
+                splineSurfaces.append(surf)
+
+        # check if this is the plasma
+        if len(splineSurfaces) == 1:
+            outerSurface = splineSurfaces[0]
+            innerSurface = None
+
+        else:
+            # the outer surface bounding box will have the larger max xy value
+            if cubit.get_bounding_box('surface', splineSurfaces[1])[4] > cubit.get_bounding_box('surface',splineSurfaces[0])[4]:
+                outerSurface = splineSurfaces[1]
+                innerSurface = splineSurfaces[0]
+
+            else:
+                outerSurface = splineSurfaces[0]
+                innerSurface = splineSurfaces[1]
+
+        return innerSurface, outerSurface
+
+    def merge_layer_surfaces(): # assumes that components dict is ordered from inside to out
+        """Merge surfaces based on surface ids rather than imprinting/merging all
+        """
+
+        # tracks the surface id of the outer surface of the previous layer
+        lastOuterSurface = None
+
+        for name in components.keys():
+
+            # get volume ID for layer
+
+            vol_id = components[name]['vol_id']
+
+            # get the inner and outer surface IDs of the current layer
+            innerSurface, outerSurface = orient_spline_surfaces(vol_id)
+
+            # wait to merge until the next layer if the plasma is included
+            # store surface to be merged for next loop
+            if name == 'plasma':
+
+                lastOuterSurface = outerSurface
+
+            # check if we are in the first layer in a build with plasma excluded
+            # store outer surface to be merged in next loop
+            elif lastOuterSurface is None:
+
+                lastOuterSurface = outerSurface
+
+            # merge inner surface with outer surface of previous layer
+            else:
+
+                cubit.cmd('merge surface '+ str(innerSurface) + ' ' + str(lastOuterSurface))
+
+                lastOuterSurface = outerSurface
+
+
     for name in components.keys():
         import_path = dir / Path(name).with_suffix('.step')
         cubit.cmd(f'import step "{import_path}" heal')
         components[name]['vol_id'] = cubit.get_last_id("volume")
 
-    cubit.cmd('imprint volume all')
-    cubit.cmd('merge volume all')
+    if export['skip_imprinting']:
+        merge_layer_surfaces()
 
+    else:
+        cubit.cmd('imprint volume all')
+        cubit.cmd('merge volume all')
+
     if export['native_meshing']:
         native_export()
+
     else:
         legacy_export()
-
+
+
 
 def exports(export, components, magnets, logger):
     """Export components.
@@ -570,6 +641,7 @@ def expand_ang(ang_list, num_ang):
     'facet_tol': None,
     'len_tol': None,
     'norm_tol': None,
+    'skip_imprinting': False,
     'anisotropic_ratio': 100,
     'deviation_angle': 5,
     'min_mesh_size': 5.0,
@@ -683,6 +755,9 @@ def parastell(
                     (float, defaults to None),
                 'norm_tol': maximum change in angle between normal vector of
                     adjacent facets (float, defaults to None),
+                'skip_imprinting': choose whether to imprint and merge all in
+                    cubit or to merge surfaces based on import order and
+                    geometry information.
                 'anisotropic_ratio': controls edge length ratio of elements
                     (float, defaults to 100.0),
                 'deviation_angle': controls deviation angle of facet from
@@ -832,6 +907,7 @@ def parastell(
                 vmec, s, tor_ext, repeat, phi_list_exp, theta_list_exp, interp,
                 cutter
             )
+
         except ValueError as e:
             logger.error(e.args[0])
             raise e
@@ -847,7 +923,7 @@ def parastell(
 
     if magnets is not None:
         magnets['vol_id'] = magnet_coils.magnet_coils(
-            magnets, (repeat + 1)*tor_ext, export['dir'], logger = logger
+            magnets, (repeat + 1)*tor_ext, export_dict['dir'], logger = logger
         )
 
     try:
@@ -858,6 +934,12 @@ def parastell(
 
     if source is not None:
         strengths = source_mesh.source_mesh(
-            vmec, source, export['dir'], logger = logger
+            vmec, source, export_dict['dir'], logger = logger
         )
         return strengths
+
+    # reset cubit to avoid issues when looping parastell
+    if export_dict['h5m_export'] == 'Cubit':
+        cubit.cmd('reset')
+
+