Merge pull request #44 from pshriwise/global_scene

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];