From 758826762d7eaf4e5db67a6ffa47d5e09b89ed9d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 17 Jan 2024 00:09:33 -0600 Subject: [PATCH 01/12] Adding the start of a simpleparticle mini-app --- CMakeLists.txt | 5 +++ include/xdg/ray.h | 4 +++ include/xdg/ray_tracing_interface.h | 3 +- include/xdg/vec3da.h | 2 +- include/xdg/xdg.h | 4 +++ src/ray_tracing_interface.cpp | 2 +- tools/CMakeLists.txt | 4 +++ tools/particle_sim.cpp | 47 +++++++++++++++++++++++++++++ 8 files changed, 67 insertions(+), 4 deletions(-) create mode 100644 tools/CMakeLists.txt create mode 100644 tools/particle_sim.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e9033e6..c20ee60 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,7 @@ option(XDG_ENABLE_MOAB "Enable support for the MOAB mesh library" ON) option(XDG_ENABLE_MFEM "Enable support for the MFEM mesh library" OFF) option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) +option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) # Set version numbers set(XDG_VERSION_MAJOR 0) @@ -80,6 +81,10 @@ if (XDG_BUILD_TESTS) add_subdirectory("tests") endif() +if (XDG_BUILD_TOOLS) + add_subdirectory("tools") +endif() + target_include_directories(xdg PUBLIC $ diff --git a/include/xdg/ray.h b/include/xdg/ray.h index c2755b2..92fcb3c 100644 --- a/include/xdg/ray.h +++ b/include/xdg/ray.h @@ -3,6 +3,7 @@ #define _XDG_RAY_H #include +#include #include "xdg/vec3da.h" @@ -13,6 +14,9 @@ namespace xdg { +// forward declaration +class TriangleRef; + // TO-DO: there should be a few more double elements here (barycentric coords) /*! Stucture that is an extension of Embree's RTCRay with diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 09ca831..9b88a95 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -30,11 +30,10 @@ class RayTracer { TreeID register_volume(const std::shared_ptr mesh_manager, MeshID volume); // Query Methods - bool point_in_volume(TreeID scene, const Position& point, const Direction* direction = nullptr, - const std::vector* exclude_primitives = nullptr); + const std::vector* exclude_primitives = nullptr) const; void ray_fire(TreeID scene, const Position& origin, diff --git a/include/xdg/vec3da.h b/include/xdg/vec3da.h index 3ea3a07..1fc3212 100644 --- a/include/xdg/vec3da.h +++ b/include/xdg/vec3da.h @@ -150,7 +150,7 @@ __forceinline std::ostream& operator <<(std::ostream &os, Vec3da const& v) { inline bool lower(const Vec3da& a, const Vec3da& b) { for (int i = 0; i < 3; i++) - if (a[i] != b[i]) + if (a[i] != b[i]) return a[i] < b[i]; return false; } diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 0a7da84..6cddba4 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -57,6 +57,10 @@ Direction surface_normal(MeshID surface, Position point, const std::vector* exclude_primitives = nullptr) const; + MeshID find_volume(const Position& point, + const Direction& direction) const; + + // Geometric Measurements double measure_volume(MeshID volume) const; double measure_surface_area(MeshID surface) const; diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 11cc08e..02c73eb 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -118,7 +118,7 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, bool RayTracer::point_in_volume(TreeID scene, const Position& point, const Direction* direction, - const std::vector* exclude_primitives) + const std::vector* exclude_primitives) const { RTCDRayHit rayhit; rayhit.ray.set_org(point); diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt new file mode 100644 index 0000000..ecc6051 --- /dev/null +++ b/tools/CMakeLists.txt @@ -0,0 +1,4 @@ + + +add_executable(particle-sim particle_sim.cpp) +target_link_libraries(particle-sim xdg) \ No newline at end of file diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp new file mode 100644 index 0000000..d5ebcf6 --- /dev/null +++ b/tools/particle_sim.cpp @@ -0,0 +1,47 @@ +#include +#include +#include +#include + +#include "xdg/mesh_manager_interface.h" +#include "xdg/moab/mesh_manager.h" +#include "xdg/triangle_ref.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +using namespace xdg; + +struct Particle { + Position r; + Direction u; + MeshID volume; + std::vector history; +}; + +int main(int argc, char** argv) { + +// create a mesh manager +std::shared_ptr mm = std::make_shared(); + +std::shared_ptr xdg = + std::make_shared(); +xdg->set_mesh_manager_interface(mm); + +std::string filename {argv[1]}; + +mm->load_file(filename); +mm->init(); + +xdg->prepare_raytracer(); + +// create a new particle +Particle p{ {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {ID_NONE}, {}}; + +// determine what volume this particle is in +p.volume = xdg->find_volume(p.r, p.u); + +std::cout << "Particle starts in volume " << p.volume << std::endl; + +return 0; + +} \ No newline at end of file From 80724a9045fc9104b831a0018c58dd3c7594e23e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 17 Jan 2024 00:13:03 -0600 Subject: [PATCH 02/12] Update class name in sim tool --- tools/particle_sim.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index d5ebcf6..5acc9cb 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -23,8 +23,7 @@ int main(int argc, char** argv) { // create a mesh manager std::shared_ptr mm = std::make_shared(); -std::shared_ptr xdg = - std::make_shared(); +std::shared_ptr xdg = std::make_shared(); xdg->set_mesh_manager_interface(mm); std::string filename {argv[1]}; From 7fe37a19efb142de8a730a4df712748aed537038 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jan 2024 11:31:07 -0600 Subject: [PATCH 03/12] Weridness --- include/xdg/ray.h | 1 + include/xdg/ray_tracing_interface.h | 9 ++++---- include/xdg/xdg.h | 27 +++++++++++------------ src/ray_tracing_interface.cpp | 7 +++--- src/triangle_intersect.cpp | 1 + src/xdg.cpp | 15 +++---------- tests/test_moab.cpp | 20 +++++++++++------ tests/test_ray_fire.cpp | 34 ++++++++++++++--------------- tools/particle_sim.cpp | 24 +++++++++++++------- 9 files changed, 71 insertions(+), 67 deletions(-) diff --git a/include/xdg/ray.h b/include/xdg/ray.h index 92fcb3c..a30c034 100644 --- a/include/xdg/ray.h +++ b/include/xdg/ray.h @@ -94,6 +94,7 @@ struct RTCDHit : RTCHit { // data members const PrimitiveRef* primitive_ref {nullptr}; //!< Pointer to the primitive reference for this hit + MeshID surface; //!< ID of the surface this hit belongs to Vec3da dNg; //!< Double precision version of the primitive normal }; diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 9b88a95..e95a631 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -35,11 +35,10 @@ class RayTracer { const Direction* direction = nullptr, const std::vector* exclude_primitives = nullptr) const; - void ray_fire(TreeID scene, - const Position& origin, - const Direction& direction, - double& distance, - const std::vector* exclude_primitives = nullptr); + std::pairray_fire(TreeID scene, + const Position& origin, + const Direction& direction, + const std::vector* exclude_primitives = nullptr); void closest(TreeID scene, const Position& origin, diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 6cddba4..2ce9f22 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -28,16 +28,19 @@ class XDG { MeshID find_volume(const Position& point, const Direction& direction) const; -bool point_in_volume(MeshID volume, - const Position& point, - const Direction* direction = nullptr, - const std::vector* exclude_primitives = nullptr) const; +inline bool point_in_volume(MeshID volume, + const Position point, + const Direction* direction = nullptr, + const std::vector* exclude_primitives = nullptr) const +{ + TreeID scene = volume_to_scene_map_.at(volume); + return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); +} -void ray_fire(MeshID volume, - const Position& origin, - const Direction& direction, - double& distance, - const std::vector* exclude_primitives = nullptr) const; +std::pair ray_fire(MeshID volume, + const Position& origin, + const Direction& direction, + const std::vector* exclude_primitives = nullptr) const; void closest(MeshID volume, const Position& origin, @@ -57,9 +60,6 @@ Direction surface_normal(MeshID surface, Position point, const std::vector* exclude_primitives = nullptr) const; - MeshID find_volume(const Position& point, - const Direction& direction) const; - // Geometric Measurements double measure_volume(MeshID volume) const; @@ -79,8 +79,8 @@ Direction surface_normal(MeshID surface, const std::shared_ptr& mesh_manager() const { return mesh_manager_; } - // Private methods + std::map volume_to_scene_map_; // ray_tracing_interface_ {std::make_shared()}; std::shared_ptr mesh_manager_ {nullptr}; - std::map volume_to_scene_map_; // surface_to_scene_map_; // surface_to_geometry_map_; // 0.0; } -void +std::pair RayTracer::ray_fire(TreeID scene, const Position& origin, const Direction& direction, - double& distance, const std::vector* exclude_primitves) { RTCDRayHit rayhit; @@ -170,9 +169,9 @@ RayTracer::ray_fire(TreeID scene, } if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID) - distance = INFTY; + return {INFTY, ID_NONE}; else - distance = rayhit.ray.dtfar; + return {rayhit.ray.dtfar, rayhit.hit.surface}; } void RayTracer::closest(TreeID scene, diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index 2174cfd..dfb8e11 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -91,6 +91,7 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { rayhit->hit.geomID = args->geomID; rayhit->hit.primID = args->primID; rayhit->hit.primitive_ref = &primitive_ref; + rayhit->hit.surface = user_data->surface_id; rayhit->hit.dNg = normal; } diff --git a/src/xdg.cpp b/src/xdg.cpp index 840d3d9..8ab257a 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -46,23 +46,14 @@ MeshID XDG::find_volume(const Position& point, return ID_NONE; } -bool XDG::point_in_volume(MeshID volume, - const Position& point, - const Direction* direction, - const std::vector* exclude_primitives) const -{ - TreeID scene = volume_to_scene_map_.at(volume); - return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); -} - -void XDG::ray_fire(MeshID volume, +std::pair +XDG::ray_fire(MeshID volume, const Position& origin, const Direction& direction, - double& distance, const std::vector* exclude_primitives) const { TreeID scene = volume_to_scene_map_.at(volume); - ray_tracing_interface()->ray_fire(scene, origin, direction, distance, exclude_primitives); + return ray_tracing_interface()->ray_fire(scene, origin, direction, exclude_primitives); } void XDG::closest(MeshID volume, diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index 6ba9eea..83d81d8 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -84,20 +84,26 @@ TEST_CASE("Test Ray Fire MOAB") Position origin {0.0, 0.0, 0.0}; Direction direction {1.0, 0.0, 0.0}; - double intersection_distance {0.0}; + std::pair intersection; - xdg->ray_fire(volume, origin, direction, intersection_distance); + intersection = xdg->ray_fire(volume, origin, direction); // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(5.0, 1e-6)); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); origin = {3.0, 0.0, 0.0}; - xdg->ray_fire(volume, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(2.0, 1e-6)); + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); origin = {-10.0, 0.0, 0.0}; - xdg->ray_fire(volume, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(15.0, 1e-6)); + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); + + origin = {0.0, 0.0, 0.0}; + TreeID scene = xdg->volume_to_scene_map_[volume]; + REQUIRE(xdg->ray_tracing_interface()->point_in_volume(scene, origin)); + REQUIRE(xdg->point_in_volume(volume, origin)); + } TEST_CASE("TEST XDG Factory Method") diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index 1bd6914..1adec46 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -23,41 +23,41 @@ TEST_CASE("Test Ray Fire Mesh Mock") Position origin {0.0, 0.0, 0.0}; Direction direction {1.0, 0.0, 0.0}; - double intersection_distance {0.0}; + std::pair intersection; // fire from the origin toward each face, ensuring that the intersection distances are correct - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(5.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(2.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); direction = {0.0, 1.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(6.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(6.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(3.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(3.0, 1e-6)); direction = {0.0, 0.0, 1.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(7.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(7.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(4.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(4.0, 1e-6)); // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct // rays should skip entering intersections and intersect with the far side of the cube origin = {-10.0, 0.0, 0.0}; direction = {1.0, 0.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(15.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); origin = {10.0, 0.0, 0.0}; direction = {-1.0, 0.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(12.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(12.0, 1e-6)); } \ No newline at end of file diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index 5acc9cb..34508c9 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -5,7 +5,6 @@ #include "xdg/mesh_manager_interface.h" #include "xdg/moab/mesh_manager.h" -#include "xdg/triangle_ref.h" #include "xdg/vec3da.h" #include "xdg/xdg.h" @@ -15,22 +14,19 @@ struct Particle { Position r; Direction u; MeshID volume; - std::vector history; + std::vector history; }; int main(int argc, char** argv) { // create a mesh manager -std::shared_ptr mm = std::make_shared(); - -std::shared_ptr xdg = std::make_shared(); -xdg->set_mesh_manager_interface(mm); +std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); +const auto& mm = xdg->mesh_manager(); std::string filename {argv[1]}; mm->load_file(filename); mm->init(); - xdg->prepare_raytracer(); // create a new particle @@ -38,9 +34,21 @@ Particle p{ {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {ID_NONE}, {}}; // determine what volume this particle is in p.volume = xdg->find_volume(p.r, p.u); - +p.volume = 1; std::cout << "Particle starts in volume " << p.volume << std::endl; +std::pair intersection; +// TreeID scene = xdg->volume_to_scene_map_[p.volume]; +// intersection = xdg->ray_tracing_interface()->ray_fire(scene, p.r, p.u, &p.history); +intersection = xdg->ray_fire(p.volume, p.r, p.u, &p.history); + +std::cout << "Intersected volume " << intersection.second << " at distance " << intersection.first << std::endl; + +p.volume = xdg->mesh_manager()->next_volume(p.volume, intersection.second); +p.r += intersection.first * p.u; + +std::cout << "Particle enters volume " << p.volume << " at position " << p.r << std::endl; + return 0; } \ No newline at end of file From ba2da255f5751bcd38466c81f6576a88ed9f72e3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jan 2024 11:51:27 -0600 Subject: [PATCH 04/12] Fixes? --- include/xdg/mesh_manager_interface.h | 2 +- include/xdg/xdg.h | 13 ++++--------- src/mesh_manager_interface.cpp | 2 +- src/xdg.cpp | 9 +++++++++ tests/test_moab.cpp | 2 -- tools/particle_sim.cpp | 7 +++---- 6 files changed, 18 insertions(+), 17 deletions(-) diff --git a/include/xdg/mesh_manager_interface.h b/include/xdg/mesh_manager_interface.h index 1d83cf0..01524e7 100644 --- a/include/xdg/mesh_manager_interface.h +++ b/include/xdg/mesh_manager_interface.h @@ -62,7 +62,7 @@ class MeshManager { virtual void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) = 0; - MeshID next_volume(MeshID surface, MeshID current_volume) const; + MeshID next_volume(MeshID current_volume, MeshID surface) const; // Methods MeshID next_volume_id() const; diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 2ce9f22..1232884 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -28,14 +28,10 @@ class XDG { MeshID find_volume(const Position& point, const Direction& direction) const; -inline bool point_in_volume(MeshID volume, +bool point_in_volume(MeshID volume, const Position point, const Direction* direction = nullptr, - const std::vector* exclude_primitives = nullptr) const -{ - TreeID scene = volume_to_scene_map_.at(volume); - return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); -} + const std::vector* exclude_primitives = nullptr) const; std::pair ray_fire(MeshID volume, const Position& origin, @@ -80,19 +76,18 @@ Direction surface_normal(MeshID surface, return mesh_manager_; } // Private methods - std::map volume_to_scene_map_; // ray_tracing_interface_ {std::make_shared()}; std::shared_ptr mesh_manager_ {nullptr}; + std::map volume_to_scene_map_; // surface_to_scene_map_; // surface_to_geometry_map_; //get_parent_volumes(surface); diff --git a/src/xdg.cpp b/src/xdg.cpp index 8ab257a..3eefe9c 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -33,6 +33,15 @@ std::shared_ptr XDG::create(MeshLibrary library) return xdg; } +bool XDG::point_in_volume(MeshID volume, + const Position point, + const Direction* direction, + const std::vector* exclude_primitives) const +{ + TreeID scene = volume_to_scene_map_.at(volume); + return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); +} + MeshID XDG::find_volume(const Position& point, const Direction& direction) const { diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index 83d81d8..236b41d 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -100,8 +100,6 @@ TEST_CASE("Test Ray Fire MOAB") REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); origin = {0.0, 0.0, 0.0}; - TreeID scene = xdg->volume_to_scene_map_[volume]; - REQUIRE(xdg->ray_tracing_interface()->point_in_volume(scene, origin)); REQUIRE(xdg->point_in_volume(volume, origin)); } diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index 34508c9..f00e86b 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -34,15 +34,14 @@ Particle p{ {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {ID_NONE}, {}}; // determine what volume this particle is in p.volume = xdg->find_volume(p.r, p.u); -p.volume = 1; + std::cout << "Particle starts in volume " << p.volume << std::endl; std::pair intersection; -// TreeID scene = xdg->volume_to_scene_map_[p.volume]; -// intersection = xdg->ray_tracing_interface()->ray_fire(scene, p.r, p.u, &p.history); + intersection = xdg->ray_fire(p.volume, p.r, p.u, &p.history); -std::cout << "Intersected volume " << intersection.second << " at distance " << intersection.first << std::endl; +std::cout << "Intersected surface " << intersection.second << " at distance " << intersection.first << std::endl; p.volume = xdg->mesh_manager()->next_volume(p.volume, intersection.second); p.r += intersection.first * p.u; From de5d177069dfd98a5511ada6f66e4e1ac59a9f5b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 18 Jan 2024 12:49:43 -0600 Subject: [PATCH 05/12] Beefing up the particle class a bit --- tools/particle_sim.cpp | 100 ++++++++++++++++++++++++++++++++--------- 1 file changed, 78 insertions(+), 22 deletions(-) diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index f00e86b..11bd76c 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -11,10 +11,68 @@ using namespace xdg; struct Particle { - Position r; - Direction u; - MeshID volume; - std::vector history; + + Particle(std::shared_ptr xdg) : xdg_(xdg) {} + + void initialize() { + // TODO: replace with sampling + r_ = {0.0, 0.0, 0.0}; + u_ = {1.0, 0.0, 0.0}; + + volume_ = xdg_->find_volume(r_, u_); + } + + void surf_dist() { + surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); + if (surface_intersection_.second == ID_NONE) { + std::cerr << "Particle lost in volume " << volume_ << std::endl; + alive_ = false; + return; + } + std::cout << "Intersected surface " << surface_intersection_.second << " at distance " << surface_intersection_.first << std::endl; + } + + void sample_collision_distance() { collision_distance_ = INFTY; } + + void collide() { + n_events_++; + std::cout << "Particle collides with material at position " << r_ << std::endl; + } + + void advance() { + if (collision_distance_ < surface_intersection_.first) { + r_ += collision_distance_ * u_; + std::cout << "Particle collides with material at position " << r_ << std::endl; + } else { + std::cout << "Particle advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; + r_ += surface_intersection_.first * u_; + } + } + + void cross_surface() { + n_events_++; + + volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); + std::cout << "Particle enters volume " << volume_ << std::endl; + if (volume_ == ID_NONE) { + alive_ = false; + return; + } + } + + // Data Members + std::shared_ptr xdg_; + + Position r_; + Direction u_; + MeshID volume_ {ID_NONE}; + std::vector history_{}; + + std::pair surface_intersection_ {INFTY, ID_NONE}; + double collision_distance_ {INFTY}; + + int32_t n_events_ {0}; + bool alive_ {true}; }; int main(int argc, char** argv) { @@ -30,24 +88,22 @@ mm->init(); xdg->prepare_raytracer(); // create a new particle -Particle p{ {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {ID_NONE}, {}}; - -// determine what volume this particle is in -p.volume = xdg->find_volume(p.r, p.u); - -std::cout << "Particle starts in volume " << p.volume << std::endl; - -std::pair intersection; - -intersection = xdg->ray_fire(p.volume, p.r, p.u, &p.history); - -std::cout << "Intersected surface " << intersection.second << " at distance " << intersection.first << std::endl; - -p.volume = xdg->mesh_manager()->next_volume(p.volume, intersection.second); -p.r += intersection.first * p.u; - -std::cout << "Particle enters volume " << p.volume << " at position " << p.r << std::endl; +Particle p(xdg); + +const int max_events {10}; +p.initialize(); +while (p.n_events_ < max_events) { + p.surf_dist(); + // terminate for leakage + if (!p.alive_) break; + p.sample_collision_distance(); + p.advance(); + if (p.surface_intersection_.first < p.collision_distance_) + p.cross_surface(); + else + p.collide(); + if (!p.alive_) break; +} return 0; - } \ No newline at end of file From e9803eecf6ab2c4547868f1aa2ab2e7bfb6ba5a5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 30 Jan 2024 15:12:29 -0600 Subject: [PATCH 06/12] Particle tool current --- CMakeLists.txt | 6 +- include/xdg/constants.h | 10 +++ include/xdg/mesh_manager_interface.h | 7 +- include/xdg/moab/mesh_manager.h | 5 -- include/xdg/ray.h | 9 ++- include/xdg/ray_tracing_interface.h | 2 + include/xdg/vec3da.h | 9 +++ src/mesh_manager_interface.cpp | 26 +++++++ src/moab/mesh_manager.cpp | 12 --- src/ray_tracing_interface.cpp | 1 + src/triangle_intersect.cpp | 4 +- src/xdg.cpp | 2 +- tests/mesh_mock.h | 10 --- tests/test_point_in_volume.cpp | 2 +- tools/particle_sim.cpp | 108 ++++++++++++++++----------- 15 files changed, 134 insertions(+), 79 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c20ee60..9817b8f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ find_package(MOAB REQUIRED HINTS ${MOAB_DIR}) endif() # use Embree for CPU ray tracing -find_package(embree 3.6.1 REQUIRED) +find_package(embree 4.0.0 REQUIRED) if (NOT ${EMBREE_VERSION} VERSION_GREATER 3.6.0) message(FATAL_ERROR "XDG requires Embree v3.6.1 or higher.") endif() @@ -77,6 +77,10 @@ else() target_compile_definitions(xdg PUBLIC XDG_EMBREE3) endif() +if (${CMAKE_BUILD_TYPE} MATCHES "Debug") + target_compile_definitions(xdg PUBLIC XDG_DEBUG) +endif() + if (XDG_BUILD_TESTS) add_subdirectory("tests") endif() diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 2de9986..0f59d15 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -9,6 +9,16 @@ namespace xdg { constexpr double INFTY {std::numeric_limits::max()}; +#ifdef XDG_DEBUG + // from Embree, if the floating point value of Tfar is larger than this value, + // it is considered overflow by the internal hit verification function. This + // is only enabled when Embree is compiled in debug mode. I think a + // corresponding behavior makes sense to keep here for now + constexpr double INFTYF {1.844E18f}; +#else + constexpr double INFTYF {std::numeric_limits::max()}; +#endif + // Whether information pertains to a surface or volume enum class GeometryType { SURFACE = 2, diff --git a/include/xdg/mesh_manager_interface.h b/include/xdg/mesh_manager_interface.h index 01524e7..fe14552 100644 --- a/include/xdg/mesh_manager_interface.h +++ b/include/xdg/mesh_manager_interface.h @@ -74,8 +74,11 @@ class MeshManager { // Metadata methods virtual void parse_metadata() = 0; - virtual Property get_volume_property(MeshID volume, PropertyType type) const = 0; - virtual Property get_surface_property(MeshID surface, PropertyType type) const = 0; + bool volume_has_property(MeshID volume, PropertyType type) const; + bool surface_has_property(MeshID surface, PropertyType type) const; + + Property get_volume_property(MeshID volume, PropertyType type) const; + Property get_surface_property(MeshID surface, PropertyType type) const; // Accessors const std::vector& volumes() const { return volumes_; } diff --git a/include/xdg/moab/mesh_manager.h b/include/xdg/moab/mesh_manager.h index 5de4c0d..dba048f 100644 --- a/include/xdg/moab/mesh_manager.h +++ b/include/xdg/moab/mesh_manager.h @@ -81,11 +81,6 @@ class MOABMeshManager : public MeshManager { // Metadata void parse_metadata() override; - // TODO: move to mesh_manger_interface??? - Property get_volume_property(MeshID volume, PropertyType type) const override; - - Property get_surface_property(MeshID surface, PropertyType type) const override; - // Other MeshLibrary mesh_library() const override {return MeshLibrary::MOAB; } diff --git a/include/xdg/ray.h b/include/xdg/ray.h index a30c034..f4d213c 100644 --- a/include/xdg/ray.h +++ b/include/xdg/ray.h @@ -27,7 +27,7 @@ struct RTCDRay: RTCRay { RTCDRay() { this->tnear = 0.0; - this->tfar = INFTY; + this->tfar = INFTYF; } //! \brief Set both the single and double precision versions of the ray origin @@ -68,7 +68,7 @@ struct RTCDRay: RTCRay { //! \brief Set both the single and double precision versions of the ray max distance void set_tfar(double d) { - tfar = d; + tfar = std::min(d, INFTYF); dtfar = d; } @@ -90,6 +90,9 @@ struct RTCDHit : RTCHit { RTCDHit() { this->geomID = RTC_INVALID_GEOMETRY_ID; this->primID = RTC_INVALID_GEOMETRY_ID; + this->Ng_x = 0.0; + this->Ng_y = 0.0; + this->Ng_z = 0.0; } // data members @@ -120,7 +123,7 @@ struct RTCDPointQuery : RTCPointQuery { //! \brief Set both the single and double precision versions of the query radius void set_radius(double rad) { - radius = rad; + radius = std::min(rad, INFTYF); dradius = rad; } diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index e95a631..b8e70b9 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -11,6 +11,8 @@ #include "xdg/primitive_ref.h" #include "xdg/geometry_data.h" + + namespace xdg { diff --git a/include/xdg/vec3da.h b/include/xdg/vec3da.h index 1fc3212..c082c2e 100644 --- a/include/xdg/vec3da.h +++ b/include/xdg/vec3da.h @@ -160,6 +160,15 @@ using Vertex = Vec3da; using Position = Vec3da; using Direction = Vec3da; +inline Direction rand_dir() { + double theta = drand48() * 2.0 * M_PI; + double u = (1.0 - drand48()) - 1.0; + double phi = acos(u); + return Direction(sin(phi) * cos(theta), sin(phi) * sin(theta), cos(phi)).normalize(); + +} + + } // end namespace xdg #endif diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index ed65f06..cc5e77e 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -40,6 +40,32 @@ MeshID MeshManager::next_surface_id() const return *std::max_element(surfaces().begin(), surfaces().end()) + 1; } +bool +MeshManager::volume_has_property(MeshID volume, PropertyType type) const +{ + return volume_metadata_.count({volume, type}) > 0; +} + +bool +MeshManager::surface_has_property(MeshID surface, PropertyType type) const +{ + return surface_metadata_.count({surface, type}) > 0; +} + +Property +MeshManager::get_volume_property(MeshID volume, PropertyType type) const +{ + return volume_metadata_.at({volume, type}); +} + +Property +MeshManager::get_surface_property(MeshID surface, PropertyType type) const +{ + if (surface_metadata_.count({surface, type}) == 0) + return {PropertyType::BOUNDARY_CONDITION, "transmission"}; + return surface_metadata_.at({surface, type}); +} + MeshID MeshManager::next_volume(MeshID current_volume, MeshID surface) const { auto parent_vols = this->get_parent_volumes(surface); diff --git a/src/moab/mesh_manager.cpp b/src/moab/mesh_manager.cpp index bb76d49..19be979 100644 --- a/src/moab/mesh_manager.cpp +++ b/src/moab/mesh_manager.cpp @@ -299,18 +299,6 @@ MOABMeshManager::get_volume_surfaces(MeshID volume) const return surface_ids; } -Property -MOABMeshManager::get_volume_property(MeshID volume, PropertyType type) const -{ - return volume_metadata_.at({volume, type}); -} - -Property -MOABMeshManager::get_surface_property(MeshID surface, PropertyType type) const -{ - return surface_metadata_.at({surface, type}); -} - void MOABMeshManager::parse_metadata() { diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 0c28427..7031f5f 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -79,6 +79,7 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, double max_distance = std::sqrt(dx*dx + dy*dy + dz*dz); double bump = max_distance * std::pow(10, -std::numeric_limits::digits10); + bump = std::max(bump, 1e-03); // create a new geometry for each surface int buffer_start = 0; diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index dfb8e11..66e4d9d 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -87,12 +87,14 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { // zero-out barycentric coords rayhit->hit.u = 0.0; rayhit->hit.v = 0.0; + rayhit->hit.Ng_x = 0.0; + rayhit->hit.Ng_y = 0.0; + rayhit->hit.Ng_z = 0.0; // set the hit information rayhit->hit.geomID = args->geomID; rayhit->hit.primID = args->primID; rayhit->hit.primitive_ref = &primitive_ref; rayhit->hit.surface = user_data->surface_id; - rayhit->hit.dNg = normal; } diff --git a/src/xdg.cpp b/src/xdg.cpp index 3eefe9c..f0f0c3b 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -96,7 +96,7 @@ Direction XDG::surface_normal(MeshID surface, const std::vector* exclude_primitives) const { MeshID element; - if (exclude_primitives != nullptr) { + if (exclude_primitives != nullptr && exclude_primitives->size() > 0) { element = exclude_primitives->back(); } else { auto surface_vols = mesh_manager()->get_parent_volumes(surface); diff --git a/tests/mesh_mock.h b/tests/mesh_mock.h index 927c99e..14fc9b1 100644 --- a/tests/mesh_mock.h +++ b/tests/mesh_mock.h @@ -120,16 +120,6 @@ class MeshMock : public MeshManager { fatal_error("MockMesh does not support parse_metadata()"); } - virtual Property get_volume_property(MeshID volume, PropertyType type) const override { - fatal_error("MockMesh does not support get_volume_property()"); - return Property(); - } - - virtual Property get_surface_property(MeshID surface, PropertyType type) const override { - fatal_error("MockMesh does not support get_surface_property()"); - return Property(); - } - // Other virtual MeshLibrary mesh_library() const override { return MeshLibrary::INTERNAL; } diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index 0c00a97..c674b8c 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -27,7 +27,7 @@ TEST_CASE("Test Point in Volume") REQUIRE(result == false); // test a point just inside the positive x boundary - point = {5.0 - 1e-06, 0.0, 0.0}; + point = {4.0 - 1e-06, 0.0, 0.0}; result = rti->point_in_volume(volume_tree, point); REQUIRE(result == true); diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index 11bd76c..f2da608 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -10,48 +10,63 @@ using namespace xdg; +static double mfp {1.0}; + struct Particle { - Particle(std::shared_ptr xdg) : xdg_(xdg) {} +Particle(std::shared_ptr xdg) : xdg_(xdg) {} - void initialize() { - // TODO: replace with sampling - r_ = {0.0, 0.0, 0.0}; - u_ = {1.0, 0.0, 0.0}; +void initialize() { + // TODO: replace with sampling + r_ = {0.0, 0.0, 0.0}; + u_ = {1.0, 0.0, 0.0}; - volume_ = xdg_->find_volume(r_, u_); - } + volume_ = xdg_->find_volume(r_, u_); +} - void surf_dist() { - surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); - if (surface_intersection_.second == ID_NONE) { - std::cerr << "Particle lost in volume " << volume_ << std::endl; - alive_ = false; - return; - } - std::cout << "Intersected surface " << surface_intersection_.second << " at distance " << surface_intersection_.first << std::endl; +void surf_dist() { + surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); + if (surface_intersection_.second == ID_NONE) { + std::cerr << "Particle lost in volume " << volume_ << std::endl; + alive_ = false; + return; } + std::cout << "Intersected surface " << surface_intersection_.second << " at distance " << surface_intersection_.first << std::endl; +} - void sample_collision_distance() { collision_distance_ = INFTY; } +void sample_collision_distance() { + collision_distance_ = -std::log(1.0 - drand48()) / mfp; +} - void collide() { - n_events_++; - std::cout << "Particle collides with material at position " << r_ << std::endl; - } +void collide() { + n_events_++; + u_ = rand_dir(); + std::cout << "Particle collides with material at position " << r_ << ", new direction is " << u_ << std::endl; + history_.clear(); +} - void advance() { - if (collision_distance_ < surface_intersection_.first) { - r_ += collision_distance_ * u_; - std::cout << "Particle collides with material at position " << r_ << std::endl; - } else { - std::cout << "Particle advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; - r_ += surface_intersection_.first * u_; - } +void advance() +{ + if (collision_distance_ < surface_intersection_.first) { + r_ += collision_distance_ * u_; + std::cout << "Particle collides with material at position " << r_ << std::endl; + } else { + r_ += surface_intersection_.first * u_; + std::cout << "Particle advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; } +} - void cross_surface() { - n_events_++; - +void cross_surface() +{ + n_events_++; + // check for the surface boundary condition + if (xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION).value == "reflecting") { + std::cout << "Particle reflects off surface " << surface_intersection_.second << std::endl; + Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); + u_ = u_ - 2.0 * dot(u_, normal) * normal; + u_ = u_.normalize(); + if (history_.size() > 0) history_ = {history_.back()}; // reset to last intersection + } else { volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); std::cout << "Particle enters volume " << volume_ << std::endl; if (volume_ == ID_NONE) { @@ -59,20 +74,21 @@ struct Particle { return; } } +} - // Data Members - std::shared_ptr xdg_; +// Data Members +std::shared_ptr xdg_; - Position r_; - Direction u_; - MeshID volume_ {ID_NONE}; - std::vector history_{}; +Position r_; +Direction u_; +MeshID volume_ {ID_NONE}; +std::vector history_{}; - std::pair surface_intersection_ {INFTY, ID_NONE}; - double collision_distance_ {INFTY}; +std::pair surface_intersection_ {INFTY, ID_NONE}; +double collision_distance_ {INFTY}; - int32_t n_events_ {0}; - bool alive_ {true}; +int32_t n_events_ {0}; +bool alive_ {true}; }; int main(int argc, char** argv) { @@ -85,14 +101,15 @@ std::string filename {argv[1]}; mm->load_file(filename); mm->init(); +mm->parse_metadata(); xdg->prepare_raytracer(); // create a new particle Particle p(xdg); -const int max_events {10}; +const int max_events {100}; p.initialize(); -while (p.n_events_ < max_events) { +while (true) { p.surf_dist(); // terminate for leakage if (!p.alive_) break; @@ -103,6 +120,11 @@ while (p.n_events_ < max_events) { else p.collide(); if (!p.alive_) break; + + if (p.n_events_ > max_events) { + std::cout << "Maximum number of events (" << max_events << ") reached" << std::endl; + break; + } } return 0; From f75ce77c23742fb807e03b6176a2890e854c6e9b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 15 Feb 2024 17:24:39 -0600 Subject: [PATCH 07/12] Initializing ray mask as empty --- include/xdg/ray.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/xdg/ray.h b/include/xdg/ray.h index f4d213c..4d0dee1 100644 --- a/include/xdg/ray.h +++ b/include/xdg/ray.h @@ -28,6 +28,7 @@ struct RTCDRay: RTCRay { RTCDRay() { this->tnear = 0.0; this->tfar = INFTYF; + this->mask = -1; } //! \brief Set both the single and double precision versions of the ray origin @@ -97,7 +98,7 @@ struct RTCDHit : RTCHit { // data members const PrimitiveRef* primitive_ref {nullptr}; //!< Pointer to the primitive reference for this hit - MeshID surface; //!< ID of the surface this hit belongs to + MeshID surface {ID_NONE}; //!< ID of the surface this hit belongs to Vec3da dNg; //!< Double precision version of the primitive normal }; From 04ff05f4a947387f5f7810c78daf042f311c7303 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 15 Feb 2024 17:25:06 -0600 Subject: [PATCH 08/12] Storing shared pointers in case memory locations are reallocated --- include/xdg/ray_tracing_interface.h | 12 ++++++------ src/ray_tracing_interface.cpp | 15 +++++++++------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index b8e70b9..6a01b48 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -37,10 +37,10 @@ class RayTracer { const Direction* direction = nullptr, const std::vector* exclude_primitives = nullptr) const; - std::pairray_fire(TreeID scene, - const Position& origin, - const Direction& direction, - const std::vector* exclude_primitives = nullptr); + std::pair ray_fire(TreeID scene, + const Position& origin, + const Direction& direction, + const std::vector* exclude_primitives = nullptr); void closest(TreeID scene, const Position& origin, @@ -59,7 +59,7 @@ class RayTracer { // Accessors int num_registered_scenes() const { return scenes_.size(); } - const GeometryUserData& geometry_data(MeshID surface) const { return user_data_map_.at(surface_to_geometry_map_.at(surface)); } + const std::shared_ptr& geometry_data(MeshID surface) const { return user_data_map_.at(surface_to_geometry_map_.at(surface)); } // Data members private: @@ -74,7 +74,7 @@ class RayTracer { RTCScene gloabal_scene_; // Internal Embree Mappings - std::unordered_map user_data_map_; + std::unordered_map> user_data_map_; // Internal parameters double numerical_precision_ {1e-3}; diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 7031f5f..e745c23 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -80,6 +80,7 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, double max_distance = std::sqrt(dx*dx + dy*dy + dz*dz); double bump = max_distance * std::pow(10, -std::numeric_limits::digits10); bump = std::max(bump, 1e-03); + bump = 1.0; // create a new geometry for each surface int buffer_start = 0; @@ -90,16 +91,17 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, unsigned int embree_surface = rtcAttachGeometry(volume_scene, surface_geometry); this->surface_to_geometry_map_[surface] = surface_geometry; - GeometryUserData surface_data; - surface_data.surface_id = surface; - surface_data.mesh_manager = mesh_manager.get(); - surface_data.prim_ref_buffer = tri_ref_ptr + buffer_start; + std::shared_ptr surface_data = std::make_shared(); + surface_data->surface_id = surface; + surface_data->mesh_manager = mesh_manager.get(); + surface_data->prim_ref_buffer = tri_ref_ptr + buffer_start; user_data_map_[surface_geometry] = surface_data; - rtcSetGeometryUserData(surface_geometry, &user_data_map_[surface_geometry]); + // TODO: This could be a problem if user_data_map_ is reallocated? + rtcSetGeometryUserData(surface_geometry, user_data_map_[surface_geometry].get()); for (int i = 0; i < surface_triangles.size(); ++i) { - auto& triangle_ref = surface_data.prim_ref_buffer[i]; + auto& triangle_ref = surface_data->prim_ref_buffer[i]; // triangle_ref.embree_surface = embree_surface; } buffer_start += surface_triangles.size(); @@ -129,6 +131,7 @@ bool RayTracer::point_in_volume(TreeID scene, rayhit.ray.orientation = HitOrientation::ANY; rayhit.ray.set_tfar(INFTY); rayhit.ray.set_tnear(0.0); + if (exclude_primitives != nullptr) rayhit.ray.exclude_primitives = exclude_primitives; { From b56dcb91e838f78aa3cda844e1b29f1407b76d36 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 15 Feb 2024 17:33:40 -0600 Subject: [PATCH 09/12] Adding particle IDs. Setting random seed --- tools/particle_sim.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index f2da608..6d8db83 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -14,7 +14,7 @@ static double mfp {1.0}; struct Particle { -Particle(std::shared_ptr xdg) : xdg_(xdg) {} +Particle(std::shared_ptr xdg, uint32_t id) : xdg_(xdg), id_(id) {} void initialize() { // TODO: replace with sampling @@ -27,7 +27,7 @@ void initialize() { void surf_dist() { surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); if (surface_intersection_.second == ID_NONE) { - std::cerr << "Particle lost in volume " << volume_ << std::endl; + std::cerr << "Particle " << id_ << " lost in volume " << volume_ << std::endl; alive_ = false; return; } @@ -41,7 +41,7 @@ void sample_collision_distance() { void collide() { n_events_++; u_ = rand_dir(); - std::cout << "Particle collides with material at position " << r_ << ", new direction is " << u_ << std::endl; + std::cout << "Particle " << id_ << " collides with material at position " << r_ << ", new direction is " << u_ << std::endl; history_.clear(); } @@ -49,10 +49,10 @@ void advance() { if (collision_distance_ < surface_intersection_.first) { r_ += collision_distance_ * u_; - std::cout << "Particle collides with material at position " << r_ << std::endl; + std::cout << "Particle " << id_ << " collides with material at position " << r_ << std::endl; } else { r_ += surface_intersection_.first * u_; - std::cout << "Particle advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; + std::cout << "Particle " << id_ << " advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; } } @@ -61,14 +61,14 @@ void cross_surface() n_events_++; // check for the surface boundary condition if (xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION).value == "reflecting") { - std::cout << "Particle reflects off surface " << surface_intersection_.second << std::endl; + std::cout << "Particle " << id_ << " reflects off surface " << surface_intersection_.second << std::endl; Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); u_ = u_ - 2.0 * dot(u_, normal) * normal; u_ = u_.normalize(); if (history_.size() > 0) history_ = {history_.back()}; // reset to last intersection } else { volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); - std::cout << "Particle enters volume " << volume_ << std::endl; + std::cout << "Particle " << id_ << " enters volume " << volume_ << std::endl; if (volume_ == ID_NONE) { alive_ = false; return; @@ -78,7 +78,7 @@ void cross_surface() // Data Members std::shared_ptr xdg_; - +uint32_t id_ {0}; Position r_; Direction u_; MeshID volume_ {ID_NONE}; @@ -93,6 +93,8 @@ bool alive_ {true}; int main(int argc, char** argv) { +srand48(42); + // create a mesh manager std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); const auto& mm = xdg->mesh_manager(); @@ -105,7 +107,7 @@ mm->parse_metadata(); xdg->prepare_raytracer(); // create a new particle -Particle p(xdg); +Particle p(xdg, 1); const int max_events {100}; p.initialize(); From 9b2248aa0ec022b7fcaf1880579b226f7f0e4996 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 24 Feb 2024 14:15:12 -0600 Subject: [PATCH 10/12] Allowing embree version range --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9817b8f..a5a0e93 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ find_package(MOAB REQUIRED HINTS ${MOAB_DIR}) endif() # use Embree for CPU ray tracing -find_package(embree 4.0.0 REQUIRED) +find_package(embree 3.0.0...4.0.0 REQUIRED) if (NOT ${EMBREE_VERSION} VERSION_GREATER 3.6.0) message(FATAL_ERROR "XDG requires Embree v3.6.1 or higher.") endif() From 46802d78bf8a65b2f2f276fc60ca5e9025f1dba6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 24 Feb 2024 14:15:29 -0600 Subject: [PATCH 11/12] Improving particle simulation output --- tools/particle_sim.cpp | 83 +++++++++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 30 deletions(-) diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index 6d8db83..ade24da 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -3,6 +3,7 @@ #include #include +#include "xdg/error.h" #include "xdg/mesh_manager_interface.h" #include "xdg/moab/mesh_manager.h" #include "xdg/vec3da.h" @@ -14,7 +15,13 @@ static double mfp {1.0}; struct Particle { -Particle(std::shared_ptr xdg, uint32_t id) : xdg_(xdg), id_(id) {} +Particle(std::shared_ptr xdg, uint32_t id, bool verbose=true) : verbose_(verbose), xdg_(xdg), id_(id) {} + +template +void log (const std::string& msg, const Params&... fmt_args) { + if (!verbose_) return; + write_message(msg, fmt_args...); +} void initialize() { // TODO: replace with sampling @@ -26,12 +33,17 @@ void initialize() { void surf_dist() { surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); + if (surface_intersection_.first == 0.0) { + fatal_error("Particle {} stuck at position ({}, {}, {}) on surfacce {}", id_, r_.x, r_.y, r_.z, surface_intersection_.second); + alive_ = false; + return; + } if (surface_intersection_.second == ID_NONE) { std::cerr << "Particle " << id_ << " lost in volume " << volume_ << std::endl; alive_ = false; return; } - std::cout << "Intersected surface " << surface_intersection_.second << " at distance " << surface_intersection_.first << std::endl; + log("Intersected surface {} at distance {} ", surface_intersection_.second, surface_intersection_.first); } void sample_collision_distance() { @@ -41,7 +53,7 @@ void sample_collision_distance() { void collide() { n_events_++; u_ = rand_dir(); - std::cout << "Particle " << id_ << " collides with material at position " << r_ << ", new direction is " << u_ << std::endl; + 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(); } @@ -49,10 +61,11 @@ void advance() { if (collision_distance_ < surface_intersection_.first) { r_ += collision_distance_ * u_; - std::cout << "Particle " << id_ << " collides with material at position " << r_ << std::endl; + log("Particle {} collides with material at position ({}, {}, {}) ", id_, r_.x, r_.y, r_.z); + } else { r_ += surface_intersection_.first * u_; - std::cout << "Particle " << id_ << " advances to surface " << surface_intersection_.second << " at position " << r_ << std::endl; + log("Particle {} advances to surface {} at position ({}, {}, {}) ", id_, surface_intersection_.second, r_.x, r_.y, r_.z); } } @@ -61,14 +74,15 @@ void cross_surface() n_events_++; // check for the surface boundary condition if (xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION).value == "reflecting") { - std::cout << "Particle " << id_ << " reflects off surface " << surface_intersection_.second << std::endl; + log("Particle {} reflects off surface {}", id_, surface_intersection_.second); + Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); u_ = u_ - 2.0 * dot(u_, normal) * normal; u_ = u_.normalize(); if (history_.size() > 0) history_ = {history_.back()}; // reset to last intersection } else { volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); - std::cout << "Particle " << id_ << " enters volume " << volume_ << std::endl; + log("Particle {} enters volume {}", id_, volume_); if (volume_ == ID_NONE) { alive_ = false; return; @@ -77,6 +91,7 @@ void cross_surface() } // Data Members +bool verbose_ {true}; std::shared_ptr xdg_; uint32_t id_ {0}; Position r_; @@ -107,27 +122,35 @@ mm->parse_metadata(); xdg->prepare_raytracer(); // create a new particle -Particle p(xdg, 1); - -const int max_events {100}; -p.initialize(); -while (true) { - p.surf_dist(); - // terminate for leakage - if (!p.alive_) break; - p.sample_collision_distance(); - p.advance(); - if (p.surface_intersection_.first < p.collision_distance_) - p.cross_surface(); - else - p.collide(); - if (!p.alive_) break; - - if (p.n_events_ > max_events) { - std::cout << "Maximum number of events (" << max_events << ") reached" << std::endl; - break; - } -} -return 0; -} \ No newline at end of file +const int n_particles {100}; + +const int max_events {1000}; + +bool verbose = false; + + for (int i = 0; i < n_particles; i++) { + write_message("Starting particle {}", i); + Particle p(xdg, i, verbose); + p.initialize(); + while (true) { + p.surf_dist(); + // terminate for leakage + if (!p.alive_) break; + p.sample_collision_distance(); + p.advance(); + if (p.surface_intersection_.first < p.collision_distance_) + p.cross_surface(); + else + p.collide(); + if (!p.alive_) break; + + if (p.n_events_ > max_events) { + write_message("Maximum number of events ({}) reached", max_events); + break; + } + } + } + + return 0; +} From a472a3a040576480816355cda44b6aa83a5a0446 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 28 Feb 2024 12:31:28 -0600 Subject: [PATCH 12/12] Halting particle sim program on failed ray fire or zero distance intersection --- include/xdg/ray_tracing_interface.h | 2 +- include/xdg/xdg.h | 2 +- src/ray_tracing_interface.cpp | 3 ++- src/xdg.cpp | 2 +- tools/particle_sim.cpp | 22 ++++++++++++++++++---- 5 files changed, 23 insertions(+), 8 deletions(-) diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 6a01b48..756b555 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -40,7 +40,7 @@ class RayTracer { std::pair ray_fire(TreeID scene, const Position& origin, const Direction& direction, - const std::vector* exclude_primitives = nullptr); + std::vector* const exclude_primitives = nullptr); void closest(TreeID scene, const Position& origin, diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 1232884..197f223 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -36,7 +36,7 @@ bool point_in_volume(MeshID volume, std::pair ray_fire(MeshID volume, const Position& origin, const Direction& direction, - const std::vector* exclude_primitives = nullptr) const; + std::vector* 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 index e745c23..53aa23d 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -150,7 +150,7 @@ std::pair RayTracer::ray_fire(TreeID scene, const Position& origin, const Direction& direction, - const std::vector* exclude_primitves) + std::vector* 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 index f0f0c3b..c8c9493 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -59,7 +59,7 @@ std::pair XDG::ray_fire(MeshID volume, const Position& origin, const Direction& direction, - const std::vector* exclude_primitives) const + std::vector* 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 index ade24da..57e3985 100644 --- 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,6 +52,7 @@ 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(); @@ -59,6 +60,7 @@ void collide() { 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);