Make Cubit installation optional

README.md

@@ -5,24 +5,23 @@
 [![Build status](https://github.com/svalinn/parastell/actions/workflows/docker_publish.yml/badge.svg)](https://github.com/svalinn/parastell/actions/workflows/docker_publish.yml)
 ---
 
-Open-source Python package featuring a parametric 3-D CAD modeling toolset for stellarator fusion devices with additional neutronics support. ParaStell uses plasma equilibrium VMEC data and a user-defined radial build to model in-vessel components of varying thickness in low-fidelity. Furthermore, coil filament point-locus data and a user-defined cross-section are used to model magnet coils. Additional neutronics support includes the use of VMEC data and MOAB to generate tetrahedral neutron source definitions and Coreform Cubit to generate DAGMC geometries for use in Monte Carlo radiation transport software. In addition, an option is included to generate a tetrahedral mesh of the magnets using Coreform Cubit for use in Monte Carlo mesh tallies. A neutron wall-loading utility is included that uses OpenMC to fire rays from a ParaStell neutron source mesh onto a ParaStell first wall CAD geometry.
+Open-source Python package featuring a parametric 3-D CAD modeling toolset for stellarator fusion devices with additional neutronics support. ParaStell uses plasma equilibrium VMEC data and a user-defined radial build to model in-vessel components of varying thickness in low-fidelity. Furthermore, coil filament point-locus data and a user-defined cross-section are used to model magnet coils. Additional neutronics support includes the generation of tetrahedral neutron source definitions and DAGMC geometries for use in Monte Carlo radiation transport software. In addition, an option is included to generate tetrahedral meshes of in-vessel components and magnets using Coreform Cubit for use in Monte Carlo mesh tallies. A neutron wall-loading utility is included that uses OpenMC to fire rays from a ParaStell neutron source mesh onto a ParaStell first wall CAD geometry.
 
 ![Example model](ParaStell-Example.png)
 
 ## Dependencies
 ParaStell depends on:
 
-- [CadQuery](https://cadquery.readthedocs.io/en/latest/installation.html)
 - [PyStell](https://github.com/aaroncbader/pystell_uw)
+- [CadQuery](https://cadquery.readthedocs.io/en/latest/installation.html)
+- [PyDAGMC](https://github.com/svalinn/pydagmc)
 - [MOAB](https://bitbucket.org/fathomteam/moab/src/master/)
-- [Coreform Cubit](https://coreform.com/products/downloads/)
 - [CAD-to-DAGMC](https://github.com/fusion-energy/cad_to_dagmc)
 - [OpenMC](https://github.com/openmc-dev/openmc)
 - [NumPy](https://numpy.org/install/)
 - [SciPy](https://scipy.org/install/)
-- [scikit-learn](https://scikit-learn.org/stable/install.html)
 - [PyYAML](https://pyyaml.org/wiki/PyYAMLDocumentation)
-- [PyDAGMC](https://github.com/svalinn/pydagmc)
+- [Coreform Cubit](https://coreform.com/products/downloads/) (optional)
 
 ## Install ParaStell
 Download and extract the ParaStell repository:
@@ -52,7 +51,7 @@ conda activate parastell_env
 Alternatively, view `INSTALL.md` for instructions on manually installing these Python dependencies using mamba.
 
 ### Install Coreform Cubit
-Download and install the latest version from [Coreform's Website](https://coreform.com/products/downloads/), then add the `/Coreform-Cubit-[version]/bin/` directory to your `PYTHONPATH` by adding a line similar to the following to your `.bashrc` file:
+To make use of ParaStell's Cubit functionality, download and install the latest version from [Coreform's Website](https://coreform.com/products/downloads/), then add the `/Coreform-Cubit-[version]/bin/` directory to your `PYTHONPATH` by adding a line similar to the following to your `.bashrc` file:
 
 ```bash
 export PYTHONPATH=$PYTHONPATH:$HOME/Coreform-Cubit-[version]/bin/

environment.yml

@@ -6,7 +6,6 @@ dependencies:
   - pip
   - numpy=1.26.4
   - scipy
-  - scikit-learn
   - cadquery
   - moab=5.5.1[build=nompi_tempest_*]
   - cad_to_dagmc >=0.8.2

parastell/cubit_io.py → parastell/cubit_utils.py

@@ -1,13 +1,14 @@
 from pathlib import Path
 import subprocess
 
-import cubit
-
 initialized = False
 
 
 def init_cubit():
     """Initializes Coreform Cubit with the DAGMC plugin."""
+    global cubit
+    import cubit
+
     global initialized
 
     if not initialized:
@@ -63,8 +64,6 @@ def export_step_cubit(filename, export_dir=""):
         export_dir (str): directory to which to export the STEP output file
             (defaults to empty string).
     """
-    init_cubit()
-
     export_path = Path(export_dir) / Path(filename).with_suffix(".step")
     cubit.cmd(f'export step "{export_path}" overwrite')
 
@@ -94,8 +93,6 @@ def export_cub5(filename, export_dir=""):
         export_dir (str): directory to which to export the cub5 output file
             (defaults to empty string).
     """
-    init_cubit()
-
     export_path = Path(export_dir) / Path(filename).with_suffix(".cub5")
     cubit.cmd(f'save cub5 "{export_path}" overwrite')
 
@@ -111,8 +108,6 @@ def export_mesh_cubit(filename, export_dir="", delete_upon_export=True):
         delete_upon_export (bool): delete the mesh from the Cubit instance
             after exporting. Prevents inclusion of mesh in future exports.
     """
-    init_cubit()
-
     exo_path = Path(export_dir) / Path(filename).with_suffix(".exo")
     h5m_path = Path(export_dir) / Path(filename).with_suffix(".h5m")
 
@@ -162,8 +157,6 @@ def export_dagmc_cubit(
         delete_upon_export (bool): delete the mesh from the Cubit instance
             after exporting. Prevents inclusion of mesh in future exports.
     """
-    init_cubit()
-
     cubit.cmd(
         f"set trimesher coarse on ratio {anisotropic_ratio} "
         f"angle {deviation_angle}"
@@ -196,3 +189,120 @@ def import_geom_to_cubit(filename, import_dir=""):
     filename = Path(filename)
     vol_id = importers[filename.suffix](filename, import_dir)
     return vol_id
+
+
+def make_material_block(mat_tag, block_id, vol_id_str):
+    """Issue commands to make a material block in Cubit.
+
+    Arguments:
+       mat_tag (str) : name of material block
+       block_id (int) : block number
+       vol_id_str (str) : space-separated list of volume ids
+    """
+    cubit.cmd("set duplicate block elements off")
+    cubit.cmd(f'create material "{mat_tag}" property_group ' '"CUBIT-ABAQUS"')
+    cubit.cmd(f"block {block_id} add volume {vol_id_str}")
+    cubit.cmd(f'block {block_id} material "{mat_tag}"')
+
+
+def imprint_and_merge():
+    """Imprints and merges all volumes in Cubit."""
+    cubit.cmd("imprint volume all")
+    cubit.cmd("merge volume all")
+
+
+def orient_spline_surfaces(volume_id):
+    """Extracts the inner and outer surface IDs for a given ParaStell in-vessel
+    component volume in Coreform Cubit.
+
+    Arguments:
+        volume_id (int): Cubit volume ID.
+
+    Returns:
+        inner_surface_id (int): Cubit ID of in-vessel component inner surface.
+        outer_surface_id (int): Cubit ID of in-vessel component outer surface.
+    """
+
+    surfaces = cubit.get_relatives("volume", volume_id, "surface")
+
+    spline_surfaces = []
+    for surface in surfaces:
+        if cubit.get_surface_type(surface) == "spline surface":
+            spline_surfaces.append(surface)
+
+    if len(spline_surfaces) == 1:
+        outer_surface_id = spline_surfaces[0]
+        inner_surface_id = None
+    else:
+        # The outer surface bounding box will have the larger maximum XY value
+        if (
+            cubit.get_bounding_box("surface", spline_surfaces[1])[4]
+            > cubit.get_bounding_box("surface", spline_surfaces[0])[4]
+        ):
+            outer_surface_id = spline_surfaces[1]
+            inner_surface_id = spline_surfaces[0]
+        else:
+            outer_surface_id = spline_surfaces[0]
+            inner_surface_id = spline_surfaces[1]
+
+    return inner_surface_id, outer_surface_id
+
+
+def merge_surfaces(surface_1, surface_2):
+    """Merges two surfaces in Cubit.
+
+    Arguments:
+        surface_1 (int): Cubit ID of one surface to be merged.
+        surface_2 (int): Cubit ID of the other surface to be merged.
+    """
+    cubit.cmd(f"merge surface {surface_1} {surface_2}")
+
+
+def mesh_volume_auto_factor(volume_ids, mesh_size=5.0):
+    """Meshes a volume in Cubit using automatically calculated interval sizes.
+
+    Arguments:
+        volume_ids (iterable of int): Cubit IDs of volumes to be meshed.
+        mesh_size (float): controls the size of the mesh. Takes values between
+            1.0 (finer) and 10.0 (coarser) (optional, defaults to 5.0).
+    """
+    volume_ids_str = " ".join(str(id) for id in volume_ids)
+
+    cubit.cmd(f"volume {volume_ids_str} scheme tetmesh")
+    cubit.cmd(f"volume {volume_ids_str} size auto factor {mesh_size}")
+    cubit.cmd(f"mesh volume {volume_ids_str}")
+
+
+def mesh_volume_skeleton(
+    volume_ids, min_size=20.0, max_size=50.0, max_gradient=1.5
+):
+    """Meshes a volume in Cubit using the skeleton sizing function.
+
+    Arguments:
+        volume_ids (iterable of int): Cubit IDs of volumes to be meshed.
+        min_size (float): minimum size of mesh elements (defaults to 20.0).
+        max_size (float): maximum size of mesh elements (defaults to 50.0).
+        max_gradient (float): maximum transition in mesh element size
+            (defaults to 1.5).
+    """
+    volume_ids_str = " ".join(str(id) for id in volume_ids)
+
+    cubit.cmd(f"volume {volume_ids_str} scheme tetmesh")
+    cubit.cmd(
+        f"volume {volume_ids_str} sizing function type skeleton min_size "
+        f"{min_size} max_size {max_size} max_gradient {max_gradient} "
+        "min_num_layers_3d 1 min_num_layers_2d 1 min_num_layers_1d 1"
+    )
+    cubit.cmd(f"mesh volume {volume_ids_str}")
+
+
+def get_last_id(entity):
+    """Returns the ID of the most recently created entity of a given type.
+
+    Arguments:
+        entity (str): the type of entity for which the ID should be retrieved.
+
+    Returns:
+        (int): the ID of the most recently created entity of the given type.
+    """
+    return cubit.get_last_id(entity)

parastell/invessel_build.py

@@ -4,14 +4,19 @@
 import numpy as np
 from scipy.interpolate import RegularGridInterpolator
 
-import cubit
 import cadquery as cq
 import pystell.read_vmec as read_vmec
 import dagmc
 from pymoab import core, types
 
 from . import log
-from . import cubit_io
+from .cubit_utils import (
+    import_step_cubit,
+    export_mesh_cubit,
+    orient_spline_surfaces,
+    merge_surfaces,
+    mesh_volume_auto_factor,
+)
 from .utils import (
     normalize,
     expand_list,
@@ -49,43 +54,6 @@ def create_moab_tris_from_verts(corners, mbc, reverse=False):
     return [tri_1, tri_2]
 
 
-def orient_spline_surfaces(volume_id):
-    """Extracts the inner and outer surface IDs for a given ParaStell in-vessel
-    component volume in Coreform Cubit.
-
-    Arguments:
-        volume_id (int): Cubit volume ID.
-
-    Returns:
-        inner_surface_id (int): Cubit ID of in-vessel component inner surface.
-        outer_surface_id (int): Cubit ID of in-vessel component outer surface.
-    """
-
-    surfaces = cubit.get_relatives("volume", volume_id, "surface")
-
-    spline_surfaces = []
-    for surface in surfaces:
-        if cubit.get_surface_type(surface) == "spline surface":
-            spline_surfaces.append(surface)
-
-    if len(spline_surfaces) == 1:
-        outer_surface_id = spline_surfaces[0]
-        inner_surface_id = None
-    else:
-        # The outer surface bounding box will have the larger maximum XY value
-        if (
-            cubit.get_bounding_box("surface", spline_surfaces[1])[4]
-            > cubit.get_bounding_box("surface", spline_surfaces[0])[4]
-        ):
-            outer_surface_id = spline_surfaces[1]
-            inner_surface_id = spline_surfaces[0]
-        else:
-            outer_surface_id = spline_surfaces[0]
-            inner_surface_id = spline_surfaces[1]
-
-    return inner_surface_id, outer_surface_id
-
-
 class InVesselBuild(object):
     """Parametrically models fusion stellarator in-vessel components using
     plasma equilibrium VMEC data and a user-defined radial build.
@@ -490,15 +458,13 @@ def merge_layer_surfaces(self):
             if prev_outer_surface_id is None:
                 prev_outer_surface_id = outer_surface_id
             else:
-                cubit.cmd(
-                    f"merge surface {inner_surface_id} {prev_outer_surface_id}"
-                )
+                merge_surfaces(inner_surface_id, prev_outer_surface_id)
                 prev_outer_surface_id = outer_surface_id
 
     def import_step_cubit(self):
         """Imports STEP files from in-vessel build into Coreform Cubit."""
         for name, data in self.radial_build.radial_build.items():
-            vol_id = cubit_io.import_step_cubit(name, self.export_dir)
+            vol_id = import_step_cubit(name, self.export_dir)
             data["vol_id"] = vol_id
 
     def export_step(self, export_dir=""):
@@ -544,19 +510,18 @@ def export_component_mesh(
         Arguments:
             components (array of strings): array containing the name
                 of the in-vessel components to be meshed.
-            mesh_size (int): controls the size of the mesh. Takes values
-                between 1 (finer) and 10 (coarser) (optional, defaults to 5).
+            mesh_size (float): controls the size of the mesh. Takes values
+                between 1.0 (finer) and 10.0 (coarser) (optional, defaults to
+                5.0).
             import_dir (str): directory containing the STEP file of
                 the in-vessel component (optional, defaults to empty string).
             export_dir (str): directory to which to export the h5m
                 output file (optional, defaults to empty string).
         """
         for component in components:
-            vol_id = cubit_io.import_step_cubit(component, import_dir)
-            cubit.cmd(f"volume {vol_id} scheme tetmesh")
-            cubit.cmd(f"volume {vol_id} size auto factor {mesh_size}")
-            cubit.cmd(f"mesh volume {vol_id}")
-            cubit_io.export_mesh_cubit(
+            vol_id = import_step_cubit(component, import_dir)
+            mesh_volume_auto_factor([vol_id], mesh_size=mesh_size)
+            export_mesh_cubit(
                 filename=component,
                 export_dir=export_dir,
                 delete_upon_export=True,