Merge pull request svalinn#192 from svalinn/bugfix/cad-to-dagmc

Fix open-source DAGMC workflow

parastell/parastell.py

@@ -1,10 +1,9 @@
 import argparse
-import yaml
 from pathlib import Path
 
 import cubit
-import numpy as np
 import pystell.read_vmec as read_vmec
+import cadquery as cq
 import cad_to_dagmc
 
 from . import log
@@ -441,21 +440,23 @@ def build_cad_to_dagmc_model(self):
             "Building DAGMC neutronics model via CAD-to-DAGMC..."
         )
 
-        self.dagmc_model = cad_to_dagmc.CadToDagmc()
+        solids = []
+        self._material_tags = []
 
         if self.invessel_build:
-            for solid, mat_tag in zip(
-                *self.invessel_build.extract_solids_and_mat_tags()
-            ):
-                self.dagmc_model.add_cadquery_object(
-                    solid, material_tags=[mat_tag]
-                )
+            ivb_solids, ivb_material_tags = (
+                self.invessel_build.extract_solids_and_mat_tags()
+            )
+            solids.extend(ivb_solids)
+            self._material_tags.extend(ivb_material_tags)
 
         if self.magnet_set:
-            for solid in self.magnet_set.coil_solids:
-                self.dagmc_model.add_cadquery_object(
-                    solid, material_tags=[self.magnet_set.mat_tag]
-                )
+            solids.extend(self.magnet_set.coil_solids)
+            self._material_tags.extend(
+                [self.magnet_set.mat_tag] * len(self.magnet_set.coil_solids)
+            )
+
+        self._geometry = cq.Compound.makeCompound(solids)
 
     def export_cad_to_dagmc(
         self,
@@ -482,11 +483,28 @@ def export_cad_to_dagmc(
         )
 
         export_path = Path(export_dir) / Path(filename).with_suffix(".h5m")
-        self.dagmc_model.export_dagmc_h5m_file(
-            filename=str(export_path),
-            min_mesh_size=min_mesh_size,
-            max_mesh_size=max_mesh_size,
-            imprint=False,
+
+        gmsh_obj = cad_to_dagmc.init_gmsh()
+
+        _, volumes = cad_to_dagmc.get_volumes(
+            gmsh_obj, self._geometry, method="in memory"
+        )
+
+        cad_to_dagmc.mesh_brep(
+            gmsh_obj, min_mesh_size=min_mesh_size, max_mesh_size=max_mesh_size
+        )
+
+        vertices, triangles_by_solid_and_by_face = (
+            cad_to_dagmc.mesh_to_vertices_and_triangles(volumes)
+        )
+
+        gmsh_obj.finalize()
+
+        cad_to_dagmc.vertices_to_h5m(
+            vertices,
+            triangles_by_solid_and_by_face,
+            self._material_tags,
+            h5m_filename=export_path,
         )
 
 

tests/test_parastell.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import dagmc
 
 import parastell.parastell as ps
 from parastell.cubit_io import create_new_cubit_instance
@@ -33,6 +34,13 @@ def remove_files():
         Path.unlink("step_export.log")
 
 
+def check_surfaces_and_volumes(filename, num_surfaces_exp, num_volumes_exp):
+    dagmc_model = dagmc.DAGModel(str(Path(filename).with_suffix(".h5m")))
+
+    assert len(dagmc_model.surfaces) == num_surfaces_exp
+    assert len(dagmc_model.volumes) == num_volumes_exp
+
+
 @pytest.fixture
 def stellarator():
 
@@ -119,6 +127,10 @@ def test_parastell(stellarator):
     component_volume_id_exp = 2
     magnet_volume_ids_exp = list(range(3, 4))
     filename_exp = "dagmc"
+    # Each in-vessel component (2 present) gives 3 unique surfaces; each magnet
+    # (1 present) gives 4 surfaces
+    num_surfaces_exp = 10
+    num_volumes_exp = 3
 
     stellarator.build_cubit_model()
 
@@ -142,12 +154,17 @@ def test_parastell(stellarator):
     assert Path(filename_exp).with_suffix(".h5m").exists()
     assert Path(filename_exp).with_suffix(".cub5").exists()
 
+    check_surfaces_and_volumes(filename_exp, num_surfaces_exp, num_volumes_exp)
+
     remove_files()
 
     stellarator.build_cad_to_dagmc_model()
     stellarator.export_cad_to_dagmc(min_mesh_size=50, max_mesh_size=100)
 
     assert Path(filename_exp).with_suffix(".h5m").exists()
+
+    check_surfaces_and_volumes(filename_exp, num_surfaces_exp, num_volumes_exp)
+
     # Test with custom magnet geometry
 
     create_new_cubit_instance()
@@ -188,11 +205,15 @@ def test_parastell(stellarator):
     assert Path(filename_exp).with_suffix(".h5m").exists()
     assert Path(filename_exp).with_suffix(".cub5").exists()
 
+    check_surfaces_and_volumes(filename_exp, num_surfaces_exp, num_volumes_exp)
+
     remove_files()
 
     stellarator.build_cad_to_dagmc_model()
     stellarator.export_cad_to_dagmc(min_mesh_size=50, max_mesh_size=100)
 
     assert Path(filename_exp).with_suffix(".h5m").exists()
 
+    check_surfaces_and_volumes(filename_exp, num_surfaces_exp, num_volumes_exp)
+
     remove_files()