svalinn · gonuke · Mar 21, 2023 · Feb 11, 2023 · Feb 12, 2023 · Feb 13, 2023
diff --git a/doc/CHANGELOG.rst b/doc/CHANGELOG.rst
@@ -9,6 +9,7 @@ Next version
 
 **Changed:**
 
+   * Improvements/corrections to graveyard capabilities (#855)
    * Using multi stage Dockerfile to reduce the number of Dockerfile (#813)
    * Adding safe folder to allow CI to compile DAGMC (#814)
    * Correction to CMake variable name in OpenMC install instructions (#817)

diff --git a/src/dagmc/DagMC.cpp b/src/dagmc/DagMC.cpp
@@ -261,6 +261,8 @@ ErrorCode DagMC::get_graveyard_group(EntityHandle& graveyard_group) {
 }
 
 ErrorCode DagMC::remove_graveyard() {
+  if (!has_graveyard()) return MB_SUCCESS;
+
   ErrorCode rval;
 
   EntityHandle graveyard_group;
@@ -274,12 +276,8 @@ ErrorCode DagMC::remove_graveyard() {
   Range sets_to_delete, ents_to_delete, verts_to_delete;
   sets_to_delete.insert(graveyard_group);
 
-  // check for bounding box trees on the model
-#ifdef DOUBLE_DOWN
-  bool trees_exist = ray_tracer->has_bvh();
-#else
-  bool trees_exist = geom_tool()->have_obb_tree();
-#endif
+  bool trees_exist = has_acceleration_datastructures();
+
   // get the graveyard volume
   Range graveyard_vols;
   rval =
@@ -290,7 +288,7 @@ ErrorCode DagMC::remove_graveyard() {
   // get the implicit complement, it's children will need updating
   EntityHandle implicit_complement = 0;
   rval = geom_tool()->get_implicit_complement(implicit_complement);
-  if (rval != MB_ENTITY_NOT_FOUND || rval != MB_SUCCESS) {
+  if (rval != MB_ENTITY_NOT_FOUND && rval != MB_SUCCESS) {
     MB_CHK_SET_ERR(rval, "Could not get the implicit complement");
   }
 
@@ -391,33 +389,52 @@ ErrorCode DagMC::create_graveyard(bool overwrite) {
   */
   ErrorCode rval;
 
-  // remove existing graveyard
+  // remove existing graveyard if overwrite is true
   if (overwrite) {
     remove_graveyard();
   }
 
-  // if a graveyard already exists at this point,
+  // if a graveyard already exists and we aren't overwriting it,
   // report an error
   if (has_graveyard()) {
     MB_CHK_SET_ERR(MB_FAILURE, "Graveyard already exists");
   }
 
-  // currently relying on the BVH as looping over all vertices may
-  // be too expensive
-  if (!has_acceleration_datastructures()) {
-    MB_CHK_SET_ERR(MB_FAILURE, "Graveyard creation attempted without BVH");
-  }
-
+  // create a bounding box for all volumes
   BBOX box;
-  for (int i = 0; i < num_entities(3); i++) {
-    // get the bounding box of the volume
-    moab::EntityHandle vol = this->entity_by_index(3, i + 1);
-    double vmin[3], vmax[3];
-    rval = this->getobb(vol, vmin, vmax);  // this method name is a misnomer
-    MB_CHK_SET_ERR(rval, "Failed to get volume bounding box");
-    // update the global bounding box
-    box.update(vmin);
-    box.update(vmax);
+
+  // if there are no acceleration data structures present, build
+  // the bounding box using the vertex coordinates
+  if (!has_acceleration_datastructures()) {
+    // get the vertices of every surface
+    for (int i = 0; i < num_entities(2); i++) {
+      // get the bounding box of the volume
+      moab::EntityHandle surf = this->entity_by_index(2, i + 1);
+      moab::Range vertices;
+      rval = this->moab_instance()->get_entities_by_type(surf, moab::MBVERTEX,
+                                                         vertices);
+      MB_CHK_SET_ERR(rval, "Failed to get surface vertices");
+      double coords[3];
+      for (auto vertex : vertices) {
+        // get each vertex coordinate and update box
+        rval = this->moab_instance()->get_coords(&vertex, 1, coords);
+        MB_CHK_SET_ERR(rval, "Failed to get vertex coordinates");
+        box.update(coords);
+      }
+    }
+    // if there acceleration data structures exist, use those for
+    // a faster bounding box build
+  } else {
+    for (int i = 0; i < num_entities(3); i++) {
+      // get the bounding box of the volume
+      moab::EntityHandle vol = this->entity_by_index(3, i + 1);
+      double vmin[3], vmax[3];
+      rval = this->getobb(vol, vmin, vmax);  // this method name is a misnomer
+      MB_CHK_SET_ERR(rval, "Failed to get volume bounding box");
+      // update the global bounding box
+      box.update(vmin);
+      box.update(vmax);
+    }
   }
 
   if (!box.valid()) {
@@ -471,9 +488,21 @@ ErrorCode DagMC::create_graveyard(bool overwrite) {
   box.expand(10.0 * numerical_precision());
 
   // tear down the implicit complement tree
-  EntityHandle implicit_complement;
+  EntityHandle implicit_complement = 0;
   rval = geom_tool()->get_implicit_complement(implicit_complement);
-  MB_CHK_SET_ERR(rval, "Failed to get the implicit complement");
+  if (rval != MB_ENTITY_NOT_FOUND && rval != MB_SUCCESS) {
+    MB_CHK_SET_ERR(rval, "Could not get the implicit complement");
+  }
+  // create the implicit complement if it doesn't exist at this point
+  // the code below that inserts the graveyard into the implicit complement can
+  // be run without changing the model at this point
+  if (!implicit_complement) {
+    rval = setup_impl_compl();
+    MB_CHK_SET_ERR(rval, "Failed to create the implicit complement.");
+    rval = geom_tool()->get_implicit_complement(implicit_complement);
+    MB_CHK_SET_ERR(rval,
+                   "Failed to get implicit complement right after creation");
+  }
 
   EntityHandle inner_surface;
   rval = box_to_surf(box.lower, box.upper, inner_surface);
@@ -515,7 +544,8 @@ ErrorCode DagMC::create_graveyard(bool overwrite) {
                  "Failed to create the graveyard parent-child relationship");
 
   // set the surface senses (all triangles have outward normals so this should
-  // be FORWARD wrt the graveyard volume)
+  // be FORWARD wrt the graveyard volume and REVERSE wrt the implicit
+  // complement)
   EntityHandle outer_senses[2] = {volume_set, implicit_complement};
   rval = MBI->tag_set_data(sense_tag(), &outer_surface, 1, outer_senses);
   MB_CHK_SET_ERR(rval, "Failed to set graveyard surface senses");
@@ -526,21 +556,24 @@ ErrorCode DagMC::create_graveyard(bool overwrite) {
   rval = geom_tool()->find_geomsets();
   MB_CHK_SET_ERR(rval, "Failed to update the geometry sets");
 
-  // delete the implicit complement tree (but not the surface trees)
-  rval = remove_bvh(implicit_complement, true);
-  MB_CHK_SET_ERR(rval, "Failed to delete the implicit complement tree");
-
   // create BVH for both the new implicit complement and the new graveyard
   // volume
-  rval = build_bvh(volume_set);
-  MB_CHK_SET_ERR(
-      rval,
-      "Failed to build accel. data structure for the new graveyard volume");
-
-  rval = build_bvh(implicit_complement);
-  MB_CHK_SET_ERR(
-      rval,
-      "Failed to build accel. data structure for the new implicit complement");
+  if (has_acceleration_datastructures()) {
+    // delete the implicit complement tree (but not the surface trees)
+    rval = remove_bvh(implicit_complement, true);
+    MB_CHK_SET_ERR(rval, "Failed to delete the implicit complement tree");
+
+    // build the BVH for the new graveyard volume
+    rval = build_bvh(volume_set);
+    MB_CHK_SET_ERR(
+        rval,
+        "Failed to build accel. data structure for the new graveyard volume");
+    // re-build the BVH for the implicit complement
+    rval = build_bvh(implicit_complement);
+    MB_CHK_SET_ERR(rval,
+                   "Failed to build accel. data structure for the new implicit "
+                   "complement");
+  }
 
   // re-initialize indices
   rval = setup_indices();

diff --git a/src/dagmc/tests/dagmc_graveyard_test.cpp b/src/dagmc/tests/dagmc_graveyard_test.cpp
@@ -253,3 +253,160 @@ TEST_F(DagmcGraveyardTest, dagmc_graveyard_test_trelis_file) {
   EXPECT_EQ(model_sets.size(), ending_sets.size());
   EXPECT_EQ(0, subtract(model_sets, ending_sets).size());
 }
+
+TEST_F(DagmcGraveyardTest, dagmc_graveyard_test_trelis_file_no_bvh) {
+  // create DAGMC instance
+  std::unique_ptr<DagMC> DAG(new DagMC());
+  // load the trelis test file
+  ErrorCode rval = DAG->load_file(trelis_file.c_str());
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  ////////////////////////////////////////////////
+  // BVH creation would normally occur here.    //
+  // It is intentionally skipped for this test. //
+  ////////////////////////////////////////////////
+
+  rval = DAG->setup_indices();
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // collect starting sets to make sure we end up with the same thing at the end
+  Range starting_vertices;
+  rval = DAG->moab_instance()->get_entities_by_type(0, MBVERTEX,
+                                                    starting_vertices);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range starting_triangles;
+  rval =
+      DAG->moab_instance()->get_entities_by_type(0, MBTRI, starting_triangles);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range starting_sets;
+  rval =
+      DAG->moab_instance()->get_entities_by_type(0, MBENTITYSET, starting_sets);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // a graveyard is already present in this model
+  EXPECT_TRUE(DAG->has_graveyard());
+
+  int n_groups = DAG->num_entities(4);
+  int n_vols = DAG->num_entities(3);
+  int n_surfs = DAG->num_entities(2);
+  // the DagMC class only tracks relevant geometry sets (surfaces and volumes)
+  int n_curves = DAG->geom_tool()->num_ents_of_dim(1);
+  int n_geom_verts = DAG->geom_tool()->num_ents_of_dim(0);
+
+  // remove original graveyard
+  rval = DAG->remove_graveyard();
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // geometric sets removed:
+  // (geometric vertices - 16, curves - 24, surfaces - 12, volumes - 1, groups -
+  // 1)
+  EXPECT_EQ(16, n_geom_verts - DAG->geom_tool()->num_ents_of_dim(0));
+  EXPECT_EQ(24, n_curves - DAG->geom_tool()->num_ents_of_dim(1));
+  EXPECT_EQ(12, n_surfs - DAG->num_entities(2));
+  EXPECT_EQ(1, n_vols - DAG->num_entities(3));
+  EXPECT_EQ(1, n_groups - DAG->num_entities(4));
+
+  DAG->moab_instance()->write_file("test.h5m");
+  // set of vertices, triangles, and sets without the original graveyard volume
+  Range model_vertices;
+  rval =
+      DAG->moab_instance()->get_entities_by_type(0, MBVERTEX, model_vertices);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range model_triangles;
+  rval = DAG->moab_instance()->get_entities_by_type(0, MBTRI, model_triangles);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range model_sets;
+  rval = DAG->moab_instance()->get_entities_by_type(0, MBENTITYSET, model_sets);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // make sure the right number of mesh elements were removed
+  // sixteen vertices removed (corners of two cuboid volumes)
+  EXPECT_EQ(starting_vertices.size() - 16, model_vertices.size());
+  // twenty-four triangles removed (two per face for two cuboid volumes)
+  EXPECT_EQ(starting_triangles.size() - 24, model_triangles.size());
+
+  // update number of surfaces and volumes before creating the graveyard
+  n_vols = DAG->num_entities(3);
+  n_surfs = DAG->num_entities(2);
+
+  EXPECT_FALSE(DAG->has_graveyard());
+
+  // create graveyard with overwrite on when graveyard isn't present
+  rval = DAG->create_graveyard(true);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // the implicit complement will be created here too
+  // so there will be two more volumes than before
+  EXPECT_EQ(n_vols + 2, DAG->num_entities(3));
+  EXPECT_EQ(n_surfs + 2, DAG->num_entities(2));
+
+  EXPECT_TRUE(DAG->has_graveyard());
+
+  rval = DAG->remove_graveyard();
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  // here we've removed the graveyard, but the implicit complement remains
+  // so there will be one fewer volumes than the last check
+  EXPECT_EQ(n_vols + 1, DAG->num_entities(3));
+  EXPECT_EQ(n_surfs, DAG->num_entities(2));
+
+  EXPECT_FALSE(DAG->has_graveyard());
+
+  rval = DAG->create_graveyard();
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  EXPECT_TRUE(DAG->has_graveyard());
+
+  // test overwrite capability
+  bool overwrite_graveyard = true;
+  rval = DAG->create_graveyard(true);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  EXPECT_EQ(n_vols + +2, DAG->num_entities(3));
+  EXPECT_EQ(n_surfs + 2, DAG->num_entities(2));
+
+  // this should fail, graveyard already exists
+  // and overwrite isn't specified
+  rval = DAG->create_graveyard();
+  EXPECT_EQ(MB_FAILURE, rval);
+
+  rval = DAG->remove_graveyard();
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  EXPECT_FALSE(DAG->has_graveyard());
+
+  // checks to make sure we didn't accumulate any new data on the mesh
+  Range ending_vertices;
+  rval =
+      DAG->moab_instance()->get_entities_by_type(0, MBVERTEX, ending_vertices);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range ending_triangles;
+  rval = DAG->moab_instance()->get_entities_by_type(0, MBTRI, ending_triangles);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  Range ending_sets;
+  rval =
+      DAG->moab_instance()->get_entities_by_type(0, MBENTITYSET, ending_sets);
+  EXPECT_EQ(MB_SUCCESS, rval);
+
+  EXPECT_EQ(model_vertices.size(), ending_vertices.size());
+  EXPECT_EQ(0, subtract(model_vertices, ending_vertices).size());
+
+  EXPECT_EQ(model_triangles.size(), ending_triangles.size());
+  EXPECT_EQ(0, subtract(model_triangles, ending_triangles).size());
+
+  // here we've added the implict complement before setting the model sets
+  // so there should be one extra set and it will be equal to the implicit
+  // complement set
+  EXPECT_EQ(model_sets.size() + 1, ending_sets.size());
+  EXPECT_EQ(1, subtract(ending_sets, model_sets).size());
+  EntityHandle implicit_complement;
+  rval = DAG->geom_tool()->get_implicit_complement(implicit_complement);
+  EXPECT_EQ(MB_SUCCESS, rval);
+  EXPECT_EQ(implicit_complement, subtract(ending_sets, model_sets)[0]);
+}