Faster find_volume with global scene

include/double_down/RTI.hpp

@@ -248,6 +248,10 @@ class RayTracingInterface {
   moab::ErrorCode
   createBVH(moab::EntityHandle volume);
 
+  //! \brief Create a global scene containing all geometries (surfaces) in the
+  //! model
+  void create_global_scene();
+
   //! \brief Deletes the BVH for the volume if present.
   //! \param volume MOAB EntityHandle of the volume.
   void
@@ -342,6 +346,7 @@ class RayTracingInterface {
   std::unordered_map<moab::EntityHandle, RTCScene> scene_map; //!< Mapping from MOAB volume EntityHandle's to Embree Scenes.
   double numerical_precision {1E-3}; //!< Numerical precision for triangle intersections.
   double overlap_thickness {0.0}; //!< Allowed overlap thickness for self-intersecting volumes.
+  RTCScene global_scene {nullptr};
   RTCDevice g_device {nullptr}; //!< Embree device object.
 };
 

include/double_down/ray.h

@@ -9,7 +9,7 @@
 
 namespace double_down {
 
-enum RayFireType { RF, PIV, ACCUM };
+enum RayFireType { RF, PIV, ACCUM, FV };
 
 // TO-DO: there should be a few more double elements here (barycentric coords)
 

src/RTI.cpp

@@ -50,6 +50,8 @@ moab::ErrorCode RayTracingInterface::init()
     createBVH(vol);
   } // end volume loop
 
+  create_global_scene();
+
   return moab::MB_SUCCESS;
 }
 
@@ -150,8 +152,6 @@ RayTracingInterface::get_bbox(moab::EntityHandle volume,
   return rval;
 }
 
-
-
 moab::ErrorCode
 RayTracingInterface::allocateTriangleBuffer(moab::EntityHandle volume)
 {
@@ -193,7 +193,7 @@ moab::ErrorCode RayTracingInterface::createBVH(moab::EntityHandle volume)
   // add new scene to our vector
   RTCScene scene = rtcNewScene(g_device);
   rtcSetSceneFlags(scene, RTC_SCENE_FLAG_ROBUST);
-  rtcSetSceneBuildQuality(scene,RTC_BUILD_QUALITY_HIGH);
+  rtcSetSceneBuildQuality(scene, RTC_BUILD_QUALITY_HIGH);
 
   // add this scene to the map from MOAB volumes to Embree scenes
   scene_map[volume] = scene;
@@ -309,6 +309,26 @@ moab::ErrorCode RayTracingInterface::createBVH(moab::EntityHandle volume)
     return moab::MB_SUCCESS;
 }
 
+void RayTracingInterface::create_global_scene() {
+
+  // re-create the scene if one is already present
+  if (global_scene) {
+    rtcReleaseScene(global_scene);
+  }
+
+  global_scene = rtcNewScene(g_device);
+  rtcSetSceneFlags(global_scene, RTC_SCENE_FLAG_ROBUST);
+  rtcSetSceneBuildQuality(global_scene, RTC_BUILD_QUALITY_HIGH);
+
+  // loop over all geometries
+  for (auto entry : data_ptr_map) {
+    auto geom = entry.first;
+    rtcAttachGeometry(global_scene, geom);
+  }
+
+  rtcCommitScene(global_scene);
+}
+
 void RayTracingInterface::deleteBVH(moab::EntityHandle volume)
 {
 
@@ -700,9 +720,59 @@ moab::ErrorCode
 RayTracingInterface::find_volume(const double xyz[3],
                                  moab::EntityHandle& volume,
                                  const double* uvw) {
+  moab::ErrorCode rval;
+
+  const double huge_val = std::numeric_limits<double>::max();
+  double dist_limit = huge_val;
+
+  MBRayHit rayhit;
+
+  MBRay& mbray = rayhit.ray;
+  mbray.set_org(xyz);
+  if (uvw) mbray.set_dir(uvw);
+  else mbray.set_dir({0.7071, 0.7071, 0.0});
+
+  // initial ray information
+  mbray.tnear = 0.0;
+  mbray.set_len(dist_limit);
+  // use normal ray fire for exiting intersection
+  mbray.rf_type = RayFireType::FV;
+  mbray.orientation = 0;
+  mbray.mask = -1;
+
+  // initial hit information
+  MBHit& mbhit = rayhit.hit;
+  mbhit.geomID = RTC_INVALID_GEOMETRY_ID;
+  mbhit.primID = RTC_INVALID_GEOMETRY_ID;
+
+  // fire ray
+  {
+    rtcIntersect1(global_scene, (RTCRayHit*)&rayhit);
+    mbhit.Ng_x = -mbhit.Ng_x;
+    mbhit.Ng_y = -mbhit.Ng_y;
+    mbhit.Ng_z = -mbhit.Ng_z;
+  }
+
+  Vec3da ray_dir(mbray.ddir);
+  Vec3da tri_norm(rayhit.hit.dNg[0], rayhit.hit.dNg[1], rayhit.hit.dNg[2]);
+
+  if (mbhit.geomID != RTC_INVALID_GEOMETRY_ID) {
+    // get the volumes on either side of this geometry
+    moab::EntityHandle fwd, bwd;
+    rval = GTT->get_surface_senses(mbhit.surf_handle, fwd, bwd);
+    MB_CHK_SET_ERR_CONT(rval, "Failed to get parent meshsets");
+    if (dot(ray_dir, tri_norm) > 0.0) {
+      volume = fwd;
+    } else {
+      volume = bwd;
+    }
+    return MB_SUCCESS;
+  }
+
+  // if the ray misses for some reason, fall back on the linear search method
   int result = 0;
   moab::Range vols;
-  moab::ErrorCode rval = GTT->get_gsets_by_dimension(3, vols);
+  rval = GTT->get_gsets_by_dimension(3, vols);
   if (rval != moab::MB_SUCCESS) return rval;
   for (moab::EntityHandle vol : vols) {
     rval = point_in_volume(vol, xyz, result, uvw);

src/primitives.cpp

@@ -88,7 +88,7 @@ void DblTriIntersectFunc(RTCIntersectFunctionNArguments* args) {
     normal.normalize();
 
     // flip the triangle normal if the sense is reversed
-    if( -1 == this_tri.sense ) normal *= -1;
+    if (-1 == this_tri.sense && ray.rf_type != RayFireType::FV) normal *= -1;
 
     // set the double precision normal of the hit
     hit.dNg[0] = normal[0];

test/CMakeLists.txt

@@ -13,6 +13,12 @@ target_link_libraries(test_rf ${MOAB_LIBRARIES} ${EMBREE_LIBRARY} dd)
 set_target_properties(test_rf PROPERTIES BUILD_RPATH "${MOAB_LIBRARY_DIRS}")
 set_target_properties(test_rf PROPERTIES BUILD_RPATH_USE_LINK_PATH TRUE)
 
+add_executable(test_find_vol test_find_vol.cpp test_utils.cpp)
+target_include_directories(test_find_vol PUBLIC ${MOAB_INCLUDE_DIRS} ${CMAKE_SOURCE_DIR}/include/)
+target_link_libraries(test_find_vol ${MOAB_LIBRARIES} ${EMBREE_LIBRARY} dd)
+set_target_properties(test_find_vol PROPERTIES BUILD_RPATH "${MOAB_LIBRARY_DIRS}")
+set_target_properties(test_find_vol PROPERTIES BUILD_RPATH_USE_LINK_PATH TRUE)
+
 add_executable(test_build test_build.cpp test_utils.cpp)
 target_include_directories(test_build PUBLIC ${MOAB_INCLUDE_DIRS} ${CMAKE_SOURCE_DIR}/include/)
 target_link_libraries(test_build ${MOAB_LIBRARIES} ${EMBREE_LIBRARY} dd)
@@ -28,5 +34,6 @@ set_target_properties(test_closest PROPERTIES BUILD_RPATH_USE_LINK_PATH TRUE)
 
 add_test(NAME test_moab_write COMMAND test_mb)
 add_test(NAME test_rti_rf COMMAND test_rf)
+add_test(NAME test_find_volume COMMAND test_find_vol)
 add_test(NAME test_rti_build COMMAND test_build)
 add_test(NAME test_rti_closest COMMAND test_closest)

test/test_find_vol.cpp

@@ -0,0 +1,59 @@
+#include <cstdlib>
+#include <iostream>
+#include <math.h>
+
+#include "test_utils.hpp"
+
+#include "double_down/RTI.hpp"
+
+using namespace double_down;
+
+int main() {
+
+  // create new MOAB instance
+  moab::Interface* MBI = new moab::Core();
+
+  moab::ErrorCode rval;
+
+  rval = MBI->load_file("../test_files/sphere.h5m");
+  MB_CHK_SET_ERR(rval, "Failed to load test file");
+
+  std::unique_ptr<RayTracingInterface> RTI(new RayTracingInterface(MBI));
+
+  rval = RTI->init();
+  MB_CHK_SET_ERR(rval, "Failed to initialize the RTI.");
+
+  moab::Range vols;
+  rval = RTI->get_vols(vols);
+  MB_CHK_SET_ERR(rval, "Failed to get volumes from the RTI.");
+
+  moab::EntityHandle sphere_vol = vols[0];
+
+  // fire a test ray
+  double org[3] = {0.0, 0.0, 0.0};
+  double dir[3] = {1.0, 0.0, 0.0};
+
+  double dist = 0.0;
+  moab::EntityHandle surf;
+  for (int i = 0; i < 1000; i++) {
+    std::cout << org[0] << " " << org[1] << " " << org[2] << std::endl;
+    moab::EntityHandle volume = 0;
+    RTI->find_volume(org, volume, dir);
+    if (volume != sphere_vol) { return 1; }
+
+    // test this without a direction supplied
+    volume = 0;
+    RTI->find_volume(org, volume);
+    if (volume != sphere_vol) { return 1; }
+
+    // sample a random point in the sphere volume
+    double mu = -1.0 + 2.0 * (double)std::rand() / (double)RAND_MAX;
+    double phi = 2.0 * M_PI * (double)std::rand() / (double)RAND_MAX;
+    double r = 9.9 * (double)std::rand() / (double)RAND_MAX;
+    org[0] = sqrt(1-mu*mu)*cos(phi)*r;
+    org[1] = sqrt(1-mu*mu)*sin(phi)*r;
+    org[2] = r*mu;
+  }
+
+  return 0;
+}