diff --git a/parastell/cubit_io.py b/parastell/cubit_io.py index 4a984ad..602fef0 100644 --- a/parastell/cubit_io.py +++ b/parastell/cubit_io.py @@ -51,6 +51,23 @@ def import_step_cubit(filename, import_dir): return vol_id +def import_cub5_cubit(filename, import_dir): + """Imports cub5 file with Coreform Cubit with default import settings. + Arguments: + filename (str): name of cub5 input file. + import_dir (str): directory from which to import cub5 file. + Returns: + vol_id (int): Cubit volume ID of imported CAD solid. + """ + init_cubit() + import_path = Path(import_dir) / Path(filename).with_suffix(".cub5") + cubit.cmd( + f'import cubit "{import_path}" nofreesurfaces attributes_on separate_bodies' + ) + vol_id = cubit.get_last_id("volume") + return vol_id + + def export_step_cubit(filename, export_dir=""): """Export CAD solid as a STEP file via Coreform Cubit. diff --git a/parastell/invessel_build.py b/parastell/invessel_build.py index 078517d..aaf0e60 100644 --- a/parastell/invessel_build.py +++ b/parastell/invessel_build.py @@ -231,6 +231,7 @@ def generate_components(self): build by cutting the interior surface solid from the outer surface solid for a given component. """ + self.cubit_volumes = False self._logger.info( "Constructing CadQuery objects for in-vessel components..." ) @@ -275,23 +276,46 @@ def merge_layer_surfaces(self): overlaps between magnet volumes and in-vessel components will not be merged in this workflow. """ - # Tracks the surface id of the outer surface of the previous layer - prev_outer_surface_id = None - - for data in self.radial_build.radial_build.values(): + if self.cubit_volumes: + layers = list(self.radial_build.radial_build.values()) + for layer, next_layer in zip( + reversed(layers[0:-1]), reversed(layers[1:]) + ): + print(layer) + cubit.cmd( + f"imprint volume {layer['vol_id']} {next_layer['vol_id']}" + ) + cubit.cmd( + f"merge volume {layer['vol_id']} {next_layer['vol_id']}" + ) + else: + # Tracks the surface id of the outer surface of the previous layer + prev_outer_surface_id = None - inner_surface_id, outer_surface_id = orient_spline_surfaces( - data["vol_id"] - ) + for data in self.radial_build.radial_build.values(): - # Conditionally skip merging (first iteration only) - 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}" + inner_surface_id, outer_surface_id = orient_spline_surfaces( + data["vol_id"] ) - prev_outer_surface_id = outer_surface_id + + # Conditionally skip merging (first iteration only) + 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}" + ) + prev_outer_surface_id = outer_surface_id + + def import_geom_cubit(self): + if self.cubit_volumes: + self.import_cub5_cubit() + else: + self.import_step_cubit() + + def import_cub5_cubit(self): + """Import cub5 files from in-vessel build into Coreform Cubit""" + cubit_io.import_cub5_cubit("invessel_build.cub5", self.export_dir) def import_step_cubit(self): """Imports STEP files from in-vessel build into Coreform Cubit.""" @@ -316,6 +340,18 @@ def export_step(self, export_dir=""): ) cq.exporters.export(component, str(export_path)) + def save_cub5(self, filename="invessel_build.cub5", export_dir=""): + """Save current state of cubit model. + + Arguments: + filename (str): name to save file as. + export_dir (str): directory in which to save the cub5 file. + """ + self.export_dir = export_dir + cubit.cmd( + f'save cub5 "{Path(self.export_dir) / Path(filename).with_suffix(".cub5")}" overwrite' + ) + def export_cad_to_dagmc(self, dagmc_filename="dagmc", export_dir=""): """Exports DAGMC neutronics H5M file of ParaStell in-vessel components via CAD-to-DAGMC. @@ -346,6 +382,29 @@ def export_cad_to_dagmc(self, dagmc_filename="dagmc", export_dir=""): model.export_dagmc_h5m_file(filename=str(export_path)) + def build_cubit_volumes(self): + self.cubit_volumes = True + cubit_io.init_cubit() + self._logger.info("Creating faceted volumes in Coreform Cubit...") + for surface_name, surface in self.Surfaces.items(): + surface.build_cubit_vertices_from_loci() + surface.build_cubit_surface() + # TODO handle when full device is specified, rather than one period + surface.close_cubit_surface_ends() + surface.build_cubit_volume() + self.radial_build.radial_build[surface_name][ + "vol_id" + ] = surface.volume_id + + layers = list(self.radial_build.radial_build.values()) + for layer, next_layer in zip( + reversed(layers[0:-1]), reversed(layers[1:]) + ): + cubit.cmd( + "remove overlap volume " + f"{layer['vol_id']} {next_layer['vol_id']} modify larger" + ) + class Surface(object): """An object representing a surface formed by lofting across a set of @@ -414,6 +473,58 @@ def get_loci(self): """Returns the set of point-loci defining the ribs in the surface.""" return np.array([rib.rib_loci for rib in self.Ribs]) + def build_cubit_vertices_from_loci(self): + vertex_ids = [] + for loop in self.get_loci(): + loop_vert_ids = [] + for point in loop: + cubit.cmd(f"create vertex {point[0]} {point[1]} {point[2]}") + vertex_id = cubit.get_last_id("vertex") + loop_vert_ids.append(vertex_id) + vertex_ids.append(loop_vert_ids) + self.vertex_ids = np.array(vertex_ids) + + def build_cubit_surface(self): + self.surface_ids = [] + for loop_index, loop in enumerate(self.vertex_ids[0:-1, :]): + for point_index, _ in enumerate(loop[0:-1]): + ul = self.vertex_ids[loop_index, point_index] + ll = self.vertex_ids[loop_index + 1, point_index] + ur = self.vertex_ids[loop_index, point_index + 1] + lr = self.vertex_ids[loop_index + 1, point_index + 1] + cubit.cmd(f"create surface vertex {ul} {ll} {ur}") + self.surface_ids.append(cubit.get_last_id("surface")) + cubit.cmd(f"create surface vertex {lr} {ll} {ur}") + self.surface_ids.append(cubit.get_last_id("surface")) + surface_id_str = " ".join( + [str(surf_id) for surf_id in self.surface_ids] + ) + cubit.cmd(f"unite surface {surface_id_str}") + self.surface_ids = [cubit.get_last_id("surface")] + + def close_cubit_surface_ends(self): + start_cap_ids = " ".join( + [str(vertex) for vertex in self.vertex_ids[0][0:-1]] + ) + end_cap_ids = " ".join( + [str(vertex) for vertex in self.vertex_ids[-1][0:-1]] + ) + cubit.cmd(f"create surface vertex {start_cap_ids}") + self.surface_ids.append(cubit.get_last_id("surface")) + cubit.cmd(f"create surface vertex {end_cap_ids}") + self.surface_ids.append(cubit.get_last_id("surface")) + + def build_cubit_volume(self): + surface_id_str = " ".join( + [str(surf_id) for surf_id in self.surface_ids] + ) + cubit.cmd(f"create volume surface {surface_id_str}") + # When creating a surface, it actually makes a weird volume, so when + # combining surfaces into a volume, naturally it leaves gaps in the + # id space. Compressing here gives consecutively id'ed volumes + cubit.cmd("compress") + self.volume_id = cubit.get_last_id("volume") + class Rib(object): """An object representing a curve formed by interpolating a spline through diff --git a/parastell/parastell.py b/parastell/parastell.py index e73628d..faa047a 100644 --- a/parastell/parastell.py +++ b/parastell/parastell.py @@ -91,6 +91,7 @@ def construct_invessel_build( wall_s, radial_build, split_chamber=False, + cadquery_geometry=True, **kwargs, ): """Construct InVesselBuild class object. @@ -129,7 +130,9 @@ def construct_invessel_build( scrape-off layer definition for 'chamber', add an item with a 'chamber' key and desired 'thickness_matrix' value to the radial_build dictionary. - + cadquery_geometry (bool): If True, use CadQuery to build spline + based geometry. If False, use Coreform Cubit to build facet + based geometry. Optional attributes: plasma_mat_tag (str): alternate DAGMC material tag to use for plasma. If none is supplied, 'Vacuum' will be used (defaults to @@ -167,7 +170,10 @@ def construct_invessel_build( self.invessel_build.populate_surfaces() self.invessel_build.calculate_loci() - self.invessel_build.generate_components() + if cadquery_geometry: + self.invessel_build.generate_components() + else: + self.invessel_build.build_cubit_volumes() def export_invessel_build( self, export_cad_to_dagmc=False, dagmc_filename="dagmc", export_dir="" @@ -184,7 +190,10 @@ def export_invessel_build( export_dir (str): directory to which to export the output files (optional, defaults to empty string). """ - self.invessel_build.export_step(export_dir=export_dir) + if not self.invessel_build.cubit_volumes: + self.invessel_build.export_step(export_dir=export_dir) + else: + self.invessel_build.save_cub5(export_dir=export_dir) if export_cad_to_dagmc: self.invessel_build.export_cad_to_dagmc( @@ -374,7 +383,7 @@ def build_cubit_model(self, skip_imprint=False, legacy_faceting=True): cubit_io.init_cubit() if self.invessel_build: - self.invessel_build.import_step_cubit() + self.invessel_build.import_geom_cubit() if self.magnet_set: self.magnet_set.import_step_cubit() @@ -679,7 +688,7 @@ def parastell(): nwl_required_keys = ["toroidal_angles", "poloidal_angles", "wall_s"] nwl_build = {} - for key in nwl_keys: + for key in nwl_required_keys: nwl_build[key] = invessel_build[key] nwl_build["radial_build"] = {}