Merge pull request #935 from pshriwise/ipc-in-da-back

Ensure implicit complement is placed at the back of the DAGMC index space

doc/CHANGELOG.rst

@@ -10,6 +10,7 @@ Next version
 **Changed:**
 
   * Update hdf5 to v1.14.3 from v1.10.4 (#931 #933)
+  * Ensure implicit complement handle is placed at the back of DAGMC volume indices (#935)
 
 v3.2.3
 ====================

src/dagmc/DagMC.cpp

@@ -894,12 +894,6 @@ ErrorCode DagMC::build_indices(Range& surfs, Range& vols) {
     logger.message("Volumes or Surfaces not found");
     return MB_ENTITY_NOT_FOUND;
   }
-  setOffset = std::min(*surfs.begin(), *vols.begin());
-  // surf/vol offsets are just first handles
-  EntityHandle tmp_offset = std::max(surfs.back(), vols.back());
-
-  // set size
-  entIndices.resize(tmp_offset - setOffset + 1);
 
   // store surf/vol handles lists (surf/vol by index) and
   // index by handle lists
@@ -911,8 +905,10 @@ ErrorCode DagMC::build_indices(Range& surfs, Range& vols) {
   *(iter++) = 0;
   std::copy(surfs.begin(), surfs.end(), iter);
   int idx = 1;
-  for (Range::iterator rit = surfs.begin(); rit != surfs.end(); ++rit)
-    entIndices[*rit - setOffset] = idx++;
+  for (auto surf_handle : surf_handles()) {
+    if (surf_handle == 0) continue;
+    entIndices[surf_handle] = idx++;
+  }
 
   vol_handles().resize(vols.size() + 1);
   iter = vol_handles().begin();
@@ -922,9 +918,32 @@ ErrorCode DagMC::build_indices(Range& surfs, Range& vols) {
   // (iter++) thereafter.
   *(iter++) = 0;
   std::copy(vols.begin(), vols.end(), iter);
+
+  // Ensure the implicit complement volume is placed at the back of this vector
+  //   Many codes iterate over the DAGMC volumes by index and all explicit
+  //   volumes should be checked before the implicit complement
+  EntityHandle implicit_complement{0};
+  rval = geom_tool()->get_implicit_complement(implicit_complement);
+  if (rval == MB_SUCCESS && implicit_complement != 0) {
+    auto it = std::find(vol_handles().begin(), vol_handles().end(),
+                        implicit_complement);
+    if (it != vol_handles().end()) {
+      vol_handles().erase(it);
+    } else {
+      logger.message(
+          "Could not find the implicit complement in the volume handles "
+          "vector");
+      return MB_FAILURE;
+    }
+    // insert the implicit complement at the end of the vector
+    vol_handles().push_back(implicit_complement);
+  }
+
   idx = 1;
-  for (Range::iterator rit = vols.begin(); rit != vols.end(); ++rit)
-    entIndices[*rit - setOffset] = idx++;
+  for (auto vol_handle : vol_handles()) {
+    if (vol_handle == 0) continue;
+    entIndices[vol_handle] = idx++;
+  }
 
   // get group handles
   Range groups;

src/dagmc/DagMC.hpp

@@ -8,6 +8,7 @@
 #include <memory>
 #include <stdexcept>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "DagMCVersion.hpp"
@@ -499,10 +500,8 @@ class DagMC {
  private:
   /** store some lists indexed by handle */
   std::vector<EntityHandle> entHandles[5];
-  /** lowest-valued handle among entity sets representing surfs and vols */
-  EntityHandle setOffset;
-  /** entity index (contiguous 1-N indices); indexed like rootSets */
-  std::vector<int> entIndices;
+  /** surface and volume mapping from EntitiyHandle to DAGMC index */
+  std::unordered_map<EntityHandle, int> entIndices;
   /** corresponding geometric entities; also indexed like rootSets */
   std::vector<RefEntity*> geomEntities;
 
@@ -571,8 +570,8 @@ inline EntityHandle DagMC::entity_by_index(int dimension, int index) const {
 }
 
 inline int DagMC::index_by_handle(EntityHandle handle) const {
-  assert(handle - setOffset < entIndices.size());
-  return entIndices[handle - setOffset];
+  assert(entIndices.count(handle) > 0);
+  return entIndices.at(handle);
 }
 
 inline unsigned int DagMC::num_entities(int dimension) const {

src/dagmc/tests/CMakeLists.txt

@@ -12,6 +12,7 @@ dagmc_install_test(dagmc_pointinvol_test cpp)
 dagmc_install_test(dagmc_rayfire_test    cpp)
 dagmc_install_test(dagmc_simple_test     cpp)
 dagmc_install_test(dagmc_graveyard_test  cpp)
+dagmc_install_test(dagmc_ipc_index_test cpp)
 
 dagmc_install_test_file(test_dagmc.h5m)
 dagmc_install_test_file(test_dagmc_impl.h5m)

src/dagmc/tests/dagmc_ipc_index_test.cpp

@@ -0,0 +1,79 @@
+#include <gtest/gtest.h>
+
+#include <filesystem>
+#include <iostream>
+#include <memory>
+
+#include "DagMC.hpp"
+#include "moab/Core.hpp"
+#include "moab/Interface.hpp"
+
+static std::string simple_file = "test_dagmc.h5m";
+
+class DagmcIPCPositionTest : public ::testing::Test {
+ protected:
+  virtual void SetUp() {}
+  virtual void TearDown() {
+    if (std::filesystem::exists("tmp.h5m")) {
+      std::filesystem::remove("tmp.h5m");
+    }
+  }
+};
+
+using namespace moab;
+
+// This test exists to ensure that the IPC is correctly positioned in the model
+// See GitHub issue #934 for more information
+// The test loads a known, working DAGMC file, adds a bunch of meshsets to it,
+// and then writes a temporary file. The temporary file is then read back in
+// and the IPC is checked to ensure it is in the correct position. The condision
+// being tested is only triggered during a file read when EntityHandle sequences
+// are generated and accessed internally and cannot be reproduced
+// using the external-facing MOAB API.
+TEST_F(DagmcIPCPositionTest, dagmc_implicit_complement_position_test) {
+  // create a DAGMC instance
+  std::unique_ptr<DagMC> dagmc = std::make_unique<DagMC>();
+
+  // load the geometry test file
+  ErrorCode rval = dagmc->load_file(simple_file.c_str());
+
+  // add one million mesh sets to this file to ensure the threshold of
+  // meshsets in a file is surpassed to trigger this condition
+  // NOTE: far fewer meshsets are needed to trigger this condition, but
+  //       one million was chosen to be conservative should the internal
+  //       allocation size in MOAB change again in the future
+  EntityHandle tmp;
+  for (int i = 0; i < 1E6; i++) {
+    rval = dagmc->moab_instance()->create_meshset(MESHSET_SET, tmp);
+    ASSERT_EQ(rval, MB_SUCCESS);
+  }
+
+  // write a temprary file
+  dagmc->moab_instance()->write_file("tmp.h5m");
+
+  // create a second DAGMC instance and read the temporary file
+  std::unique_ptr<DagMC> dagmc2 = std::make_unique<DagMC>();
+  dagmc2->load_file("tmp.h5m");
+
+  // perform necessary DAGMC setup
+  dagmc2->geom_tool()->find_geomsets();
+  dagmc2->setup_impl_compl();
+  dagmc2->setup_indices();
+
+  // get the implicit complement handle, it should exist after the explicit call
+  // for its creation above
+  EntityHandle ipc = 0;
+  rval = dagmc2->geom_tool()->get_implicit_complement(ipc);
+  ASSERT_EQ(rval, MB_SUCCESS);
+  ASSERT_NE(ipc, 0);
+
+  // there are 3 volumes in the original model and we've added the implicit
+  // complement
+  int num_vols = dagmc2->num_entities(3);
+  ASSERT_EQ(num_vols, 4);
+
+  // Make sure the IPC handle index is highest
+  // Reminder: DAGMC indieces are 1-based
+  std::cout << "IPC index: " << dagmc2->index_by_handle(ipc) << std::endl;
+  ASSERT_EQ(dagmc2->index_by_handle(ipc), 4);
+}