Merge pull request #30 from Edgar-21/add_cubit_meshing

.gitignore

@@ -135,3 +135,4 @@ dmypy.json
 *.nc
 *.jou
 *.ods
+*.vtk

ExampleScript.py

@@ -18,7 +18,7 @@
                 [5, 5, 5, 5, 5, 5, 5, 5, 5],
                 [5, 5, 5, 5, 5, 5, 5, 5, 5]
             ]
-        },
+                },
         'breeder': {
             'thickness_matrix': [
                 [100, 100, 30, 10, 10, 10, 30, 100, 100],
@@ -76,9 +76,9 @@
                 [15, 15, 15, 15, 15, 15, 15, 15, 15]
             ],
             'h5m_tag': 'vac_vessel'
-        }
     }
 }
+}
 # Define number of periods in stellarator plasma
 num_periods = 4
 # Define number of periods to generate
@@ -116,6 +116,9 @@
     'facet_tol': 1,
     'len_tol': 5,
     'norm_tol': None,
+    'native_meshing': False, #choose whether to use native cubit meshing v2023.11+ or legacy DAGMC workflow
+    'anisotropic_ratio': 100,
+    'deviation_angle': 5,
     # Note the following export parameters are used only for Gmsh H5M exports
     'min_mesh_size': 5.0,
     'max_mesh_size': 20.0,

parastell.py

@@ -46,6 +46,13 @@ def cubit_export(components, export, magnets):
                     (float, defaults to None),
                 'norm_tol': maximum change in angle between normal vector of
                     adjacent facets (float, defaults to None),
+                'native_meshing': choose native or legacy faceting for DAGMC export
+                    (bool, defaults to True),
+                'anisotropic_ratio': controls edge length ratio of elements
+                    (float, defaults to 100.0),
+                'deviation_angle': controls deviation angle of facet from surface, i.e.
+                    lower deviation angle => more elements in areas with higher curvature
+                    (float, defaults to 5.0),
                 'min_mesh_size': minimum mesh element size for Gmsh export
                     (float, defaults to 5.0),
                 'max_mesh_size': maximum mesh element size for Gmsh export
@@ -76,10 +83,86 @@ def cubit_export(components, export, magnets):
             For a rectangular cross-section, the list format is
             ['rectangle' (str), width (float, cm), thickness (float, cm)]
     """
-    # Extract Cubit export parameters
-    facet_tol = export['facet_tol']
-    len_tol = export['len_tol']
-    norm_tol = export['norm_tol']
+    def legacy_export():
+
+        # Conditionally assign magnet material group
+        if magnets is not None:
+            magnet_h5m_tag = magnets['h5m_tag']
+            cubit.cmd(
+                f'group "mat:{magnet_h5m_tag}" add volume '
+                + " ".join(str(i) for i in magnets['vol_id'])
+            )
+
+        # Assign material groups
+        for comp in components.values():
+            cubit.cmd(f'group "mat:{comp["h5m_tag"]}" add volume {comp["vol_id"]}')
+
+        # Extract Cubit export parameters
+        facet_tol = export['facet_tol']
+        len_tol = export['len_tol']
+        norm_tol = export['norm_tol']
+
+        # Initialize tolerance strings for export statement as empty strings
+        facet_tol_str = ''
+        len_tol_str = ''
+        norm_tol_str = ''
+
+        # Conditionally fill tolerance strings
+        if facet_tol is not None:
+            facet_tol_str = f'faceting_tolerance {facet_tol}'
+        if len_tol is not None:
+            len_tol_str = f'length_tolerance {len_tol}'
+        if norm_tol is not None:
+            norm_tol_str = f'normal_tolerance {norm_tol}'
+
+        # DAGMC export
+        cubit.cmd(
+            f'export dagmc "dagmc.h5m" {facet_tol_str} {len_tol_str} {norm_tol_str}'
+            f' make_watertight'
+        )
+
+    def native_export():
+        #extract Cubit export parameters
+        anisotropic_ratio = export['anisotropic_ratio']
+        deviation_angle = export['deviation_angle']
+
+        # create materials for native cubit meshing
+        for comp in components.values():
+            cubit.cmd(f'create material "{comp["h5m_tag"]}" property_group "CUBIT-ABAQUS"')
+
+        # assign components to blocks
+        for comp in components.values():
+            cubit.cmd('set duplicate block elements off')
+            cubit.cmd("block " + str(comp['vol_id']) + " add volume " + str(comp['vol_id']))
+
+        # assign materials to blocks
+        for comp in components.values():
+            cubit.cmd("block " + str(comp['vol_id']) + " material " + ''.join(("\'",comp['h5m_tag'],"\'")))
+
+
+        if magnets is not None: #conditionally assign material to magnets
+
+            magnet_h5m_tag = magnets['h5m_tag']
+
+            # create magnet material
+            cubit.cmd(f'create material "{magnet_h5m_tag}" property_group "CUBIT-ABAQUS')
+
+            # assign magnets to block
+            block_number = min(magnets['vol_id'])
+            for vol in magnets['vol_id']:
+                cubit.cmd('set duplicate block elements off')
+                cubit.cmd("block " + str(block_number) + " add volume " + str(vol))
+
+            # assign magnet material to block
+            cubit.cmd("block " + str(block_number) + " material " + ''.join(("\'",magnet_h5m_tag,"\'")))
+
+        #mesh the model
+        cubit.cmd("set trimesher coarse on ratio " + str(anisotropic_ratio) + " angle " + str(deviation_angle))
+        cubit.cmd("surface all scheme trimesh")
+        cubit.cmd("mesh surface all")
+
+        #export dagmc file
+        cubit.cmd(f'export cf_dagmc "{cwd + "/dagmc.h5m"}" overwrite')
 
     # Get current working directory
     cwd = os.getcwd()
@@ -93,36 +176,13 @@ def cubit_export(components, export, magnets):
     cubit.cmd('imprint volume all')
     cubit.cmd('merge volume all')
 
-    # Conditionally assign magnet material group
-    if magnets is not None:
-        magnet_h5m_tag = magnets['h5m_tag']
-        cubit.cmd(
-            f'group "mat:{magnet_h5m_tag}" add volume '
-            + " ".join(str(i) for i in magnets['vol_id'])
-        )
-
-    # Assign material groups
-    for comp in components.values():
-        cubit.cmd(f'group "mat:{comp["h5m_tag"]}" add volume {comp["vol_id"]}')
-
-    # Initialize tolerance strings for export statement as empty strings
-    facet_tol_str = ''
-    len_tol_str = ''
-    norm_tol_str = ''
-
-    # Conditionally fill tolerance strings
-    if facet_tol is not None:
-        facet_tol_str = f'faceting_tolerance {facet_tol}'
-    if len_tol is not None:
-        len_tol_str = f'length_tolerance {len_tol}'
-    if norm_tol is not None:
-        norm_tol_str = f'normal_tolerance {norm_tol}'
-
-    # DAGMC export
-    cubit.cmd(
-        f'export dagmc "dagmc.h5m" {facet_tol_str} {len_tol_str} {norm_tol_str}'
-        f' make_watertight'
-    )
+    if export['native_meshing']: #export using native cubit meshing capabilities v2023.11+
+        native_export()
+
+    else: #export with legacy dagmc workflow
+        legacy_export()
+
+
 
 
 def exports(export, components, magnets, logger):
@@ -212,7 +272,7 @@ def exports(export, components, magnets, logger):
             cq.exporters.export(comp['solid'], name + '.step')
 
         # Conditionally export tetrahedral meshing
-        if magnets['meshing']:
+        if magnets is not None and magnets['meshing']:
             # Assign export paths
             cwd = os.getcwd()
             base_name = 'coil_mesh'
@@ -547,6 +607,9 @@ def expand_ang(ang_list, num_ang):
     'facet_tol': None,
     'len_tol': None,
     'norm_tol': None,
+    'native_meshing': True,
+    'anisotropic_ratio': 100,
+    'deviation_angle': 5,
     'min_mesh_size': 5.0,
     'max_mesh_size': 20.0,
     'volume_atol': 0.00001,
@@ -646,6 +709,13 @@ def parastell(
                     (float, defaults to None),
                 'norm_tol': maximum change in angle between normal vector of
                     adjacent facets (float, defaults to None),
+                'native_meshing': choose native or legacy faceting for DAGMC export
+                    (bool, defaults to True),
+                'anisotropic_ratio': controls edge length ratio of elements
+                    (float, defaults to 100.0),
+                'deviation_angle': controls deviation angle of facet from surface, i.e.
+                    lower deviation angle => more elements in areas with higher curvature
+                    (float, defaults to 5.0),
                 'min_mesh_size': minimum mesh element size for Gmsh export
                     (float, defaults to 5.0),
                 'max_mesh_size': maximum mesh element size for Gmsh export