Halting particle sim program on failed ray fire or zero distance inte…

…rsection
pshriwise · Feb 28, 2024 · a472a3a · a472a3a
1 parent 46802d7
commit a472a3a
Show file tree

Hide file tree

Showing 5 changed files with 23 additions and 8 deletions.
diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h
@@ -40,7 +40,7 @@ class RayTracer {
   std::pair<double, MeshID> ray_fire(TreeID scene,
                                      const Position& origin,
                                      const Direction& direction,
-                                     const std::vector<MeshID>* exclude_primitives = nullptr);
+                                     std::vector<MeshID>* const exclude_primitives = nullptr);
 
   void closest(TreeID scene,
                const Position& origin,

diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h
@@ -36,7 +36,7 @@ bool point_in_volume(MeshID volume,
 std::pair<double, MeshID> ray_fire(MeshID volume,
                                    const Position& origin,
                                    const Direction& direction,
-                                   const std::vector<MeshID>* exclude_primitives = nullptr) const;
+                                   std::vector<MeshID>* const exclude_primitives = nullptr) const;
 
 void closest(MeshID volume,
               const Position& origin,

diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp
@@ -150,7 +150,7 @@ std::pair<double, MeshID>
 RayTracer::ray_fire(TreeID scene,
                     const Position& origin,
                     const Direction& direction,
-                    const std::vector<MeshID>* exclude_primitves)
+                    std::vector<MeshID>* const exclude_primitves)
 {
   RTCDRayHit rayhit;
   // set ray data
@@ -175,6 +175,7 @@ RayTracer::ray_fire(TreeID scene,
   if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID)
     return {INFTY, ID_NONE};
   else
+    if (exclude_primitves) exclude_primitves->push_back(rayhit.hit.primitive_ref->primitive_id);
     return {rayhit.ray.dtfar, rayhit.hit.surface};
 }
 

diff --git a/src/xdg.cpp b/src/xdg.cpp
@@ -59,7 +59,7 @@ std::pair<double, MeshID>
 XDG::ray_fire(MeshID volume,
               const Position& origin,
               const Direction& direction,
-              const std::vector<MeshID>* exclude_primitives) const
+              std::vector<MeshID>* const exclude_primitives) const
 {
   TreeID scene = volume_to_scene_map_.at(volume);
   return ray_tracing_interface()->ray_fire(scene, origin, direction, exclude_primitives);

diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp
@@ -39,7 +39,7 @@ void surf_dist() {
     return;
   }
   if (surface_intersection_.second == ID_NONE) {
-    std::cerr << "Particle " << id_ << " lost in volume " << volume_ << std::endl;
+    fatal_error("Particle {} lost in volume {}", id_, volume_);
     alive_ = false;
     return;
   }
@@ -52,13 +52,15 @@ void sample_collision_distance() {
 
 void collide() {
   n_events_++;
+  log("Event {} for particle {}", n_events_, id_);
   u_ = rand_dir();
   log("Particle {} collides with material at position ({}, {}, {}), new direction is ({}, {}, {})", id_, r_.x, r_.y, r_.z, u_.z, u_.y, u_.z);
   history_.clear();
 }
 
 void advance()
 {
+  log("Comparing surface intersection distance {} to collision distance {}", surface_intersection_.first, collision_distance_);
   if (collision_distance_ < surface_intersection_.first) {
     r_ += collision_distance_ * u_;
     log("Particle {} collides with material at position ({}, {}, {}) ", id_, r_.x, r_.y, r_.z);
@@ -72,14 +74,26 @@ void advance()
 void cross_surface()
 {
   n_events_++;
+  log("Event {} for particle {}", n_events_, id_);
   // check for the surface boundary condition
   if (xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION).value == "reflecting") {
     log("Particle {} reflects off surface {}", id_, surface_intersection_.second);
+    log("Direction before reflection: ({}, {}, {})", u_.x, u_.y, u_.z);
 
     Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_);
-    u_ = u_ - 2.0 * dot(u_, normal) * normal;
+    log("Normal to surface: ({}, {}, {})", normal.x, normal.y, normal.z);
+
+    double proj = dot(normal, u_);
+    double mag = normal.length();
+    normal = normal * (2.0 * proj/mag);
+    u_ = u_ - normal;
     u_ = u_.normalize();
-    if (history_.size() > 0) history_ = {history_.back()}; // reset to last intersection
+    log("Direction after reflection: ({}, {}, {})", u_.x, u_.y, u_.z);
+    // reset to last intersection
+    if (history_.size() > 0) {
+      log("Resetting particle history to last intersection");
+      history_ = {history_.back()};
+    }
   } else {
     volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second);
     log("Particle {} enters volume {}", id_, volume_);
@@ -127,7 +141,7 @@ const int n_particles {100};
 
 const int max_events {1000};
 
-bool verbose = false;
+bool verbose = true;
 
  for (int i = 0; i < n_particles; i++) {
  write_message("Starting particle {}", i);