From 181b6e9341d4818455cb9f9fec8f6308b687ad6b Mon Sep 17 00:00:00 2001 From: Edgar-21 Date: Sat, 16 Nov 2024 15:12:26 -0600 Subject: [PATCH 1/4] building faceted surfaces in cubit working --- parastell/invessel_build.py | 45 +++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/parastell/invessel_build.py b/parastell/invessel_build.py index 078517d..c20d662 100644 --- a/parastell/invessel_build.py +++ b/parastell/invessel_build.py @@ -346,6 +346,51 @@ def export_cad_to_dagmc(self, dagmc_filename="dagmc", export_dir=""): model.export_dagmc_h5m_file(filename=str(export_path)) + def build_volumes_cubit(self): + cubit_io.init_cubit() + for surface_name, surface in self.Surfaces.items(): + print(surface) + vertex_ids = [] + for loop in surface.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) + surface_ids = [] + vertex_ids = np.array(vertex_ids) + for loop_index, loop in enumerate(vertex_ids[0:-1, :]): + for point_index, point in enumerate(loop[0:-1]): + ul = vertex_ids[loop_index, point_index] + ll = vertex_ids[loop_index + 1, point_index] + ur = vertex_ids[loop_index, point_index + 1] + lr = vertex_ids[loop_index + 1, point_index + 1] + cubit.cmd(f"create surface vertex {ul} {ll} {ur}") + surface_ids.append(cubit.get_last_id("surface")) + cubit.cmd(f"create surface vertex {lr} {ll} {ur}") + surface_ids.append(cubit.get_last_id("surface")) + start_cap_ids = " ".join( + [str(vertex) for vertex in vertex_ids[0][0:-1]] + ) + print(start_cap_ids) + end_cap_ids = " ".join( + [str(vertex) for vertex in vertex_ids[-1][0:-1]] + ) + cubit.cmd(f"create surface vertex {start_cap_ids}") + surface_ids.append(cubit.get_last_id("surface")) + cubit.cmd(f"create surface vertex {end_cap_ids}") + surface_ids.append(cubit.get_last_id("surface")) + + surface_id_str = " ".join( + [str(surf_id) for surf_id in surface_ids] + ) + cubit.cmd(f"create volume surface {surface_id_str}") + + cubit.cmd('save cub5 "all_surfaces_faceted.cub5" overwrite') + class Surface(object): """An object representing a surface formed by lofting across a set of From f49277f323d1be90ec60aeb2bef4c88b69e4e76b Mon Sep 17 00:00:00 2001 From: Edgar-21 Date: Sat, 16 Nov 2024 19:32:49 -0600 Subject: [PATCH 2/4] dagmc model export working --- parastell/invessel_build.py | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/parastell/invessel_build.py b/parastell/invessel_build.py index c20d662..6e7f472 100644 --- a/parastell/invessel_build.py +++ b/parastell/invessel_build.py @@ -349,7 +349,7 @@ def export_cad_to_dagmc(self, dagmc_filename="dagmc", export_dir=""): def build_volumes_cubit(self): cubit_io.init_cubit() for surface_name, surface in self.Surfaces.items(): - print(surface) + print(surface_name) vertex_ids = [] for loop in surface.get_loci(): loop_vert_ids = [] @@ -379,17 +379,48 @@ def build_volumes_cubit(self): end_cap_ids = " ".join( [str(vertex) for vertex in vertex_ids[-1][0:-1]] ) + surface_id_str = " ".join( + [str(surf_id) for surf_id in surface_ids] + ) + surface_ids = [] + # for whatever reason, uniting this surface first makes the volume + # creation go much faster (~10x) + cubit.cmd(f"unite surface {surface_id_str}") + surface_ids.append(cubit.get_last_id("surface")) cubit.cmd(f"create surface vertex {start_cap_ids}") surface_ids.append(cubit.get_last_id("surface")) cubit.cmd(f"create surface vertex {end_cap_ids}") surface_ids.append(cubit.get_last_id("surface")) - surface_id_str = " ".join( [str(surf_id) for surf_id in surface_ids] ) cubit.cmd(f"create volume surface {surface_id_str}") - + cubit.cmd("compress") + volume_id = cubit.get_last_id("volume") + self.radial_build.radial_build[surface_name]["vol_id"] = volume_id cubit.cmd('save cub5 "all_surfaces_faceted.cub5" overwrite') + # remove overlap + print("removing overlap") + 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" + ) + cubit.cmd('save cub5 "overlap_removed.cub5" overwrite') + # merge layers + 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']}") + + cubit.cmd('save cub5 "merged.cub5" overwrite') class Surface(object): From dafae7d5ec07385c4370515a25acc36b24d40450 Mon Sep 17 00:00:00 2001 From: Edgar-21 Date: Sun, 17 Nov 2024 23:02:55 -0600 Subject: [PATCH 3/4] re-organize to benefit from oo, integrate with parastell.py --- parastell/invessel_build.py | 171 +++++++++++++++++++----------------- parastell/parastell.py | 14 ++- 2 files changed, 102 insertions(+), 83 deletions(-) diff --git a/parastell/invessel_build.py b/parastell/invessel_build.py index 6e7f472..5a06f48 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,36 @@ 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_step_cubit(self): """Imports STEP files from in-vessel build into Coreform Cubit.""" @@ -346,61 +360,20 @@ def export_cad_to_dagmc(self, dagmc_filename="dagmc", export_dir=""): model.export_dagmc_h5m_file(filename=str(export_path)) - def build_volumes_cubit(self): + 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(): - print(surface_name) - vertex_ids = [] - for loop in surface.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) - surface_ids = [] - vertex_ids = np.array(vertex_ids) - for loop_index, loop in enumerate(vertex_ids[0:-1, :]): - for point_index, point in enumerate(loop[0:-1]): - ul = vertex_ids[loop_index, point_index] - ll = vertex_ids[loop_index + 1, point_index] - ur = vertex_ids[loop_index, point_index + 1] - lr = vertex_ids[loop_index + 1, point_index + 1] - cubit.cmd(f"create surface vertex {ul} {ll} {ur}") - surface_ids.append(cubit.get_last_id("surface")) - cubit.cmd(f"create surface vertex {lr} {ll} {ur}") - surface_ids.append(cubit.get_last_id("surface")) - start_cap_ids = " ".join( - [str(vertex) for vertex in vertex_ids[0][0:-1]] - ) - print(start_cap_ids) - end_cap_ids = " ".join( - [str(vertex) for vertex in vertex_ids[-1][0:-1]] - ) - surface_id_str = " ".join( - [str(surf_id) for surf_id in surface_ids] - ) - surface_ids = [] - # for whatever reason, uniting this surface first makes the volume - # creation go much faster (~10x) - cubit.cmd(f"unite surface {surface_id_str}") - surface_ids.append(cubit.get_last_id("surface")) - cubit.cmd(f"create surface vertex {start_cap_ids}") - surface_ids.append(cubit.get_last_id("surface")) - cubit.cmd(f"create surface vertex {end_cap_ids}") - surface_ids.append(cubit.get_last_id("surface")) - surface_id_str = " ".join( - [str(surf_id) for surf_id in surface_ids] - ) - cubit.cmd(f"create volume surface {surface_id_str}") - cubit.cmd("compress") - volume_id = cubit.get_last_id("volume") - self.radial_build.radial_build[surface_name]["vol_id"] = volume_id - cubit.cmd('save cub5 "all_surfaces_faceted.cub5" overwrite') - # remove overlap - print("removing overlap") + 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:]) @@ -409,18 +382,6 @@ def build_volumes_cubit(self): "remove overlap volume " f"{layer['vol_id']} {next_layer['vol_id']} modify larger" ) - cubit.cmd('save cub5 "overlap_removed.cub5" overwrite') - # merge layers - 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']}") - - cubit.cmd('save cub5 "merged.cub5" overwrite') class Surface(object): @@ -490,6 +451,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..b7be29f 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="" @@ -373,7 +379,7 @@ def build_cubit_model(self, skip_imprint=False, legacy_faceting=True): else: cubit_io.init_cubit() - if self.invessel_build: + if self.invessel_build and not self.invessel_build.cubit_volumes: self.invessel_build.import_step_cubit() if self.magnet_set: @@ -679,7 +685,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"] = {} From 48c02548437060d5014c79ba3d278000e17db5ce Mon Sep 17 00:00:00 2001 From: Edgar-21 Date: Sun, 17 Nov 2024 23:40:38 -0600 Subject: [PATCH 4/4] fix double magnet import, import cub5 to avoid conflicts --- parastell/cubit_io.py | 17 +++++++++++++++++ parastell/invessel_build.py | 22 ++++++++++++++++++++++ parastell/parastell.py | 9 ++++++--- 3 files changed, 45 insertions(+), 3 deletions(-) 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 5a06f48..aaf0e60 100644 --- a/parastell/invessel_build.py +++ b/parastell/invessel_build.py @@ -307,6 +307,16 @@ def merge_layer_surfaces(self): ) 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.""" for name, data in self.radial_build.radial_build.items(): @@ -330,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. diff --git a/parastell/parastell.py b/parastell/parastell.py index b7be29f..faa047a 100644 --- a/parastell/parastell.py +++ b/parastell/parastell.py @@ -190,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( @@ -379,8 +382,8 @@ def build_cubit_model(self, skip_imprint=False, legacy_faceting=True): else: cubit_io.init_cubit() - if self.invessel_build and not self.invessel_build.cubit_volumes: - self.invessel_build.import_step_cubit() + if self.invessel_build: + self.invessel_build.import_geom_cubit() if self.magnet_set: self.magnet_set.import_step_cubit()