From cb6e6da8957dc5db480dd26d169856acee6dddbd Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 09:50:44 -0600 Subject: [PATCH 01/30] Account for rotation in hexagon-flat and annulus-flat geometries --- .../OptixCSP/src/core/CspElement.cpp | 13 ++- .../OptixCSP/src/shaders/GeometryDataST.h | 93 +++++++++++-------- .../OptixCSP/src/shaders/intersection.cu | 32 ++++--- 3 files changed, 82 insertions(+), 56 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 7ade3585..2e0b6302 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -316,19 +316,27 @@ GeometryDataST CspElement::toDeviceGeometryData() const if (aperture_type == ApertureType::HEXAGON) { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d vx = rotation_matrix.get_x_basis(); + Vec3d vy = rotation_matrix.get_y_basis(); + ApertureHexagon hex = static_cast(*m_aperture); float s = hex.get_side_length(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); if (surface_type == SurfaceType::FLAT) { - GeometryDataST::Hexagon_Flat hex(o, n, s); + GeometryDataST::Hexagon_Flat hex(o, n, vx, vy, s); geometry_data.setHexagon_Flat(hex); } } if (aperture_type == ApertureType::ANNULUS) { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d vx = rotation_matrix.get_x_basis(); + Vec3d vy = rotation_matrix.get_y_basis(); + ApertureAnnulus anf = static_cast(*m_aperture); float radius_in = anf.get_radius_inner(); float radius_out = anf.get_radius_outer(); @@ -337,7 +345,8 @@ GeometryDataST CspElement::toDeviceGeometryData() const float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); if (surface_type == SurfaceType::FLAT) { - GeometryDataST::Annulus_Flat anf(o, n, radius_in, radius_out, arc); + GeometryDataST::Annulus_Flat anf(o, n, vx, vy, + radius_in, radius_out, arc); geometry_data.setAnnulus_Flat(anf); } } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index 217e64dd..da48d6b4 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -25,7 +25,7 @@ namespace OptixCSP QUADRILATERAL_FLAT = 6, CIRCLE_FLAT = 7, HEXAGON_FLAT = 8, - ANNULUS_FLAT = 9 + ANNULUS_FLAT = 9 }; struct Parallelogram @@ -40,10 +40,10 @@ namespace OptixCSP this->v2 *= 1.0f / dot(v2, v2); plane = make_float4(normal, d); } - float4 plane; - float3 v1; - float3 v2; - float3 anchor; + float4 plane; // plane equation: (normal, dot(anchor, normal)) + float3 v1; // edge vector 1, stored as v1/dot(v1,v1) + float3 v2; // edge vector 2, stored as v2/dot(v2,v2) + float3 anchor; // corner point of the parallelogram }; // same as parallelogram, however defined with different attributes @@ -58,12 +58,12 @@ namespace OptixCSP plane = make_float4(normal, d); } - float4 plane; - float3 center; - float3 x; - float3 y; - float width; - float height; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // center in global coordinates + float3 x; // local x axis unit vector + float3 y; // local y axis unit vector + float width; // full width along x + float height; // full height along y }; struct Cylinder_Y @@ -75,11 +75,11 @@ namespace OptixCSP assert(dot(base_x, base_z) < 1e-3f); } - float3 center; - float radius; - float half_height; - float3 base_x; // x axis of the cylinder - float3 base_z; // z axis of the cylinder + float3 center; // center of the cylinder in global coordinates + float radius; // radius + float half_height; // half the height along the local y axis + float3 base_x; // x axis of the cylinder + float3 base_z; // z axis of the cylinder }; struct Rectangle_Parabolic @@ -96,13 +96,13 @@ namespace OptixCSP plane = make_float4(normal, d); } - float4 plane; - float3 v1; - float3 v2; - float3 anchor; + float4 plane; // plane equation of the base rectangle: (normal, dot(anchor, normal)) + float3 v1; // edge vector 1, stored as v1/dot(v1,v1) + float3 v2; // edge vector 2, stored as v2/dot(v2,v2) + float3 anchor; // corner point of the base rectangle // float3 focus; - float curv_x; - float curv_y; + float curv_x; // curvature along local x axis + float curv_y; // curvature along local y axis }; struct Triangle_Flat @@ -116,8 +116,8 @@ namespace OptixCSP } float3 v0; // base vertex float3 e1, e2; // edges - float3 normal; - float d; // plane distance + float3 normal; // normal unit vector + float d; // plane distance }; struct Quadrilateral_Flat @@ -142,40 +142,51 @@ namespace OptixCSP Circle_Flat(const float3 &origin, const float3 &normal, const float &radius) : r(radius), center(origin) { - plane = make_float4(normalize(normal), dot(center, normal)); + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); } - float4 plane; - float3 center; - float r; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float r; // radius }; struct Hexagon_Flat { Hexagon_Flat() = default; - Hexagon_Flat(const float3 &origin, const float3 &normal, const float &side_length) - : s(side_length), center(origin) + Hexagon_Flat(const float3 &origin, const float3 &normal, + const float3 &x_ax, const float3 &y_ax, + const float &side_length) + : center(origin), x_axis(x_ax), y_axis(y_ax), s(side_length) { - plane = make_float4(normalize(normal), dot(center, normal)); + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); } - float4 plane; - float3 center; - float s; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float3 x_axis; // unit vector + float3 y_axis; // unit vector + float s; // side length }; struct Annulus_Flat { Annulus_Flat() = default; Annulus_Flat(const float3 &origin, const float3 &normal, + const float3 &x_ax, const float3 &y_ax, const float &r_inner, const float &r_outer, const float &arc) - : center(origin), ri(r_inner), ro(r_outer), arc(arc) + : center(origin), x_axis(x_ax), y_axis(y_ax), + ri(r_inner), ro(r_outer), arc(arc) { - plane = make_float4(normalize(normal), dot(center, normal)); + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); } - float4 plane; - float3 center; - float ri; - float ro; - float arc; // Arc angle in radians with x-axis in the middle + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float3 x_axis; // local x axis unit vector (arc is centered about this axis) + float3 y_axis; // local y axis unit vector + float ri; // inner radius + float ro; // outer radius + float arc; // total arc angle in radians, centered on x_axis }; GeometryDataST() = default; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index b9688d39..b225dc1e 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -624,6 +624,9 @@ extern "C" __global__ void __intersection__circle_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { + // Since everything is in global coordinates (e.g., the ray intersection coordinates), + // and the circle is rotationally symmetric, we don't need to worry about + // any rotation of the circle in local coordinates float3 p = ray_orig + ray_dir * t; float d = length(p - circ.center); if (d <= circ.r) @@ -652,34 +655,35 @@ extern "C" __global__ void __intersection__hexagon_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { - // TODO: Need to adjust for possible rotation... bool is_in = false; float3 p = ray_orig + ray_dir * t - hex.center; - // float d = length(p - circ.center); - float s = hex.s; - float xl = 0.5f * s; - float yl = 0.5f * sqrtf(3.0f) * s; - if (-xl <= p.x && p.x <= xl && -yl <= p.y && p.y <= yl) + // Project onto the local x and y axes which are unit vectors + const float px = dot(p, hex.x_axis); + const float py = dot(p, hex.y_axis); + const float s = hex.s; + const float xl = 0.5f * s; + const float yl = 0.5f * sqrtf(3.0f) * s; + if (-xl <= px && px <= xl && -yl <= py && py <= yl) { // Center is_in = true; } - else if (-s <= p.x && p.x < xl) + else if (-s <= px && px < xl) { // Left side - float y1 = sqrtf(3.0f) * (p.x + s); + float y1 = sqrtf(3.0f) * (px + s); float y2 = -y1; - if (y2 <= p.y && p.y <= y1) + if (y2 <= py && py <= y1) { is_in = true; } } - else if (xl < p.x && p.x <= s) + else if (xl < px && px <= s) { // Right side - float y1 = sqrtf(3.0f) * (p.x - s); + float y1 = sqrtf(3.0f) * (px - s); float y2 = -y1; - if (y1 <= p.y && p.y <= y2) + if (y1 <= py && py <= y2) { is_in = true; } @@ -716,7 +720,9 @@ extern "C" __global__ void __intersection__annulus_flat() float d = length(p); if (anf.ri <= d && d <= anf.ro) { - float theta = atan2f(p.y, p.x); + float px = dot(p, anf.x_axis); + float py = dot(p, anf.y_axis); + float theta = atan2f(py, px); if (fabsf(theta) <= 0.5f * anf.arc) { optixReportIntersection(t, From de9984c8560f86d88d4ca020a84949c8c6a91ac6 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 10:09:12 -0600 Subject: [PATCH 02/30] Update tests for rotated elements --- .../geometry_intersection_test.cpp | 130 ++++++++++++++---- 1 file changed, 100 insertions(+), 30 deletions(-) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index b1dd3f36..c02bcfba 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -14,7 +14,8 @@ const uint_fast64_t NRAYS = 10000; element_id set_default_sd(SimulationData &sd, surface_ptr surf, - aperture_ptr ap) + aperture_ptr ap, + double rotation) { sd.clear(); @@ -35,6 +36,7 @@ element_id set_default_sd(SimulationData &sd, el->set_front_optical_properties(el_optics); el->set_back_optical_properties(el_optics); el->set_name("el"); + el->set_zrot(rotation); // Add element to stage element_id id = sd.add_element(el); @@ -44,8 +46,8 @@ element_id set_default_sd(SimulationData &sd, element_ptr stop = make_element(); double xlb, xub, ylb, yub; ap->bounding_box(xlb, xub, ylb, yub); - const double sx = std::max(fabs(xlb), fabs(xub)) + 1.0; - const double sy = std::max(fabs(ylb), fabs(yub)) + 1.0; + const double sx = std::max(fabs(xlb), fabs(xub)) + 2.0; + const double sy = std::max(fabs(ylb), fabs(yub)) + 2.0; stop->set_origin(0, 0, Z_BACKSTOP); stop->set_aim_vector(0, 0, 100); stop->set_surface(make_surface()); @@ -67,11 +69,12 @@ TEST(OptixRunner, FlatRectangle) { const double XL = 10.0; const double YL = 5.0; + const double ROT_DEG = -10.0; auto surf = make_surface(); auto aper = make_aperture(XL, YL); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -86,6 +89,8 @@ TEST(OptixRunner, FlatRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -96,20 +101,22 @@ TEST(OptixRunner, FlatRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { // We hit the test element. Check that the height is as expected EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; // And that we are in the aperture - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { // We hit the back stop element. Check that the height is as expected EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); // And that we are not in the aperture. - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -117,11 +124,12 @@ TEST(OptixRunner, FlatRectangle) TEST(OptixRunner, FlatEquilateralTriangle) { const double d = 4.0; + const double ROT_DEG = 90.0; auto surf = make_surface(); auto aper = make_aperture(d); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -136,6 +144,8 @@ TEST(OptixRunner, FlatEquilateralTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -146,16 +156,18 @@ TEST(OptixRunner, FlatEquilateralTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -164,11 +176,12 @@ TEST(OptixRunner, FlatTriangle) { const double x1 = 0.0, x2 = 1.0, x3 = 2.0 * x2; const double y1 = 0.0, y2 = 2.0, y3 = y1; + const double ROT_DEG = 110.0; auto surf = make_surface(); auto aper = make_aperture(x1, y1, x2, y2, x3, y3); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -183,6 +196,8 @@ TEST(OptixRunner, FlatTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -193,16 +208,18 @@ TEST(OptixRunner, FlatTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -212,12 +229,13 @@ TEST(OptixRunner, FlatQuadrilateral) // Parallelogram const double x1 = 0.0, x2 = 3.0, x3 = (x2 - x1) + 1.0, x4 = x3 - x2 + x1; const double y1 = 0.0, y2 = y1, y3 = 2.0, y4 = y3; + const double ROT_DEG = -45.0; auto surf = make_surface(); auto aper = make_aperture( x1, y1, x2, y2, x3, y3, x4, y4); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -232,6 +250,8 @@ TEST(OptixRunner, FlatQuadrilateral) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -242,16 +262,18 @@ TEST(OptixRunner, FlatQuadrilateral) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -263,11 +285,12 @@ TEST(OptixRunner, ParabolaRectangle) constexpr double FX = 0.5 / CX; constexpr double FY = 0.5 / CY; const double XL = 10.0, YL = 5.0; + const double ROT_DEG = -135.0; auto surf = make_surface(FX, FY); auto aper = make_aperture(XL, YL); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -282,6 +305,8 @@ TEST(OptixRunner, ParabolaRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -292,17 +317,19 @@ TEST(OptixRunner, ParabolaRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { - const double z1 = Z_ELEM + 0.5 * CX * p1[0] * p1[0] + 0.5 * CY * p1[1] * p1[1]; + const double z1 = Z_ELEM + 0.5 * CX * lx * lx + 0.5 * CY * ly * ly; EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -311,11 +338,12 @@ TEST(OptixRunner, Cylinder) { const double R = 5.0; const double YL = 3.0; // Total cylinder length + const double ROT_DEG = 25.0; auto surf = make_surface(R); auto aper = make_aperture(2 * R, YL); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -330,6 +358,8 @@ TEST(OptixRunner, Cylinder) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -340,17 +370,19 @@ TEST(OptixRunner, Cylinder) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { - const double z1 = Z_ELEM + sqrt(R * R - p1[0] * p1[0]); + const double z1 = Z_ELEM + sqrt(R * R - lx * lx); EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } @@ -358,11 +390,12 @@ TEST(OptixRunner, Cylinder) TEST(OptixRunner, FlatCircle) { const double R = 5.0; + const double ROT_DEG = 10.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(2 * R); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -404,11 +437,12 @@ TEST(OptixRunner, FlatCircle) TEST(OptixRunner, FlatHexagon) { const double S = 5.0; + const double ROT_DEG = 30.0; auto surf = make_surface(); auto aper = make_aperture(2 * S); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -423,6 +457,8 @@ TEST(OptixRunner, FlatHexagon) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -433,30 +469,33 @@ TEST(OptixRunner, FlatHexagon) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); } } } -TEST(OptixRunner, FlatAnnulus) +TEST(OptixRunner, FlatAnnulus_FullArc) { const double R0 = 5.0; - const double R1 = 10.0; + const double R1 = 180.0; const double ARC = 2 * PI; + const double ROT_DEG = -15.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(R0, R1, ARC); SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -494,3 +533,34 @@ TEST(OptixRunner, FlatAnnulus) } } } + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * PI / 180.0); + const double sin_rot = sin(ROT_DEG * PI / 180.0); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot + p1[1] * sin_rot; + const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + } + } +} + From 059f44b1fd98b3cdf91e4e87a3550f4be0ab9620 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 10:35:42 -0600 Subject: [PATCH 03/30] Fix dumb AI error... --- .../geometry_intersection_test.cpp | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index c02bcfba..e53f3897 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -534,6 +534,29 @@ TEST(OptixRunner, FlatAnnulus_FullArc) } } +TEST(OptixRunner, FlatAnnulus_PartialArc) +{ + const double R0 = 5.0; + const double R1 = 180.0; + const double ARC = 0.5 * PI; + const double ROT_DEG = -15.0; + auto surf = make_surface(); + auto aper = make_aperture(R0, R1, ARC); + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); const double cos_rot = cos(ROT_DEG * PI / 180.0); From da7a9e246f459da34828f37db002c8ed03a60bcb Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 10:46:23 -0600 Subject: [PATCH 04/30] Fix missed conversions to float3 --- .../optix_runner/OptixCSP/src/core/CspElement.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 2e0b6302..3adcba4d 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -317,8 +317,8 @@ GeometryDataST CspElement::toDeviceGeometryData() const if (aperture_type == ApertureType::HEXAGON) { Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d vx = rotation_matrix.get_x_basis(); - Vec3d vy = rotation_matrix.get_y_basis(); + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); ApertureHexagon hex = static_cast(*m_aperture); float s = hex.get_side_length(); @@ -334,8 +334,8 @@ GeometryDataST CspElement::toDeviceGeometryData() const if (aperture_type == ApertureType::ANNULUS) { Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d vx = rotation_matrix.get_x_basis(); - Vec3d vy = rotation_matrix.get_y_basis(); + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); ApertureAnnulus anf = static_cast(*m_aperture); float radius_in = anf.get_radius_inner(); From 94d6cd82015b64227254a2afc97df8cac77a222d Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 11:09:44 -0600 Subject: [PATCH 05/30] Fix inverted coordinate transform --- .../geometry_intersection_test.cpp | 66 ++++++++++--------- 1 file changed, 34 insertions(+), 32 deletions(-) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index e53f3897..b916ff74 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -1,11 +1,13 @@ #include +#include #include #include #include #include using SolTrace::Runner::RunnerStatus; +using SolTrace::Data::D2R; const double Z_ELEM = 50.0; const double Z_BACKSTOP = Z_ELEM - 0.5 * Z_ELEM; @@ -89,8 +91,8 @@ TEST(OptixRunner, FlatRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -101,8 +103,8 @@ TEST(OptixRunner, FlatRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -144,8 +146,8 @@ TEST(OptixRunner, FlatEquilateralTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -156,8 +158,8 @@ TEST(OptixRunner, FlatEquilateralTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -196,8 +198,8 @@ TEST(OptixRunner, FlatTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -208,8 +210,8 @@ TEST(OptixRunner, FlatTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -250,8 +252,8 @@ TEST(OptixRunner, FlatQuadrilateral) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -262,8 +264,8 @@ TEST(OptixRunner, FlatQuadrilateral) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -305,8 +307,8 @@ TEST(OptixRunner, ParabolaRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -317,8 +319,8 @@ TEST(OptixRunner, ParabolaRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -358,8 +360,8 @@ TEST(OptixRunner, Cylinder) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -370,8 +372,8 @@ TEST(OptixRunner, Cylinder) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -457,8 +459,8 @@ TEST(OptixRunner, FlatHexagon) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -469,8 +471,8 @@ TEST(OptixRunner, FlatHexagon) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { @@ -559,8 +561,8 @@ TEST(OptixRunner, FlatAnnulus_PartialArc) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); - const double cos_rot = cos(ROT_DEG * PI / 180.0); - const double sin_rot = sin(ROT_DEG * PI / 180.0); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -571,8 +573,8 @@ TEST(OptixRunner, FlatAnnulus_PartialArc) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot + p1[1] * sin_rot; - const double ly = -p1[0] * sin_rot + p1[1] * cos_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { From 32f34af0912819bda9a8398f89dd337ff9491c81 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 12:23:23 -0600 Subject: [PATCH 06/30] Fix tests and make them more robust --- .../geometry_intersection_test.cpp | 190 ++++++++++++++++-- 1 file changed, 174 insertions(+), 16 deletions(-) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index b916ff74..70dbf0e4 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -6,14 +6,53 @@ #include #include -using SolTrace::Runner::RunnerStatus; using SolTrace::Data::D2R; +using SolTrace::Runner::RunnerStatus; const double Z_ELEM = 50.0; const double Z_BACKSTOP = Z_ELEM - 0.5 * Z_ELEM; const double TOL = 1e-6; const uint_fast64_t NRAYS = 10000; +// Adapted from CspElement::set_bounding_box_local. +// Computes the global axis-aligned bounding box for a local bounding box +// transformed by a z-axis rotation (in degrees) and origin offset. +static void compute_global_bounding_box( + const glm::dvec3 &lower_local, + const glm::dvec3 &upper_local, + double zrot_deg, // z-axis rotation in degrees (local-to-global) + const glm::dvec3 &origin, + glm::dvec3 &lower_global, + glm::dvec3 &upper_global) +{ + const double cr = cos(zrot_deg * D2R); + const double sr = sin(zrot_deg * D2R); + // glm is column-major: columns are (cr, sr, 0), (-sr, cr, 0), (0, 0, 1) + const glm::dmat3 rotation(cr, sr, 0.0, + -sr, cr, 0.0, + 0.0, 0.0, 1.0); + const glm::dvec3 c0 = lower_local; + const glm::dvec3 c7 = upper_local; + const glm::dvec3 corners[8] = { + rotation * glm::dvec3(c0[0], c0[1], c0[2]) + origin, + rotation * glm::dvec3(c0[0], c0[1], c7[2]) + origin, + rotation * glm::dvec3(c0[0], c7[1], c0[2]) + origin, + rotation * glm::dvec3(c7[0], c0[1], c0[2]) + origin, + rotation * glm::dvec3(c0[0], c7[1], c7[2]) + origin, + rotation * glm::dvec3(c7[0], c0[1], c7[2]) + origin, + rotation * glm::dvec3(c7[0], c7[1], c0[2]) + origin, + rotation * glm::dvec3(c7[0], c7[1], c7[2]) + origin, + }; + + lower_global = glm::dvec3(std::numeric_limits::max()); + upper_global = glm::dvec3(std::numeric_limits::lowest()); + for (const auto &c : corners) + { + lower_global = glm::min(lower_global, c); + upper_global = glm::max(upper_global, c); + } +} + element_id set_default_sd(SimulationData &sd, surface_ptr surf, aperture_ptr ap, @@ -46,10 +85,20 @@ element_id set_default_sd(SimulationData &sd, // Back stop element that is bigger than the created element so that the // testing element casts a shadow on this big thing. element_ptr stop = make_element(); - double xlb, xub, ylb, yub; + double xlb, xub, ylb, yub, zlb, zub; ap->bounding_box(xlb, xub, ylb, yub); - const double sx = std::max(fabs(xlb), fabs(xub)) + 2.0; - const double sy = std::max(fabs(ylb), fabs(yub)) + 2.0; + surf->bounding_box(xlb, xub, ylb, yub, zlb, zub); + + // Transform local AABB into global AABB accounting for zrot + glm::dvec3 lower_global, upper_global; + compute_global_bounding_box(glm::dvec3(xlb, ylb, zlb), + glm::dvec3(xub, yub, zub), + rotation, + glm::dvec3(0.0, 0.0, Z_ELEM), + lower_global, upper_global); + + const double sx = std::max(fabs(lower_global[0]), fabs(upper_global[0])) + 2.0; + const double sy = std::max(fabs(lower_global[1]), fabs(upper_global[1])) + 2.0; stop->set_origin(0, 0, Z_BACKSTOP); stop->set_aim_vector(0, 0, 100); stop->set_surface(make_surface()); @@ -75,6 +124,8 @@ TEST(OptixRunner, FlatRectangle) auto surf = make_surface(); auto aper = make_aperture(XL, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -103,7 +154,7 @@ TEST(OptixRunner, FlatRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) @@ -112,6 +163,8 @@ TEST(OptixRunner, FlatRectangle) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; // And that we are in the aperture EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { @@ -119,8 +172,15 @@ TEST(OptixRunner, FlatRectangle) EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); // And that we are not in the aperture. EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatEquilateralTriangle) @@ -130,6 +190,8 @@ TEST(OptixRunner, FlatEquilateralTriangle) auto surf = make_surface(); auto aper = make_aperture(d); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -158,20 +220,29 @@ TEST(OptixRunner, FlatEquilateralTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatTriangle) @@ -182,6 +253,8 @@ TEST(OptixRunner, FlatTriangle) auto surf = make_surface(); auto aper = make_aperture(x1, y1, x2, y2, x3, y3); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -210,20 +283,29 @@ TEST(OptixRunner, FlatTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatQuadrilateral) @@ -236,6 +318,8 @@ TEST(OptixRunner, FlatQuadrilateral) auto aper = make_aperture( x1, y1, x2, y2, x3, y3, x4, y4); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -264,20 +348,29 @@ TEST(OptixRunner, FlatQuadrilateral) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, ParabolaRectangle) @@ -291,6 +384,8 @@ TEST(OptixRunner, ParabolaRectangle) auto surf = make_surface(FX, FY); auto aper = make_aperture(XL, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -319,7 +414,7 @@ TEST(OptixRunner, ParabolaRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) @@ -327,13 +422,22 @@ TEST(OptixRunner, ParabolaRectangle) const double z1 = Z_ELEM + 0.5 * CX * lx * lx + 0.5 * CY * ly * ly; EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, Cylinder) @@ -344,6 +448,8 @@ TEST(OptixRunner, Cylinder) auto surf = make_surface(R); auto aper = make_aperture(2 * R, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -372,7 +478,7 @@ TEST(OptixRunner, Cylinder) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) @@ -380,22 +486,33 @@ TEST(OptixRunner, Cylinder) const double z1 = Z_ELEM + sqrt(R * R - lx * lx); EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatCircle) { const double R = 5.0; - const double ROT_DEG = 10.0; // Should make no difference + const double ROT_DEG = 10.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(2 * R); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -427,13 +544,22 @@ TEST(OptixRunner, FlatCircle) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + ++hits; + if (!aper->is_in(p1[0], p1[1])) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + ++misses; + if (aper->is_in(p1[0], p1[1])) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatHexagon) @@ -443,6 +569,8 @@ TEST(OptixRunner, FlatHexagon) auto surf = make_surface(); auto aper = make_aperture(2 * S); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -471,31 +599,42 @@ TEST(OptixRunner, FlatHexagon) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatAnnulus_FullArc) { const double R0 = 5.0; const double R1 = 180.0; - const double ARC = 2 * PI; + const double ARC = 360.0; const double ROT_DEG = -15.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(R0, R1, ARC); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -527,24 +666,35 @@ TEST(OptixRunner, FlatAnnulus_FullArc) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + ++hits; + if (!aper->is_in(p1[0], p1[1])) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + ++misses; + if (aper->is_in(p1[0], p1[1])) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatAnnulus_PartialArc) { const double R0 = 5.0; const double R1 = 180.0; - const double ARC = 0.5 * PI; + const double ARC = 90.0; const double ROT_DEG = -15.0; auto surf = make_surface(); auto aper = make_aperture(R0, R1, ARC); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; @@ -573,19 +723,27 @@ TEST(OptixRunner, FlatAnnulus_PartialArc) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; - const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } - From 9d7e1aea9aa1cd7685c5de5ab51315aee3337668 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 13:20:35 -0600 Subject: [PATCH 07/30] Final flat hexagon fix --- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index b225dc1e..9fc97e66 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -668,7 +668,7 @@ extern "C" __global__ void __intersection__hexagon_flat() // Center is_in = true; } - else if (-s <= px && px < xl) + else if (-s <= px && px < -xl) { // Left side float y1 = sqrtf(3.0f) * (px + s); From 7c33c4bfec1af00bd07f6fb0953aa21b3aa1a385 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 15:54:39 -0600 Subject: [PATCH 08/30] Add circle parabolic element to optix runner --- .../OptixCSP/src/core/CspElement.cpp | 11 ++ .../OptixCSP/src/core/geometry_manager.cpp | 4 + .../OptixCSP/src/core/pipeline_manager.cpp | 7 +- .../OptixCSP/src/core/soltrace_type.h | 30 ++-- .../OptixCSP/src/shaders/GeometryDataST.h | 54 +++++-- .../OptixCSP/src/shaders/Soltrace.h | 21 ++- .../OptixCSP/src/shaders/intersection.cu | 146 ++++++++++++++++-- 7 files changed, 227 insertions(+), 46 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 3adcba4d..38843054 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -307,11 +307,22 @@ GeometryDataST CspElement::toDeviceGeometryData() const float r = circ.get_radius(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { GeometryDataST::Circle_Flat heliostat(o, n, r); geometry_data.setCircle_Flat(heliostat); } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 v1 = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 v2 = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Circle_Parabolic heliostat(o, v1, v2, cx, cy, r); + } } if (aperture_type == ApertureType::HEXAGON) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index 118f69b9..6e07ac71 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -63,6 +63,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::CIRCLE_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::CIRCLE_PARABOLIC); + } else { std::stringstream ss; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp index ade2345d..7f6b40ca 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp @@ -45,7 +45,12 @@ const std::map IntersectionKernelMap = { {OpticalEntityType::QUADRILATERAL_FLAT, "__intersection__quadrilateral_flat"}, {OpticalEntityType::CIRCLE_FLAT, "__intersection__circle_flat"}, {OpticalEntityType::HEXAGON_FLAT, "__intersection__hexagon_flat"}, - {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}}; + {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}, + {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}, + {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, + {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, + {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, + {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; pipelineManager::pipelineManager(SoltraceState &state) : m_state(state) {} diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h index d98569a8..ff8703fa 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h @@ -31,12 +31,12 @@ namespace OptixCSP SurfaceType surfaceType; ApertureType apertureType; - SurfaceApertureMap() : surfaceType(SurfaceType::FLAT), - apertureType(ApertureType::RECTANGLE) - { - } + SurfaceApertureMap() : surfaceType(SurfaceType::FLAT), + apertureType(ApertureType::RECTANGLE) + { + } - SurfaceApertureMap(SurfaceType surf, ApertureType ap) : surfaceType(surf), + SurfaceApertureMap(SurfaceType surf, ApertureType ap) : surfaceType(surf), apertureType(ap) { } @@ -48,19 +48,19 @@ namespace OptixCSP return (surfaceType == map.surfaceType) && (apertureType == map.apertureType); } - bool operator<(const SurfaceApertureMap &b) const - { - return surfaceType < b.surfaceType || - (surfaceType == b.surfaceType && apertureType < b.apertureType); - } + bool operator<(const SurfaceApertureMap &b) const + { + return surfaceType < b.surfaceType || + (surfaceType == b.surfaceType && apertureType < b.apertureType); + } - friend std::ostream& operator<<(std::ostream &os, const SurfaceApertureMap &sam); + friend std::ostream &operator<<(std::ostream &os, const SurfaceApertureMap &sam); }; - inline std::ostream& operator<<(std::ostream &os, const SurfaceApertureMap &sam) + inline std::ostream &operator<<(std::ostream &os, const SurfaceApertureMap &sam) { - os << "SurfApMap -- Surface: " << static_cast(sam.surfaceType) - << " Aperture: " << static_cast(sam.apertureType); - return os; + os << "SurfApMap -- Surface: " << static_cast(sam.surfaceType) + << " Aperture: " << static_cast(sam.apertureType); + return os; } } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index da48d6b4..149745e1 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -75,11 +75,11 @@ namespace OptixCSP assert(dot(base_x, base_z) < 1e-3f); } - float3 center; // center of the cylinder in global coordinates - float radius; // radius - float half_height; // half the height along the local y axis - float3 base_x; // x axis of the cylinder - float3 base_z; // z axis of the cylinder + float3 center; // center of the cylinder in global coordinates + float radius; // radius + float half_height; // half the height along the local y axis + float3 base_x; // x axis of the cylinder + float3 base_z; // z axis of the cylinder }; struct Rectangle_Parabolic @@ -101,8 +101,8 @@ namespace OptixCSP float3 v2; // edge vector 2, stored as v2/dot(v2,v2) float3 anchor; // corner point of the base rectangle // float3 focus; - float curv_x; // curvature along local x axis - float curv_y; // curvature along local y axis + float curv_x; // curvature along local x axis + float curv_y; // curvature along local y axis }; struct Triangle_Flat @@ -145,9 +145,9 @@ namespace OptixCSP const float3 n = normalize(normal); plane = make_float4(n, dot(center, n)); } - float4 plane; // normal unit vector, dot(center, normal) - float3 center; // local origin in global coordinates - float r; // radius + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float r; // radius }; struct Hexagon_Flat @@ -173,7 +173,8 @@ namespace OptixCSP Annulus_Flat() = default; Annulus_Flat(const float3 &origin, const float3 &normal, const float3 &x_ax, const float3 &y_ax, - const float &r_inner, const float &r_outer, const float &arc) + const float &r_inner, const float &r_outer, + const float &arc) : center(origin), x_axis(x_ax), y_axis(y_ax), ri(r_inner), ro(r_outer), arc(arc) { @@ -189,6 +190,23 @@ namespace OptixCSP float arc; // total arc angle in radians, centered on x_axis }; + struct Circle_Parabolic + { + Circle_Parabolic() = default; + Circle_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, const float &r) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), radius(r) + { + } + float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; + float radius; + }; + GeometryDataST() = default; void setParallelogram(const Parallelogram &p) @@ -308,6 +326,19 @@ namespace OptixCSP return annulus_flat; } + void setCircle_Parabolic(const Circle_Parabolic &circp) + { + assert(type == UNKNOWN_TYPE); + type = CIRCLE_PARABOLIC; + circle_parabolic = circp; + } + + __host__ __device__ const Circle_Parabolic &getCircle_Parabolic() const + { + assert(type == CIRCLE_PARABOLIC); + return circle_parabolic; + } + Type type = UNKNOWN_TYPE; int32_t id = OptixCSP::kElementIdUnassigned; @@ -324,6 +355,7 @@ namespace OptixCSP Circle_Flat circle_flat; Hexagon_Flat hexagon_flat; Annulus_Flat annulus_flat; + Circle_Parabolic circle_parabolic; }; }; } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h index 1c4c87ab..5426b0e6 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h @@ -36,14 +36,19 @@ namespace OptixCSP{ }; enum OpticalEntityType : unsigned int { - RECTANGLE_FLAT = 0, - RECTANGLE_PARABOLIC = 1, - CYLINDRICAL = 2, - TRIANGLE_FLAT = 3, - QUADRILATERAL_FLAT = 4, - CIRCLE_FLAT = 5, - HEXAGON_FLAT = 6, - ANNULUS_FLAT = 7, + RECTANGLE_FLAT = 0, + RECTANGLE_PARABOLIC = 1, + CYLINDRICAL = 2, + TRIANGLE_FLAT = 3, + QUADRILATERAL_FLAT = 4, + CIRCLE_FLAT = 5, + HEXAGON_FLAT = 6, + ANNULUS_FLAT = 7, + CIRCLE_PARABOLIC = 8, + HEXAGON_PARABOLIC = 9, + TRIANGLE_PARABOLIC = 10, + ANNULUS_PARABOLIC = 11, + QUADRILATERAL_PARABOLIC = 12, NUM_OPTICAL_ENTITY_TYPES }; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 9fc97e66..58a82fdb 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -380,19 +380,22 @@ extern "C" __global__ void __intersection__rectangle_parabolic() float t2 = 0.0f; float x_hit = 0.0f; float y_hit = 0.0f; - const float eps = 1e-12f; + const float eps = 1e-6f; bool valid = false; - if (fabsf(A) < eps) { - if (fabsf(B) > eps) { + if (fabsf(A) < eps) + { + if (fabsf(B) > eps) + { t1 = -C / B; const float p1x = ox + t1 * dx; const float p1y = oy + t1 * dy; const float a1 = p1x / (L1 / 2.0f); const float a2 = p1y / (L2 / 2.0f); - if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax - && !(a1 < -1.0f || a1 > 1.0f || a2 < -1.0f || a2 > 1.0f)) { + if ((t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax) && + !(a1 < -1.0f || a1 > 1.0f || a2 < -1.0f || a2 > 1.0f)) + { t = t1; x_hit = p1x; y_hit = p1y; @@ -418,15 +421,15 @@ extern "C" __global__ void __intersection__rectangle_parabolic() const float a1_2 = p2x / (L1 / 2.0f); const float a2_2 = p2y / (L2 / 2.0f); - if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax - && !(a1_1 < -1.0f || a1_1 > 1.0f || a2_1 < -1.0f || a2_1 > 1.0f)) { + if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax && !(a1_1 < -1.0f || a1_1 > 1.0f || a2_1 < -1.0f || a2_1 > 1.0f)) + { t = t1; x_hit = p1x; y_hit = p1y; valid = true; } - else if (t2 > 0.0f && t2 >= ray_tmin && t2 <= ray_tmax - && !(a1_2 < -1.0f || a1_2 > 1.0f || a2_2 < -1.0f || a2_2 > 1.0f)) { + else if (t2 > 0.0f && t2 >= ray_tmin && t2 <= ray_tmax && !(a1_2 < -1.0f || a1_2 > 1.0f || a2_2 < -1.0f || a2_2 > 1.0f)) + { t = t2; x_hit = p2x; y_hit = p2y; @@ -435,7 +438,8 @@ extern "C" __global__ void __intersection__rectangle_parabolic() } } - if (!valid) { + if (!valid) + { return; } @@ -625,7 +629,7 @@ extern "C" __global__ void __intersection__circle_flat() if (t > ray_tmin && t < ray_tmax) { // Since everything is in global coordinates (e.g., the ray intersection coordinates), - // and the circle is rotationally symmetric, we don't need to worry about + // and the circle is rotationally symmetric, we don't need to worry about // any rotation of the circle in local coordinates float3 p = ray_orig + ray_dir * t; float d = length(p - circ.center); @@ -734,3 +738,123 @@ extern "C" __global__ void __intersection__annulus_flat() } } } + +extern "C" __global__ void __intersection__circle_parabolic() +{ + const OptixCSP::GeometryDataST::Circle_Parabolic &circp = params.geoemetry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(), ray_tmax = optixGetRayTmax(); + + const float eps = 1e-6f; + + const float3 center = circp.center; + const float cx = circp.cx; + const float cy = circp.cy; + const float r = circp.radius; + + // TODO: Compute the local ray direction and ray origin + const float3 x_ax = circp.x_axis; + const float3 y_ax = circp.y_axis; + const float3 n = normalize(cross(e2, e1)); + + float3 d = ray_orig - center; + float ox = dot(d, x_ax); + float oy = dot(d, y_ax); + float oz = dot(d, n); + + float dx = dot(ray_dir, x_ax); + float dy = dot(ray_dir, y_ax); + float dz = dot(ray_dir, n); + + float A = (0.5f * cx) * (dx * dx) + (0.5f * cy) * (dy * dy); + float B = cx * (ox * dx) + cy * (oy * dy) - dz; + float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; + + float t = 0.0f; + float t1 = 0.0f; + float t2 = 0.0f; + + if (fabsf(A) < eps) + { + t1 = -C / B; + const float p1x = ox + t1 * dx; + const float p1y = oy + t1 * dy; + const float v = sqrtf(p1x * p1x + p1y * p1y) / r; + + if (ray_tmin <= t1 && t1 <= ray_tmax && v <= 1.0f) + { + t = t1; + x_hit = p1x; + y_hit = p1y; + valid = true; + } + } + else + { + float discr = B * B - 4.0f * A * C; + if (discr >= 0.0f) + { + float sqrt_discr = sqrtf(discr); + t1 = -0.5f * (B + sqrt_discr) / A; + t2 = -0.5f * (B - sqrt_discr) / A; + + const float p1x = ox + t1 * dx; + const float p1y = oy + t1 * dy; + const float p2x = ox + t2 * dx; + const float p2y = oy + t2 * dy; + + const float v1 = sqrtf(p1x * p1x + p1y * p1y) / r; + const float v2 = sqrtf(p2x * p2x + p2y * p2y) / r; + + if (ray_tmin <= t1 && t1 <= ray_tmax && v1 <= 1.0f) + { + t = t1; + x_hit = p1x; + y_hit = p1y; + valid = true; + } + else if (ray_tmin <= t2 && t2 <= ray_tmax && v2 <= 1.0f) + { + t = t2; + x_hit = p2x; + y_hit = p2y; + valid = true; + } + } + } + + if (!valid) + { + return; + } + + // + // Compute the surface normal at the hit on the paraboloid. + // The height function is: + // f(x,y) = (curv_x/2)*x^2 + (curv_y/2)*y^2 + // so its partial derivatives are: + // f_x = curv_x * x and f_y = curv_y * y. + // Then the (unnormalized) local normal is: + // N_local = (-f_x, -f_y, 1) = ( -curv_x*x_hit, -curv_y*y_hit, 1 ). + // + float3 N_local = normalize(make_float3(-cx * x_hit, + -cy * y_hit, + 1.0f)); + // Transform the normal back to world coordinates. + float3 world_normal = normalize(N_local.x * e1 + + N_local.y * e2 + + N_local.z * n); + + // Compute the hit point in world space. + float3 world_hit = ray_orig + t * ray_dir; + + // Report the intersection. + // Here, the two reported extra attributes are the parametric coordinates (a1, a2), + // encoded as unsigned integers. + optixReportIntersection(t, 0, + __float_as_uint(world_normal.x), + __float_as_uint(world_normal.y), + __float_as_uint(world_normal.z)); +} From b01d8ec57c3c8e8ab031585d23e686a35e1bb6d6 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:05:48 -0600 Subject: [PATCH 09/30] Add test for parabolic circle --- .../geometry_intersection_test.cpp | 66 +++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index 70dbf0e4..654e6e10 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -747,3 +747,69 @@ TEST(OptixRunner, FlatAnnulus_PartialArc) << ", false positives: " << fpos << ", false negatives: " << fneg << std::endl; } + +TEST(OptixRunner, ParabolicCircle) +{ + const double R = 5.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = 0.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(2 * R); + + // z(x,y) = x^2/(2*FOCAL_X) + y^2/(2*FOCAL_Y). Maximum over the circle x^2+y^2<=R^2 + // is R^2 / (2 * min(FOCAL_X, FOCAL_Y)), achieved along the axis with lower focal length. + const double Z_MAX_OFFSET = R * R / (2.0 * std::min(FOCAL_X, FOCAL_Y)); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + if (id == test_elid) + { + // z is curved: Z_ELEM <= z <= Z_ELEM + Z_MAX_OFFSET + EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; + EXPECT_LE(p1[2], Z_ELEM + Z_MAX_OFFSET + TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + ++hits; + if (!aper->is_in(p1[0], p1[1])) ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + ++misses; + if (aper->is_in(p1[0], p1[1])) ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} From a34425832a0dd40166046822118f38209a8203b8 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:27:41 -0600 Subject: [PATCH 10/30] Build fixes --- .../optix_runner/OptixCSP/src/shaders/GeometryDataST.h | 3 ++- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 10 +++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index 149745e1..e12bb77c 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -25,7 +25,8 @@ namespace OptixCSP QUADRILATERAL_FLAT = 6, CIRCLE_FLAT = 7, HEXAGON_FLAT = 8, - ANNULUS_FLAT = 9 + ANNULUS_FLAT = 9, + CIRCLE_PARABOLIC = 10 }; struct Parallelogram diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 58a82fdb..7ebcec5a 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -345,7 +345,7 @@ extern "C" __global__ void __intersection__rectangle_parabolic() float3 e1 = rect.v1 * L1; // recovers the original direction of edge 1, unit vector float3 e2 = rect.v2 * L2; // recovers the original direction of edge 2, unit ve // The flat (undeformed) rectangle's normal is: - float3 n = normalize(cross(e2, e1)); + float3 n = normalize(cross(e1, e2)); // // Transform ray into local coordinates. @@ -741,7 +741,7 @@ extern "C" __global__ void __intersection__annulus_flat() extern "C" __global__ void __intersection__circle_parabolic() { - const OptixCSP::GeometryDataST::Circle_Parabolic &circp = params.geoemetry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); + const OptixCSP::GeometryDataST::Circle_Parabolic &circp = params.geometry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); const float3 ray_orig = optixGetWorldRayOrigin(); const float3 ray_dir = optixGetWorldRayDirection(); @@ -757,7 +757,7 @@ extern "C" __global__ void __intersection__circle_parabolic() // TODO: Compute the local ray direction and ray origin const float3 x_ax = circp.x_axis; const float3 y_ax = circp.y_axis; - const float3 n = normalize(cross(e2, e1)); + const float3 n = normalize(cross(x_ax, y_ax)); float3 d = ray_orig - center; float ox = dot(d, x_ax); @@ -776,6 +776,10 @@ extern "C" __global__ void __intersection__circle_parabolic() float t1 = 0.0f; float t2 = 0.0f; + float x_hit = 0.0f; + float y_hit = 0.0f; + bool valid = false; + if (fabsf(A) < eps) { t1 = -C / B; From 36b47d3095511ba0de61a830fe29d62c75fd9b20 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:38:09 -0600 Subject: [PATCH 11/30] Fix a few more errors --- .../OptixCSP/src/core/pipeline_manager.cpp | 10 +++++----- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp index 7f6b40ca..2dec3964 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp @@ -46,11 +46,11 @@ const std::map IntersectionKernelMap = { {OpticalEntityType::CIRCLE_FLAT, "__intersection__circle_flat"}, {OpticalEntityType::HEXAGON_FLAT, "__intersection__hexagon_flat"}, {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}, - {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}, - {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, - {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, - {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, - {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; + {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}}; + // {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, + // {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, + // {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, + // {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; pipelineManager::pipelineManager(SoltraceState &state) : m_state(state) {} diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 7ebcec5a..abceb30a 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -847,8 +847,8 @@ extern "C" __global__ void __intersection__circle_parabolic() -cy * y_hit, 1.0f)); // Transform the normal back to world coordinates. - float3 world_normal = normalize(N_local.x * e1 + - N_local.y * e2 + + float3 world_normal = normalize(N_local.x * x_ax + + N_local.y * y_ax + N_local.z * n); // Compute the hit point in world space. From e6f7f4ed8e00183f0edef22b463dd3e4a6ab6679 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:41:59 -0600 Subject: [PATCH 12/30] Fix another dumb error --- .../optix_runner/OptixCSP/src/core/CspElement.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 38843054..e6b76515 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -322,6 +322,7 @@ GeometryDataST CspElement::toDeviceGeometryData() const float cx = (float)(m_surface->get_curvature_1()); float cy = (float)(m_surface->get_curvature_2()); GeometryDataST::Circle_Parabolic heliostat(o, v1, v2, cx, cy, r); + geometry_data.setCircle_Parabolic(heliostat); } } From bbe5867471fa66178868893b9b8b75b96b87497a Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:47:57 -0600 Subject: [PATCH 13/30] Actually fix the dumb error --- .../optix_runner/OptixCSP/src/shaders/Soltrace.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h index 5426b0e6..d93362d8 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h @@ -44,11 +44,11 @@ namespace OptixCSP{ CIRCLE_FLAT = 5, HEXAGON_FLAT = 6, ANNULUS_FLAT = 7, - CIRCLE_PARABOLIC = 8, - HEXAGON_PARABOLIC = 9, - TRIANGLE_PARABOLIC = 10, - ANNULUS_PARABOLIC = 11, - QUADRILATERAL_PARABOLIC = 12, + CIRCLE_PARABOLIC = 8 + // HEXAGON_PARABOLIC = 9, + // TRIANGLE_PARABOLIC = 10, + // ANNULUS_PARABOLIC = 11, + // QUADRILATERAL_PARABOLIC = 12, NUM_OPTICAL_ENTITY_TYPES }; From 866dacacbf5c73e231f7c052131c56c0f27a6f05 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Mon, 8 Jun 2026 16:49:11 -0600 Subject: [PATCH 14/30] So much stupidity... --- .../optix_runner/OptixCSP/src/shaders/Soltrace.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h index d93362d8..c35d03d6 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h @@ -44,7 +44,7 @@ namespace OptixCSP{ CIRCLE_FLAT = 5, HEXAGON_FLAT = 6, ANNULUS_FLAT = 7, - CIRCLE_PARABOLIC = 8 + CIRCLE_PARABOLIC = 8, // HEXAGON_PARABOLIC = 9, // TRIANGLE_PARABOLIC = 10, // ANNULUS_PARABOLIC = 11, From e1e2d92ed0acd6f4af019b91853d0a4413e5cb3a Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Tue, 9 Jun 2026 08:50:54 -0600 Subject: [PATCH 15/30] Testing fixes --- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 2 +- .../optix_runner/geometry_intersection_test.cpp | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index abceb30a..e065a254 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -345,7 +345,7 @@ extern "C" __global__ void __intersection__rectangle_parabolic() float3 e1 = rect.v1 * L1; // recovers the original direction of edge 1, unit vector float3 e2 = rect.v2 * L2; // recovers the original direction of edge 2, unit ve // The flat (undeformed) rectangle's normal is: - float3 n = normalize(cross(e1, e2)); + float3 n = normalize(cross(e2, e1)); // // Transform ray into local coordinates. diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index 654e6e10..45a0b74d 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -102,7 +102,7 @@ element_id set_default_sd(SimulationData &sd, stop->set_origin(0, 0, Z_BACKSTOP); stop->set_aim_vector(0, 0, 100); stop->set_surface(make_surface()); - stop->set_aperture(make_aperture(sx, sy)); + stop->set_aperture(make_aperture(2.0 * sx, 2.0 * sy)); sd.add_element(stop); // Set parameters @@ -759,7 +759,7 @@ TEST(OptixRunner, ParabolicCircle) // z(x,y) = x^2/(2*FOCAL_X) + y^2/(2*FOCAL_Y). Maximum over the circle x^2+y^2<=R^2 // is R^2 / (2 * min(FOCAL_X, FOCAL_Y)), achieved along the axis with lower focal length. - const double Z_MAX_OFFSET = R * R / (2.0 * std::min(FOCAL_X, FOCAL_Y)); + // const double Z_MAX_OFFSET = R * R / (2.0 * std::min(FOCAL_X, FOCAL_Y)); uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; @@ -793,8 +793,10 @@ TEST(OptixRunner, ParabolicCircle) if (id == test_elid) { // z is curved: Z_ELEM <= z <= Z_ELEM + Z_MAX_OFFSET - EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; - EXPECT_LE(p1[2], Z_ELEM + Z_MAX_OFFSET + TOL * Z_ELEM) << "ray " << i; + // EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; + // EXPECT_LE(p1[2], Z_ELEM + Z_MAX_OFFSET + TOL * Z_ELEM) << "ray " << i; + double zsol = p1[0] * p1[0] / (4.0 * FOCAL_X) + p1[1] * p1[1] / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM); EXPECT_TRUE(aper->is_in(p1[0], p1[1])); ++hits; if (!aper->is_in(p1[0], p1[1])) ++fpos; From 4b9f41f101a7f5416cf0434ce7c7539d58d0ae92 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Tue, 9 Jun 2026 10:58:31 -0600 Subject: [PATCH 16/30] Move to helpers for parabolic computations; add parabolic hexagon; add tests and stubs for other methods; expect many problems... --- .../OptixCSP/src/core/geometry_manager.cpp | 4 + .../OptixCSP/src/core/pipeline_manager.cpp | 10 +- .../OptixCSP/src/shaders/GeometryDataST.h | 82 ++- .../OptixCSP/src/shaders/Soltrace.h | 8 +- .../OptixCSP/src/shaders/intersection.cu | 563 ++++++++---------- .../OptixCSP/src/shaders/materials.cu | 17 +- .../geometry_intersection_test.cpp | 371 +++++++++++- 7 files changed, 699 insertions(+), 356 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index 6e07ac71..9bc088f5 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -84,6 +84,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::HEXAGON_FLAT); } + elseif (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::HEXAGON_PARABOLIC); + } else { std::stringstream ss; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp index 2dec3964..7f6b40ca 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp @@ -46,11 +46,11 @@ const std::map IntersectionKernelMap = { {OpticalEntityType::CIRCLE_FLAT, "__intersection__circle_flat"}, {OpticalEntityType::HEXAGON_FLAT, "__intersection__hexagon_flat"}, {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}, - {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}}; - // {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, - // {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, - // {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, - // {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; + {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}, + {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, + {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, + {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, + {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; pipelineManager::pipelineManager(SoltraceState &state) : m_state(state) {} diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index e12bb77c..31f22db7 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -26,7 +26,11 @@ namespace OptixCSP CIRCLE_FLAT = 7, HEXAGON_FLAT = 8, ANNULUS_FLAT = 9, - CIRCLE_PARABOLIC = 10 + CIRCLE_PARABOLIC = 10, + HEXAGON_PARABOLIC = 11, + TRIANGLE_PARABOLIC = 12, + ANNULUS_PARABOLIC = 13, + QUADRILATERAL_PARABOLIC = 14 }; struct Parallelogram @@ -208,6 +212,26 @@ namespace OptixCSP float radius; }; + struct Hexagon_Parabolic + { + Hexagon_Parabolic() = default; + Hexagon_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, const float &side_len) : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), s(side_len) + { + } + float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; + float s; + }; + + struct Triangle_Parabolic{}; + struct Annulus_Parabolic{}; + struct Quadrilateral_Parabolic{}; + GeometryDataST() = default; void setParallelogram(const Parallelogram &p) @@ -340,6 +364,58 @@ namespace OptixCSP return circle_parabolic; } + void setHexagon_Parabolic(const Hexagon_Parabolic &hexp) + { + assert(type == UNKNOWN_TYPE); + type = HEXAGON_PARABOLIC; + hexagon_parabolic = hexp; + } + + __host__ __device__ const Hexagon_Parabolic &getHexagon_Parabolic() const + { + assert(type == HEXAGON_PARABOLIC); + return hexagon_parabolic; + } + + void setTriangle_Parabolic(const Triangle_Parabolic &tp) + { + assert(type == UNKNOWN_TYPE); + type = TRIANGLE_PARABOLIC; + triangle_parabolic = tp; + } + + __host__ __device__ const Triangle_Parabolic &getTriangle_Parabolic() const + { + assert(type == TRIANGLE_PARABOLIC); + return triangle_parabolic; + } + + void setAnnulus_Parabolic(const Annulus_Parabolic &ap) + { + assert(type == UNKNOWN_TYPE); + type = ANNULUS_PARABOLIC; + annulus_parabolic = ap; + } + + __host__ __device__ const Annulus_Parabolic &getAnnulus_Parabolic() const + { + assert(type == ANNULUS_PARABOLIC); + return annulus_parabolic; + } + + void setQuadrilateral_Parabolic(const Quadrilateral_Parabolic &qp) + { + assert(type == UNKNOWN_TYPE); + type = QUADRILATERAL_PARABOLIC; + quadrilateral_parabolic = qp; + } + + __host__ __device__ const Quadrilateral_Parabolic &getQuadrilateral_Parabolic() const + { + assert(type == QUADRILATERAL_PARABOLIC); + return quadrilateral_parabolic; + } + Type type = UNKNOWN_TYPE; int32_t id = OptixCSP::kElementIdUnassigned; @@ -357,6 +433,10 @@ namespace OptixCSP Hexagon_Flat hexagon_flat; Annulus_Flat annulus_flat; Circle_Parabolic circle_parabolic; + Hexagon_Parabolic hexagon_parabolic; + Triangle_Parabolic triangle_parabolic; + Annulus_Parabolic annulus_parabolic; + Quadrilateral_Parabolic quadrilateral_parabolic; }; }; } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h index c35d03d6..5426b0e6 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h @@ -45,10 +45,10 @@ namespace OptixCSP{ HEXAGON_FLAT = 6, ANNULUS_FLAT = 7, CIRCLE_PARABOLIC = 8, - // HEXAGON_PARABOLIC = 9, - // TRIANGLE_PARABOLIC = 10, - // ANNULUS_PARABOLIC = 11, - // QUADRILATERAL_PARABOLIC = 12, + HEXAGON_PARABOLIC = 9, + TRIANGLE_PARABOLIC = 10, + ANNULUS_PARABOLIC = 11, + QUADRILATERAL_PARABOLIC = 12, NUM_OPTICAL_ENTITY_TYPES }; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index e065a254..bdcfd82c 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -303,175 +303,6 @@ extern "C" __global__ void __intersection__cylinder_y_capped() __float_as_uint(world_normal.z)); } -// For a parabolic surface rectangle aperture where -// the base (normal projection) is defined by the center and its two unit edge vectors -// In a local coordinate system (with origin at the anchor) the flat rectangle covers: -// x in [0, L1] and y in [0, L2], -// where L1 and L2 are the lengths of the original edge vectors. -// The parabolic surface is given by: -// z = (curv_x/2)*x^2 + (curv_y/2)*y^2 -// and the ray (in local coordinates) is: -// (ox,oy,oz) + t*(dx,dy,dz) -// We solve for t such that: -// oz + t*dz = (curv_x/2)*(ox+t*dx)^2 + (curv_y/2)*(oy+t*dy)^2 -// which expands into a quadratic: A*t^2 + B*t + C = 0. -// After finding the valid t, we compute the local hit (x,y) and then check that -// 0 <= x <= L1 and 0 <= y <= L2. -// Finally, we compute the surface normal from the paraboloid derivative -// f_x = curv_x * x and f_y = curv_y * y, -// so that the (unnormalized) local normal is (-f_x, -f_y, 1). -// -// The local hit point is then transformed back to world space for reporting. -extern "C" __global__ void __intersection__rectangle_parabolic() -{ - const OptixCSP::GeometryDataST::Rectangle_Parabolic &rect = params.geometry_data_array[optixGetPrimitiveIndex()].getRectangleParabolic(); - // Get ray information. - const float3 ray_orig = optixGetWorldRayOrigin(); - const float3 ray_dir = optixGetWorldRayDirection(); - const float ray_tmin = optixGetRayTmin(); - const float ray_tmax = optixGetRayTmax(); - - // - // Build the local coordinate system. - // - // We assume that the rectangle was defined with an anchor at its corner and - // two edge vectors. The stored rect.v1 and rect.v2 are the reciprocals: - // stored_v1 = original_v1 / dot(original_v1, original_v1) - // Thus, the original edge lengths are: - // Note rect.v1 had the size of 1/original_v1_length - float L1 = 1.0f / length(rect.v1); - float L2 = 1.0f / length(rect.v2); - // And the unit edge directions are: - float3 e1 = rect.v1 * L1; // recovers the original direction of edge 1, unit vector - float3 e2 = rect.v2 * L2; // recovers the original direction of edge 2, unit ve - // The flat (undeformed) rectangle's normal is: - float3 n = normalize(cross(e2, e1)); - - // - // Transform ray into local coordinates. - // The local coordinates (x,y,z) are defined such that: - // - The origin is at rect.anchor. - // - The x-axis is e1. - // - The y-axis is e2. - // - The z-axis is n. - // - // Compute the rectangle center (shifting from the lower-right corner) - float3 rect_center = rect.anchor + (L1 / 2.0f) * e1 + (L2 / 2.0f) * e2; - - float3 d = ray_orig - rect_center; - float ox = dot(d, e1); - float oy = dot(d, e2); - float oz = dot(d, n); - - float dx = dot(ray_dir, e1); - float dy = dot(ray_dir, e2); - float dz = dot(ray_dir, n); - - // Retrieve curvature parameters. - const float curv_x = rect.curv_x; - const float curv_y = rect.curv_y; - - float A = (curv_x * 0.5f) * (dx * dx) + (curv_y * 0.5f) * (dy * dy); - float B = curv_x * (ox * dx) + curv_y * (oy * dy) - dz; - float C = (curv_x * 0.5f) * (ox * ox) + (curv_y * 0.5f) * (oy * oy) - oz; - - float t = 0.0f; - float t1 = 0.0f; - float t2 = 0.0f; - float x_hit = 0.0f; - float y_hit = 0.0f; - const float eps = 1e-6f; - bool valid = false; - - if (fabsf(A) < eps) - { - if (fabsf(B) > eps) - { - t1 = -C / B; - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float a1 = p1x / (L1 / 2.0f); - const float a2 = p1y / (L2 / 2.0f); - - if ((t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax) && - !(a1 < -1.0f || a1 > 1.0f || a2 < -1.0f || a2 > 1.0f)) - { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; - } - } - } - else - { - float discr = B * B - 4.0f * A * C; - if (discr >= 0.0f) - { - float sqrt_discr = sqrtf(discr); - t1 = -0.5f * (B + sqrt_discr) / A; - t2 = -0.5f * (B - sqrt_discr) / A; - - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float p2x = ox + t2 * dx; - const float p2y = oy + t2 * dy; - const float a1_1 = p1x / (L1 / 2.0f); - const float a2_1 = p1y / (L2 / 2.0f); - const float a1_2 = p2x / (L1 / 2.0f); - const float a2_2 = p2y / (L2 / 2.0f); - - if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax && !(a1_1 < -1.0f || a1_1 > 1.0f || a2_1 < -1.0f || a2_1 > 1.0f)) - { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; - } - else if (t2 > 0.0f && t2 >= ray_tmin && t2 <= ray_tmax && !(a1_2 < -1.0f || a1_2 > 1.0f || a2_2 < -1.0f || a2_2 > 1.0f)) - { - t = t2; - x_hit = p2x; - y_hit = p2y; - valid = true; - } - } - } - - if (!valid) - { - return; - } - - // - // Compute the surface normal at the hit on the paraboloid. - // The height function is: - // f(x,y) = (curv_x/2)*x^2 + (curv_y/2)*y^2 - // so its partial derivatives are: - // f_x = curv_x * x and f_y = curv_y * y. - // Then the (unnormalized) local normal is: - // N_local = (-f_x, -f_y, 1) = ( -curv_x*x_hit, -curv_y*y_hit, 1 ). - // - float3 N_local = normalize(make_float3(-curv_x * x_hit, - -curv_y * y_hit, - 1.0f)); - // Transform the normal back to world coordinates. - float3 world_normal = normalize(N_local.x * e1 + - N_local.y * e2 + - N_local.z * n); - - // Compute the hit point in world space. - float3 world_hit = ray_orig + t * ray_dir; - - // Report the intersection. - // Here, the two reported extra attributes are the parametric coordinates (a1, a2), - // encoded as unsigned integers. - optixReportIntersection(t, 0, - __float_as_uint(world_normal.x), - __float_as_uint(world_normal.y), - __float_as_uint(world_normal.z)); -} - // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm extern "C" __device__ __inline__ float _triangle_intersect( @@ -509,51 +340,6 @@ extern "C" __device__ __inline__ float _triangle_intersect( return t; } -// // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) -// // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm -// extern "C" __global__ void __intersection__triangle_flat() -// { -// const OptixCSP::GeometryDataST::Triangle_Flat& tri = params.geometry_data_array[optixGetPrimitiveIndex()].getTriangle_Flat(); - -// const float3 ro = optixGetObjectRayOrigin(); -// const float3 rd = optixGetObjectRayDirection(); - -// //printf("Ray origin: (%f,%f,%f), direction: (%f,%f,%f)\n", ro.x, ro.y, ro.z, rd.x, rd.y, rd.z); - -// const float3 edge1 = tri.e1; -// const float3 edge2 = tri.e2; - -// const float3 pvec = cross(rd, edge2); -// const float det = dot(edge1, pvec); - -// // Backface culling + parallel rejection -// // (det must be strictly positive and not tiny) -// const float eps = 1e-8f; -// if (det <= eps) return; - -// const float inv_det = 1.0f / det; - -// const float3 tvec = ro - tri.v0; -// const float u = dot(tvec, pvec) * inv_det; -// if (u < 0.0f || u > 1.0f) return; - -// const float3 qvec = cross(tvec, edge1); -// const float v = dot(rd, qvec) * inv_det; -// if (v < 0.0f || (u + v) > 1.0f) -// return; - -// const float t = dot(edge2, qvec) * inv_det; -// if (t < optixGetRayTmin() || t > optixGetRayTmax()) return; - -// float3 world_normal = tri.normal; - -// optixReportIntersection(t, 0, -// __float_as_uint(world_normal.x), -// __float_as_uint(world_normal.y), -// __float_as_uint(world_normal.z)); - -// } - // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm extern "C" __global__ void __intersection__triangle_flat() @@ -719,7 +505,6 @@ extern "C" __global__ void __intersection__annulus_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { - // TODO: Need to adjust for possible rotation... float3 p = ray_orig + ray_dir * t - anf.center; float d = length(p); if (anf.ri <= d && d <= anf.ro) @@ -739,126 +524,284 @@ extern "C" __global__ void __intersection__annulus_flat() } } -extern "C" __global__ void __intersection__circle_parabolic() +// ----------------------------------------------------------------------- +// Shared helpers for parabolic surface intersections. +// +// All parabolic surfaces share the same quadric equation: +// z = (cx/2)*x^2 + (cy/2)*y^2 +// in a local frame (center, x_ax, y_ax, n=cross(x_ax,y_ax)). +// The three helpers below factor out the ray transform, quadratic solve, +// and normal computation. Each kernel only supplies the aperture test. +// ----------------------------------------------------------------------- + +// Transform a world-space ray into the local parabolic frame. +// Outputs the frame normal n = normalize(cross(x_ax, y_ax)) and +// the local ray origin (ox,oy,oz) and direction (dx,dy,dz). +extern "C" __device__ __inline__ void parabolic_ray_to_local( + const float3 &ray_orig, const float3 &ray_dir, + const float3 ¢er, + const float3 &x_ax, const float3 &y_ax, + float3 &n, + float &ox, float &oy, float &oz, + float &dx, float &dy, float &dz) { - const OptixCSP::GeometryDataST::Circle_Parabolic &circp = params.geometry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); + n = normalize(cross(x_ax, y_ax)); + const float3 d = ray_orig - center; + ox = dot(d, x_ax); + oy = dot(d, y_ax); + oz = dot(d, n); + dx = dot(ray_dir, x_ax); + dy = dot(ray_dir, y_ax); + dz = dot(ray_dir, n); +} - const float3 ray_orig = optixGetWorldRayOrigin(); - const float3 ray_dir = optixGetWorldRayDirection(); - const float ray_tmin = optixGetRayTmin(), ray_tmax = optixGetRayTmax(); +// Solve A*t^2 + B*t + C = 0 for the paraboloid-ray intersection and return +// up to two hits within [ray_tmin, ray_tmax], ordered by ascending t. +// Returns the number of valid hits (0, 1, or 2). +// t_out[i], lx_out[i], ly_out[i] give the ray parameter and local (x,y) of each hit. +extern "C" __device__ __inline__ int parabolic_solve( + float ox, float oy, float oz, + float dx, float dy, float dz, + float cx, float cy, + float ray_tmin, float ray_tmax, + float t_out[2], float lx_out[2], float ly_out[2]) +{ + const float A = 0.5f * cx * dx * dx + 0.5f * cy * dy * dy; + const float B = cx * ox * dx + cy * oy * dy - dz; + const float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; const float eps = 1e-6f; + int count = 0; - const float3 center = circp.center; - const float cx = circp.cx; - const float cy = circp.cy; - const float r = circp.radius; - - // TODO: Compute the local ray direction and ray origin - const float3 x_ax = circp.x_axis; - const float3 y_ax = circp.y_axis; - const float3 n = normalize(cross(x_ax, y_ax)); - - float3 d = ray_orig - center; - float ox = dot(d, x_ax); - float oy = dot(d, y_ax); - float oz = dot(d, n); + if (fabsf(A) < eps) + { + if (fabsf(B) > eps) + { + const float t = -C / B; + if (t >= ray_tmin && t <= ray_tmax) + { + t_out[0] = t; + lx_out[0] = ox + t * dx; + ly_out[0] = oy + t * dy; + count = 1; + } + } + } + else + { + const float discr = B * B - 4.0f * A * C; + if (discr >= 0.0f) + { + const float sq = sqrtf(discr); + // A > 0 (physical curvature), so ta <= tb is guaranteed. + const float ta = -0.5f * (B + sq) / A; + const float tb = -0.5f * (B - sq) / A; + if (ta >= ray_tmin && ta <= ray_tmax) + { + t_out[count] = ta; + lx_out[count] = ox + ta * dx; + ly_out[count] = oy + ta * dy; + ++count; + } + if (tb >= ray_tmin && tb <= ray_tmax) + { + t_out[count] = tb; + lx_out[count] = ox + tb * dx; + ly_out[count] = oy + tb * dy; + ++count; + } + } + } + return count; +} - float dx = dot(ray_dir, x_ax); - float dy = dot(ray_dir, y_ax); - float dz = dot(ray_dir, n); +// Compute the world-space unit normal at a parabolic surface hit. +// x_hit, y_hit : local coordinates of the hit point +// cx, cy : curvature parameters +// x_ax, y_ax : local frame unit vectors +// n : normalize(cross(x_ax, y_ax)) +extern "C" __device__ __inline__ float3 parabolic_world_normal( + float x_hit, float y_hit, + float cx, float cy, + const float3 &x_ax, const float3 &y_ax, const float3 &n) +{ + const float3 N_local = make_float3(-cx * x_hit, -cy * y_hit, 1.0f); + return N_local.x * x_ax + N_local.y * y_ax + N_local.z * n; +} - float A = (0.5f * cx) * (dx * dx) + (0.5f * cy) * (dy * dy); - float B = cx * (ox * dx) + cy * (oy * dy) - dz; - float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; +// ----------------------------------------------------------------------- - float t = 0.0f; - float t1 = 0.0f; - float t2 = 0.0f; +// Parabolic surface, rectangle aperture. +// z = (curv_x/2)*x^2 + (curv_y/2)*y^2, aperture: |x| <= L1/2, |y| <= L2/2. +extern "C" __global__ void __intersection__rectangle_parabolic() +{ + const OptixCSP::GeometryDataST::Rectangle_Parabolic &rect = + params.geometry_data_array[optixGetPrimitiveIndex()].getRectangleParabolic(); - float x_hit = 0.0f; - float y_hit = 0.0f; - bool valid = false; + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); - if (fabsf(A) < eps) + // Recover unit edge vectors and half-lengths from the stored reciprocal vectors. + // rect.v1 = original_v1 / dot(original_v1, original_v1), so |rect.v1| = 1/L1. + const float L1 = 1.0f / length(rect.v1); + const float L2 = 1.0f / length(rect.v2); + const float3 e1 = rect.v1 * L1; + const float3 e2 = rect.v2 * L2; + const float3 center = rect.anchor + (0.5f * L1) * e1 + (0.5f * L2) * e2; + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, center, e1, e2, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + rect.curv_x, rect.curv_y, + ray_tmin, ray_tmax, ts, lxs, lys); + + const float half_L1 = 0.5f * L1; + const float half_L2 = 0.5f * L2; + for (int i = 0; i < nc; ++i) { - t1 = -C / B; - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float v = sqrtf(p1x * p1x + p1y * p1y) / r; - - if (ray_tmin <= t1 && t1 <= ray_tmax && v <= 1.0f) + if (lxs[i] >= -half_L1 && lxs[i] <= half_L1 && + lys[i] >= -half_L2 && lys[i] <= half_L2) { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + rect.curv_x, rect.curv_y, + e1, e2, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; } } - else +} + +// Parabolic surface, circle aperture. +// z = (cx/2)*x^2 + (cy/2)*y^2, aperture: x^2 + y^2 <= radius^2. +extern "C" __global__ void __intersection__circle_parabolic() +{ + const OptixCSP::GeometryDataST::Circle_Parabolic &circp = + params.geometry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, circp.center, circp.x_axis, circp.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + circp.cx, circp.cy, + ray_tmin, ray_tmax, ts, lxs, lys); + + const float r2 = circp.radius * circp.radius; + for (int i = 0; i < nc; ++i) { - float discr = B * B - 4.0f * A * C; - if (discr >= 0.0f) + if (lxs[i] * lxs[i] + lys[i] * lys[i] <= r2) { - float sqrt_discr = sqrtf(discr); - t1 = -0.5f * (B + sqrt_discr) / A; - t2 = -0.5f * (B - sqrt_discr) / A; + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + circp.cx, circp.cy, + circp.x_axis, circp.y_axis, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } +} + +extern "C" __global__ void __intersection__hexagon_parabolic() +{ + const OptixCSP::GeometryDataST::Hexagon_Parabolic &hexp = + params.geometry_data_array[optixGetPrimitiveIndex()].getHexagon_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float p2x = ox + t2 * dx; - const float p2y = oy + t2 * dy; + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + hexp.center, hexp.x_axis, hexp.y_axis, + n, ox, oy, oz, dx, dy, dz); - const float v1 = sqrtf(p1x * p1x + p1y * p1y) / r; - const float v2 = sqrtf(p2x * p2x + p2y * p2y) / r; + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + hexp.cx, hexp.cy, + ray_tmin, ray_tmax, + ts, lxs, lys;) - if (ray_tmin <= t1 && t1 <= ray_tmax && v1 <= 1.0f) + for (int i = 0; i < nc; ++i) + { + bool is_in = false; + float3 p = ray_orig + ray_dir * ts[i] - hexp.center; + // Project onto the local x and y axes which are unit vectors + const float px = dot(p, hexp.x_axis); + const float py = dot(p, hexp.y_axis); + const float s = hexp.s; + const float xl = 0.5f * s; + const float yl = 0.5f * sqrtf(3.0f) * s; + if (-xl <= px && px <= xl && -yl <= py && py <= yl) + { + // Center + is_in = true; + } + else if (-s <= px && px < -xl) + { + // Left side + float y1 = sqrtf(3.0f) * (px + s); + float y2 = -y1; + if (y2 <= py && py <= y1) { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; + is_in = true; } - else if (ray_tmin <= t2 && t2 <= ray_tmax && v2 <= 1.0f) + } + else if (xl < px && px <= s) + { + // Right side + float y1 = sqrtf(3.0f) * (px - s); + float y2 = -y1; + if (y1 <= py && py <= y2) { - t = t2; - x_hit = p2x; - y_hit = p2y; - valid = true; + is_in = true; } } - } - if (!valid) - { - return; + if (is_in) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + hexp.curv_x, hexp.curv_y, + hexp.x_axis, hexp.y_axis, n); + + optixReportIntersection(t, + 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + + return; + } } +} - // - // Compute the surface normal at the hit on the paraboloid. - // The height function is: - // f(x,y) = (curv_x/2)*x^2 + (curv_y/2)*y^2 - // so its partial derivatives are: - // f_x = curv_x * x and f_y = curv_y * y. - // Then the (unnormalized) local normal is: - // N_local = (-f_x, -f_y, 1) = ( -curv_x*x_hit, -curv_y*y_hit, 1 ). - // - float3 N_local = normalize(make_float3(-cx * x_hit, - -cy * y_hit, - 1.0f)); - // Transform the normal back to world coordinates. - float3 world_normal = normalize(N_local.x * x_ax + - N_local.y * y_ax + - N_local.z * n); - - // Compute the hit point in world space. - float3 world_hit = ray_orig + t * ray_dir; - - // Report the intersection. - // Here, the two reported extra attributes are the parametric coordinates (a1, a2), - // encoded as unsigned integers. - optixReportIntersection(t, 0, - __float_as_uint(world_normal.x), - __float_as_uint(world_normal.y), - __float_as_uint(world_normal.z)); +extern "C" __global__ void __intersection__triangle_parabolic() +{ +} + +extern "C" __global__ void __intersection__annulus_parabolic() +{ +} + +extern "C" __global__ void __intersection__quadrilateral_parabolic() +{ } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu index 004dbf73..2b73f562 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu @@ -58,7 +58,7 @@ extern "C" __device__ __inline__ float3 orthonormal_vector(float3 v) // w = make_float3(0.0f, 1.0f, 0.0f); // u = cross(v, w); u = make_float3(-v.z, 0.0f, v.x); - } + } return normalize(u); } @@ -104,10 +104,13 @@ extern "C" __global__ void __closesthit__element() const bool optical_errors = params.optical_errors; // Fetch the normal vector from the hit attributes passed by OptiX - float3 object_normal = make_float3(__uint_as_float(optixGetAttribute_0()), __uint_as_float(optixGetAttribute_1()), + float3 object_normal = make_float3(__uint_as_float(optixGetAttribute_0()), + __uint_as_float(optixGetAttribute_1()), __uint_as_float(optixGetAttribute_2())); - // Transform the object-space normal to world space using OptiX built-in function - float3 world_normal = normalize(optixTransformNormalFromObjectToWorldSpace(object_normal)); + // // Transform the object-space normal to world space using OptiX built-in function + // float3 world_normal = normalize(optixTransformNormalFromObjectToWorldSpace(object_normal)); + // All reports from the intersection functions are in world coordinates + float3 world_normal = normalize(object_normal); // Compute the facing normal, which handles the direction of the normal based on the incoming ray direction const float3 ray_dir = optixGetWorldRayDirection(); @@ -140,8 +143,8 @@ extern "C" __global__ void __closesthit__element() const float transmissivity = material.transmissivity; const bool use_transmissivity = material.use_refraction; const float reflectivity = material.reflectivity; - const float normal_sigma = 1e-3f * material.slope_error; // Convert mrad to rad - const float spec_sigma = 1e-3f * material.specularity_error; // Convert mrad to rad + const float normal_sigma = 1e-3f * material.slope_error; // Convert mrad to rad + const float spec_sigma = 1e-3f * material.specularity_error; // Convert mrad to rad const uint8_t error_type = material.optical_dist; // Surface normal (macro-surface) errors @@ -288,7 +291,7 @@ extern "C" __global__ void __closesthit__element() else if (!absorbed) { // Ray hit an element but max depth was reached; count it (absorption at this depth does not count). - atomicAdd(reinterpret_cast(params.d_depth_exceeded_count), 1ULL); + atomicAdd(reinterpret_cast(params.d_depth_exceeded_count), 1ULL); } setPayload(prd); diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index 45a0b74d..489da846 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -20,7 +20,7 @@ const uint_fast64_t NRAYS = 10000; static void compute_global_bounding_box( const glm::dvec3 &lower_local, const glm::dvec3 &upper_local, - double zrot_deg, // z-axis rotation in degrees (local-to-global) + double zrot_deg, // z-axis rotation in degrees (local-to-global) const glm::dvec3 &origin, glm::dvec3 &lower_global, glm::dvec3 &upper_global) @@ -164,7 +164,8 @@ TEST(OptixRunner, FlatRectangle) // And that we are in the aperture EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { @@ -173,7 +174,8 @@ TEST(OptixRunner, FlatRectangle) // And that we are not in the aperture. EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -228,14 +230,16 @@ TEST(OptixRunner, FlatEquilateralTriangle) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -291,14 +295,16 @@ TEST(OptixRunner, FlatTriangle) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -356,14 +362,16 @@ TEST(OptixRunner, FlatQuadrilateral) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -423,14 +431,16 @@ TEST(OptixRunner, ParabolaRectangle) EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -487,14 +497,16 @@ TEST(OptixRunner, Cylinder) EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -545,14 +557,16 @@ TEST(OptixRunner, FlatCircle) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); ++hits; - if (!aper->is_in(p1[0], p1[1])) ++fpos; + if (!aper->is_in(p1[0], p1[1])) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); ++misses; - if (aper->is_in(p1[0], p1[1])) ++fneg; + if (aper->is_in(p1[0], p1[1])) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -607,14 +621,16 @@ TEST(OptixRunner, FlatHexagon) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -667,14 +683,16 @@ TEST(OptixRunner, FlatAnnulus_FullArc) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); ++hits; - if (!aper->is_in(p1[0], p1[1])) ++fpos; + if (!aper->is_in(p1[0], p1[1])) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); ++misses; - if (aper->is_in(p1[0], p1[1])) ++fneg; + if (aper->is_in(p1[0], p1[1])) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -731,14 +749,300 @@ TEST(OptixRunner, FlatAnnulus_PartialArc) EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(lx, ly)) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(lx, ly)) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicHexagon) +{ + const double S = 5.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = 30.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(2 * S); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicTriangle) +{ + const double d = 4.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = 90.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(d); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicAnnulus) +{ + const double R0 = 2.0; + const double R1 = 5.0; + const double ARC = 180.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = -15.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(R0, R1, ARC); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicQuadrilateral) +{ + // Parallelogram + const double x1 = 0.0, x2 = 3.0, x3 = (x2 - x1) + 1.0, x4 = x3 - x2 + x1; + const double y1 = 0.0, y2 = y1, y3 = 2.0, y4 = y3; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = -45.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(x1, y1, x2, y2, x3, y3, x4, y4); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); @@ -753,7 +1057,7 @@ TEST(OptixRunner, ParabolicCircle) const double R = 5.0; const double FOCAL_X = 10.0; const double FOCAL_Y = 30.0; - const double ROT_DEG = 0.0; + const double ROT_DEG = -45.0; auto surf = make_surface(FOCAL_X, FOCAL_Y); auto aper = make_aperture(2 * R); @@ -779,6 +1083,10 @@ TEST(OptixRunner, ParabolicCircle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -790,23 +1098,28 @@ TEST(OptixRunner, ParabolicCircle) EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + if (id == test_elid) { // z is curved: Z_ELEM <= z <= Z_ELEM + Z_MAX_OFFSET - // EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; + // EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; // EXPECT_LE(p1[2], Z_ELEM + Z_MAX_OFFSET + TOL * Z_ELEM) << "ray " << i; - double zsol = p1[0] * p1[0] / (4.0 * FOCAL_X) + p1[1] * p1[1] / (4.0 * FOCAL_Y) + Z_ELEM; - EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM); - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM); + EXPECT_TRUE(aper->is_in(lx, ly)); ++hits; - if (!aper->is_in(p1[0], p1[1])) ++fpos; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); ++misses; - if (aper->is_in(p1[0], p1[1])) ++fneg; + if (aper->is_in(lx, ly)) + ++fneg; } } EXPECT_GT(hits, 0u); From ce9122d6e9a414938b9f425f2bc8ef0465a411b7 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Tue, 9 Jun 2026 11:35:02 -0600 Subject: [PATCH 17/30] Build fixes --- .../optix_runner/OptixCSP/src/core/CspElement.cpp | 8 ++++++++ .../optix_runner/OptixCSP/src/core/geometry_manager.cpp | 2 +- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 6 +++--- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index e6b76515..b1aec8d4 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -341,6 +341,14 @@ GeometryDataST CspElement::toDeviceGeometryData() const GeometryDataST::Hexagon_Flat hex(o, n, vx, vy, s); geometry_data.setHexagon_Flat(hex); } + + if (surface_type == SurfaceType::PARABOLIC) + { + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Hexagon_Parabolic hex(o, vx, vy, cx, cy, s); + geometry_data.setHexagon_Parabolic(hex); + } } if (aperture_type == ApertureType::ANNULUS) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index 9bc088f5..92c86d92 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -84,7 +84,7 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::HEXAGON_FLAT); } - elseif (element->get_surface_type() == SurfaceType::PARABOLIC) + else if (element->get_surface_type() == SurfaceType::PARABOLIC) { sbt_offset = static_cast(OpticalEntityType::HEXAGON_PARABOLIC); } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index bdcfd82c..4843f095 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -739,7 +739,7 @@ extern "C" __global__ void __intersection__hexagon_parabolic() const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, hexp.cx, hexp.cy, ray_tmin, ray_tmax, - ts, lxs, lys;) + ts, lxs, lys); for (int i = 0; i < nc; ++i) { @@ -780,10 +780,10 @@ extern "C" __global__ void __intersection__hexagon_parabolic() if (is_in) { const float3 wn = parabolic_world_normal(lxs[i], lys[i], - hexp.curv_x, hexp.curv_y, + hexp.cx, hexp.cy, hexp.x_axis, hexp.y_axis, n); - optixReportIntersection(t, + optixReportIntersection(ts[i], 0, __float_as_uint(wn.x), __float_as_uint(wn.y), From 075dfaef5acfcf736da40ae106ac11e16ff2c84a Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Tue, 9 Jun 2026 12:01:46 -0600 Subject: [PATCH 18/30] Make parabolic rectangle more consistent with other rectangle aperture handling --- .../optix_runner/OptixCSP/src/core/CspElement.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index b1aec8d4..358210d0 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -224,10 +224,11 @@ GeometryDataST CspElement::toDeviceGeometryData() const if (surface_type == SurfaceType::PARABOLIC) { - Vec3d edge_x = v1 * (float)(-width); - Vec3d edge_y = v2 * (float)height; + Vec3d edge_x = v1 * (float)(width); + Vec3d edge_y = v2 * (float)(height); - Vec3d local_anchor(x_coord + width, y_coord, 0.0); + // Lower left corner + Vec3d local_anchor(x_coord, y_coord, 0.0); // float3 anchor = OptixCSP::toFloat3(m_origin - v1 * 0.5 - v2 * 0.5); Vec3d global_anchor = rotation_matrix * local_anchor + m_origin; From 6008c03da990694d47afc4ebc99ae28994daf0ee Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Wed, 10 Jun 2026 08:15:26 -0600 Subject: [PATCH 19/30] Commonize hexagon aperture check in optix runner --- .../OptixCSP/src/shaders/intersection.cu | 326 +++++++++--------- 1 file changed, 160 insertions(+), 166 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 4843f095..541579d6 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -9,12 +9,136 @@ extern "C" __constant__ OptixCSP::LaunchParams params; } +/**************** Surface Helper Functions ****************/ + +// ----------------------------------------------------------------------- +// Shared helpers for a planar (flat) surface +// +// All planar surfaces have the equation +// = 0 +// where n is the normal vector, p0 = (x0, y0, z0) is a point in the plane, +// generally the origin of the aperture. +// ----------------------------------------------------------------------- + extern "C" __device__ __inline__ float ray_distance_to_plane(float3 ro, float3 rd, float4 plane) { const float3 n = make_float3(plane); return (plane.w - dot(n, ro)) / dot(rd, n); } +// ----------------------------------------------------------------------- +// Shared helpers for parabolic surface intersections. +// +// All parabolic surfaces share the same quadric equation: +// z = (cx/2)*x^2 + (cy/2)*y^2 +// in a local frame (center, x_ax, y_ax, n=cross(x_ax,y_ax)). +// The three helpers below factor out the ray transform, quadratic solve, +// and normal computation. Each kernel only supplies the aperture test. +// ----------------------------------------------------------------------- + +// Transform a world-space ray into the local parabolic frame. +// Outputs the frame normal n = normalize(cross(x_ax, y_ax)) and +// the local ray origin (ox,oy,oz) and direction (dx,dy,dz). +extern "C" __device__ __inline__ void parabolic_ray_to_local( + const float3 &ray_orig, const float3 &ray_dir, + const float3 ¢er, + const float3 &x_ax, const float3 &y_ax, + float3 &n, + float &ox, float &oy, float &oz, + float &dx, float &dy, float &dz) +{ + n = normalize(cross(x_ax, y_ax)); + const float3 d = ray_orig - center; + ox = dot(d, x_ax); + oy = dot(d, y_ax); + oz = dot(d, n); + dx = dot(ray_dir, x_ax); + dy = dot(ray_dir, y_ax); + dz = dot(ray_dir, n); +} + +// Solve A*t^2 + B*t + C = 0 for the paraboloid-ray intersection and return +// up to two hits within [ray_tmin, ray_tmax], ordered by ascending t. +// Returns the number of valid hits (0, 1, or 2). +// t_out[i], lx_out[i], ly_out[i] give the ray parameter and local (x,y) of each hit. +extern "C" __device__ __inline__ int parabolic_solve( + float ox, float oy, float oz, + float dx, float dy, float dz, + float cx, float cy, + float ray_tmin, float ray_tmax, + float t_out[2], float lx_out[2], float ly_out[2]) +{ + const float A = 0.5f * cx * dx * dx + 0.5f * cy * dy * dy; + const float B = cx * ox * dx + cy * oy * dy - dz; + const float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; + + const float eps = 1e-6f; + int count = 0; + + if (fabsf(A) < eps) + { + if (fabsf(B) > eps) + { + const float t = -C / B; + if (t >= ray_tmin && t <= ray_tmax) + { + t_out[0] = t; + lx_out[0] = ox + t * dx; + ly_out[0] = oy + t * dy; + count = 1; + } + } + } + else + { + const float discr = B * B - 4.0f * A * C; + if (discr >= 0.0f) + { + const float sq = sqrtf(discr); + // A > 0 (physical curvature), so ta <= tb is guaranteed. + const float ta = -0.5f * (B + sq) / A; + const float tb = -0.5f * (B - sq) / A; + if (ta >= ray_tmin && ta <= ray_tmax) + { + t_out[count] = ta; + lx_out[count] = ox + ta * dx; + ly_out[count] = oy + ta * dy; + ++count; + } + if (tb >= ray_tmin && tb <= ray_tmax) + { + t_out[count] = tb; + lx_out[count] = ox + tb * dx; + ly_out[count] = oy + tb * dy; + ++count; + } + } + } + return count; +} + +// Compute the world-space unit normal at a parabolic surface hit. +// x_hit, y_hit : local coordinates of the hit point +// cx, cy : curvature parameters +// x_ax, y_ax : local frame unit vectors +// n : normalize(cross(x_ax, y_ax)) +extern "C" __device__ __inline__ float3 parabolic_world_normal( + float x_hit, float y_hit, + float cx, float cy, + const float3 &x_ax, const float3 &y_ax, const float3 &n) +{ + const float3 N_local = make_float3(-cx * x_hit, -cy * y_hit, 1.0f); + return N_local.x * x_ax + N_local.y * y_ax + N_local.z * n; +} + +// ----------------------------------------------------------------------- + + +/**************** Aperture Helper Functions ****************/ + + +/**************** Optix Intersection Functions ****************/ + extern "C" __global__ void __intersection__parallelogram() { int i = optixGetPrimitiveIndex(); @@ -430,6 +554,39 @@ extern "C" __global__ void __intersection__circle_flat() } } +extern "C" __device__ __inline__ bool hexagon_contains(float px, float py, float s) +{ + bool is_in = false; + const float xl = 0.5f * s; + const float yl = 0.5f * sqrtf(3.0f) * s; + if (-xl <= px && px <= xl && -yl <= py && py <= yl) + { + // Center + is_in = true; + } + else if (-s <= px && px < -xl) + { + // Left side + float y1 = sqrtf(3.0f) * (px + s); + float y2 = -y1; + if (y2 <= py && py <= y1) + { + is_in = true; + } + } + else if (xl < px && px <= s) + { + // Right side + float y1 = sqrtf(3.0f) * (px - s); + float y2 = -y1; + if (y1 <= py && py <= y2) + { + is_in = true; + } + } + return is_in; +} + extern "C" __global__ void __intersection__hexagon_flat() { const OptixCSP::GeometryDataST::Hexagon_Flat &hex = params.geometry_data_array[optixGetPrimitiveIndex()].getHexagon_Flat(); @@ -445,41 +602,13 @@ extern "C" __global__ void __intersection__hexagon_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { - bool is_in = false; float3 p = ray_orig + ray_dir * t - hex.center; // Project onto the local x and y axes which are unit vectors const float px = dot(p, hex.x_axis); const float py = dot(p, hex.y_axis); const float s = hex.s; - const float xl = 0.5f * s; - const float yl = 0.5f * sqrtf(3.0f) * s; - if (-xl <= px && px <= xl && -yl <= py && py <= yl) - { - // Center - is_in = true; - } - else if (-s <= px && px < -xl) - { - // Left side - float y1 = sqrtf(3.0f) * (px + s); - float y2 = -y1; - if (y2 <= py && py <= y1) - { - is_in = true; - } - } - else if (xl < px && px <= s) - { - // Right side - float y1 = sqrtf(3.0f) * (px - s); - float y2 = -y1; - if (y1 <= py && py <= y2) - { - is_in = true; - } - } - if (is_in) + if (hexagon_contains(px, py, s)) { optixReportIntersection(t, 0, @@ -524,113 +653,6 @@ extern "C" __global__ void __intersection__annulus_flat() } } -// ----------------------------------------------------------------------- -// Shared helpers for parabolic surface intersections. -// -// All parabolic surfaces share the same quadric equation: -// z = (cx/2)*x^2 + (cy/2)*y^2 -// in a local frame (center, x_ax, y_ax, n=cross(x_ax,y_ax)). -// The three helpers below factor out the ray transform, quadratic solve, -// and normal computation. Each kernel only supplies the aperture test. -// ----------------------------------------------------------------------- - -// Transform a world-space ray into the local parabolic frame. -// Outputs the frame normal n = normalize(cross(x_ax, y_ax)) and -// the local ray origin (ox,oy,oz) and direction (dx,dy,dz). -extern "C" __device__ __inline__ void parabolic_ray_to_local( - const float3 &ray_orig, const float3 &ray_dir, - const float3 ¢er, - const float3 &x_ax, const float3 &y_ax, - float3 &n, - float &ox, float &oy, float &oz, - float &dx, float &dy, float &dz) -{ - n = normalize(cross(x_ax, y_ax)); - const float3 d = ray_orig - center; - ox = dot(d, x_ax); - oy = dot(d, y_ax); - oz = dot(d, n); - dx = dot(ray_dir, x_ax); - dy = dot(ray_dir, y_ax); - dz = dot(ray_dir, n); -} - -// Solve A*t^2 + B*t + C = 0 for the paraboloid-ray intersection and return -// up to two hits within [ray_tmin, ray_tmax], ordered by ascending t. -// Returns the number of valid hits (0, 1, or 2). -// t_out[i], lx_out[i], ly_out[i] give the ray parameter and local (x,y) of each hit. -extern "C" __device__ __inline__ int parabolic_solve( - float ox, float oy, float oz, - float dx, float dy, float dz, - float cx, float cy, - float ray_tmin, float ray_tmax, - float t_out[2], float lx_out[2], float ly_out[2]) -{ - const float A = 0.5f * cx * dx * dx + 0.5f * cy * dy * dy; - const float B = cx * ox * dx + cy * oy * dy - dz; - const float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; - - const float eps = 1e-6f; - int count = 0; - - if (fabsf(A) < eps) - { - if (fabsf(B) > eps) - { - const float t = -C / B; - if (t >= ray_tmin && t <= ray_tmax) - { - t_out[0] = t; - lx_out[0] = ox + t * dx; - ly_out[0] = oy + t * dy; - count = 1; - } - } - } - else - { - const float discr = B * B - 4.0f * A * C; - if (discr >= 0.0f) - { - const float sq = sqrtf(discr); - // A > 0 (physical curvature), so ta <= tb is guaranteed. - const float ta = -0.5f * (B + sq) / A; - const float tb = -0.5f * (B - sq) / A; - if (ta >= ray_tmin && ta <= ray_tmax) - { - t_out[count] = ta; - lx_out[count] = ox + ta * dx; - ly_out[count] = oy + ta * dy; - ++count; - } - if (tb >= ray_tmin && tb <= ray_tmax) - { - t_out[count] = tb; - lx_out[count] = ox + tb * dx; - ly_out[count] = oy + tb * dy; - ++count; - } - } - } - return count; -} - -// Compute the world-space unit normal at a parabolic surface hit. -// x_hit, y_hit : local coordinates of the hit point -// cx, cy : curvature parameters -// x_ax, y_ax : local frame unit vectors -// n : normalize(cross(x_ax, y_ax)) -extern "C" __device__ __inline__ float3 parabolic_world_normal( - float x_hit, float y_hit, - float cx, float cy, - const float3 &x_ax, const float3 &y_ax, const float3 &n) -{ - const float3 N_local = make_float3(-cx * x_hit, -cy * y_hit, 1.0f); - return N_local.x * x_ax + N_local.y * y_ax + N_local.z * n; -} - -// ----------------------------------------------------------------------- - // Parabolic surface, rectangle aperture. // z = (curv_x/2)*x^2 + (curv_y/2)*y^2, aperture: |x| <= L1/2, |y| <= L2/2. extern "C" __global__ void __intersection__rectangle_parabolic() @@ -741,43 +763,15 @@ extern "C" __global__ void __intersection__hexagon_parabolic() ray_tmin, ray_tmax, ts, lxs, lys); - for (int i = 0; i < nc; ++i) + for (int i = 0; i < nc; ++i) { - bool is_in = false; float3 p = ray_orig + ray_dir * ts[i] - hexp.center; // Project onto the local x and y axes which are unit vectors const float px = dot(p, hexp.x_axis); const float py = dot(p, hexp.y_axis); const float s = hexp.s; - const float xl = 0.5f * s; - const float yl = 0.5f * sqrtf(3.0f) * s; - if (-xl <= px && px <= xl && -yl <= py && py <= yl) - { - // Center - is_in = true; - } - else if (-s <= px && px < -xl) - { - // Left side - float y1 = sqrtf(3.0f) * (px + s); - float y2 = -y1; - if (y2 <= py && py <= y1) - { - is_in = true; - } - } - else if (xl < px && px <= s) - { - // Right side - float y1 = sqrtf(3.0f) * (px - s); - float y2 = -y1; - if (y1 <= py && py <= y2) - { - is_in = true; - } - } - if (is_in) + if (hexagon_contains(px, py, s)) { const float3 wn = parabolic_world_normal(lxs[i], lys[i], hexp.cx, hexp.cy, From 61965f039d2d6c31778caa4488d39052a484043f Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Wed, 10 Jun 2026 09:01:09 -0600 Subject: [PATCH 20/30] Initital implementation of parabolic annulus for optix runner --- .../OptixCSP/src/core/CspElement.cpp | 29 ++-- .../OptixCSP/src/core/geometry_manager.cpp | 4 + .../OptixCSP/src/shaders/GeometryDataST.h | 36 ++++- .../OptixCSP/src/shaders/intersection.cu | 136 +++++++++++++----- 4 files changed, 155 insertions(+), 50 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 358210d0..eb1ebe28 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -308,13 +308,13 @@ GeometryDataST CspElement::toDeviceGeometryData() const float r = circ.get_radius(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); - + if (surface_type == SurfaceType::FLAT) { GeometryDataST::Circle_Flat heliostat(o, n, r); geometry_data.setCircle_Flat(heliostat); } - + if (surface_type == SurfaceType::PARABOLIC) { Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix @@ -337,19 +337,20 @@ GeometryDataST CspElement::toDeviceGeometryData() const float s = hex.get_side_length(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { GeometryDataST::Hexagon_Flat hex(o, n, vx, vy, s); geometry_data.setHexagon_Flat(hex); } - if (surface_type == SurfaceType::PARABOLIC) - { - float cx = (float)(m_surface->get_curvature_1()); - float cy = (float)(m_surface->get_curvature_2()); - GeometryDataST::Hexagon_Parabolic hex(o, vx, vy, cx, cy, s); - geometry_data.setHexagon_Parabolic(hex); - } + if (surface_type == SurfaceType::PARABOLIC) + { + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Hexagon_Parabolic hex(o, vx, vy, cx, cy, s); + geometry_data.setHexagon_Parabolic(hex); + } } if (aperture_type == ApertureType::ANNULUS) @@ -364,12 +365,22 @@ GeometryDataST CspElement::toDeviceGeometryData() const float arc = anf.get_arc(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { GeometryDataST::Annulus_Flat anf(o, n, vx, vy, radius_in, radius_out, arc); geometry_data.setAnnulus_Flat(anf); } + + if (surface_type == SurfaceType::PARABOLIC) + { + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Annulus_Parabolic anp(o, vx, vy, cx, cy, + radius_in, radius_out, arc); + geometry_data.setAnnulus_Parabolic(anp); + } } geometry_data.id = this->m_id; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index 92c86d92..ca95ef29 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -46,6 +46,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::ANNULUS_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::ANNULUS_PARABOLIC); + } else { std::stringstream ss; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index 31f22db7..d9c090fd 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -216,8 +216,9 @@ namespace OptixCSP { Hexagon_Parabolic() = default; Hexagon_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, - const float &curv_x, const float &curv_y, const float &side_len) : center(origin), x_axis(x_ax), y_axis(y_ax), - cx(curv_x), cy(curv_y), s(side_len) + const float &curv_x, const float &curv_y, const float &side_len) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), s(side_len) { } float3 center; @@ -228,9 +229,34 @@ namespace OptixCSP float s; }; - struct Triangle_Parabolic{}; - struct Annulus_Parabolic{}; - struct Quadrilateral_Parabolic{}; + struct Triangle_Parabolic + { + }; + + struct Annulus_Parabolic + { + Annulus_Parabolic() = default; + Annulus_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float &r_inner, const float &r_outer, const float &arc) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), + ri(r_inner), ro(r_outer), arc(arc) + { + } + float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; + float ri; + float ro; + float arc; // in radians + }; + + struct Quadrilateral_Parabolic + { + }; GeometryDataST() = default; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 541579d6..7de9dd04 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -20,6 +20,9 @@ extern "C" // generally the origin of the aperture. // ----------------------------------------------------------------------- +// Return the ray parameter t at which the ray (ro + t*rd) intersects the plane. +// The plane is encoded as float4(nx, ny, nz, d) where nx/ny/nz is the unit normal +// and d = dot(n, p0) for any point p0 on the plane. extern "C" __device__ __inline__ float ray_distance_to_plane(float3 ro, float3 rd, float4 plane) { const float3 n = make_float3(plane); @@ -133,9 +136,57 @@ extern "C" __device__ __inline__ float3 parabolic_world_normal( // ----------------------------------------------------------------------- - /**************** Aperture Helper Functions ****************/ +// ----------------------------------------------------------------------- +// Shared helpers for hexagon apertures +// +// A regular hexagon with circumradius s (vertex-to-center distance) is +// decomposed into three vertical strips in the local (x, y) frame: +// Left cap : x in [-s, -s/2) — bounded by the two left edges +// Center : x in [-s/2, s/2] — full-height rectangular band +// Right cap : x in ( s/2, s] — bounded by the two right edges +// The flat (top/bottom) edges are horizontal at y = ±(sqrt(3)/2)*s and +// the diagonal edges satisfy |y| = sqrt(3)*(|x| - s) on each cap. +// ----------------------------------------------------------------------- + +// Return true if the point (px, py) lies within a flat-top regular hexagon +// centered at the origin with circumradius s (vertex-to-center distance). +// px and py must be expressed in the hexagon's local x/y frame. +extern "C" __device__ __inline__ bool hexagon_contains(float px, float py, float s) +{ + bool is_in = false; + const float xl = 0.5f * s; + const float yl = 0.5f * sqrtf(3.0f) * s; + if (-xl <= px && px <= xl && -yl <= py && py <= yl) + { + // Center + is_in = true; + } + else if (-s <= px && px < -xl) + { + // Left side + float y1 = sqrtf(3.0f) * (px + s); + float y2 = -y1; + if (y2 <= py && py <= y1) + { + is_in = true; + } + } + else if (xl < px && px <= s) + { + // Right side + float y1 = sqrtf(3.0f) * (px - s); + float y2 = -y1; + if (y1 <= py && py <= y2) + { + is_in = true; + } + } + return is_in; +} + +// ----------------------------------------------------------------------- /**************** Optix Intersection Functions ****************/ @@ -554,39 +605,6 @@ extern "C" __global__ void __intersection__circle_flat() } } -extern "C" __device__ __inline__ bool hexagon_contains(float px, float py, float s) -{ - bool is_in = false; - const float xl = 0.5f * s; - const float yl = 0.5f * sqrtf(3.0f) * s; - if (-xl <= px && px <= xl && -yl <= py && py <= yl) - { - // Center - is_in = true; - } - else if (-s <= px && px < -xl) - { - // Left side - float y1 = sqrtf(3.0f) * (px + s); - float y2 = -y1; - if (y2 <= py && py <= y1) - { - is_in = true; - } - } - else if (xl < px && px <= s) - { - // Right side - float y1 = sqrtf(3.0f) * (px - s); - float y2 = -y1; - if (y1 <= py && py <= y2) - { - is_in = true; - } - } - return is_in; -} - extern "C" __global__ void __intersection__hexagon_flat() { const OptixCSP::GeometryDataST::Hexagon_Flat &hex = params.geometry_data_array[optixGetPrimitiveIndex()].getHexagon_Flat(); @@ -716,13 +734,15 @@ extern "C" __global__ void __intersection__circle_parabolic() float3 n; float ox, oy, oz, dx, dy, dz; - parabolic_ray_to_local(ray_orig, ray_dir, circp.center, circp.x_axis, circp.y_axis, + parabolic_ray_to_local(ray_orig, ray_dir, + circp.center, circp.x_axis, circp.y_axis, n, ox, oy, oz, dx, dy, dz); float ts[2], lxs[2], lys[2]; const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, circp.cx, circp.cy, - ray_tmin, ray_tmax, ts, lxs, lys); + ray_tmin, ray_tmax, + ts, lxs, lys); const float r2 = circp.radius * circp.radius; for (int i = 0; i < nc; ++i) @@ -794,6 +814,50 @@ extern "C" __global__ void __intersection__triangle_parabolic() extern "C" __global__ void __intersection__annulus_parabolic() { + const OptixCSP::GeometryDataST::Annulus_Parabolic &anap = + params.geometry_data_array[optixGetPrimitiveIndex()].getAnnulus_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + anap.center, anap.x_axis, anap.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + anap.cx, anap.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + const float risq = anap.ri * anap.ri; + const float rosq = anap.ro * anap.ro; + const float theta_max = 0.5f * anap.arc; + + for (int i = 0; i < nc; ++i) + { + const float rsq = lxs[i] * lxs[i] + lys[i] * lys[i]; + if (risq <= rsq && rsq <= rosq) + { + float theta = atan2f(lys[i], lxs[i]); + if (fabsf(theta) <= theta_max) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + anap.cx, anap.cy, + anap.x_axis, anap.y_axis, + n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } + } } extern "C" __global__ void __intersection__quadrilateral_parabolic() From edc138e8bb720db1ef6d7efc128521030711f5cd Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Wed, 10 Jun 2026 09:13:17 -0600 Subject: [PATCH 21/30] Commonize annulus aperture check in optix runner --- .../OptixCSP/src/shaders/intersection.cu | 74 +++++++++++-------- 1 file changed, 42 insertions(+), 32 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 7de9dd04..df65ffe8 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -186,6 +186,30 @@ extern "C" __device__ __inline__ bool hexagon_contains(float px, float py, float return is_in; } +// ----------------------------------------------------------------------- +// Shared helpers for annulus apertures +// +// An annular sector aperture is defined in the local (x, y) frame by: +// Inner radius ri : points closer than ri to the origin are excluded +// Outer radius ro : points farther than ro from the origin are excluded +// Arc angle arc : full sweep angle (radians) centered on the +x axis; +// a point at angle theta is included when |theta| <= arc/2 +// A full annulus (no arc clipping) has arc = 2*pi. +// ----------------------------------------------------------------------- + +// Return true if the point (px, py) lies within an annular sector centered at the +// origin. The sector is defined by inner radius ri, outer radius ro, and a full +// arc angle arc (in radians) symmetric about the local x-axis (i.e., theta = 0). +// px and py must be expressed in the aperture's local x/y frame. +extern "C" __device__ __inline__ bool annulus_contains(float px, float py, float ri, float ro, float arc) +{ + const float rsq = px * px + py * py; + if (rsq < ri * ri || rsq > ro * ro) + return false; + const float theta = atan2f(py, px); + return fabsf(theta) <= 0.5f * arc; +} + // ----------------------------------------------------------------------- /**************** Optix Intersection Functions ****************/ @@ -653,20 +677,15 @@ extern "C" __global__ void __intersection__annulus_flat() if (t > ray_tmin && t < ray_tmax) { float3 p = ray_orig + ray_dir * t - anf.center; - float d = length(p); - if (anf.ri <= d && d <= anf.ro) + const float px = dot(p, anf.x_axis); + const float py = dot(p, anf.y_axis); + if (annulus_contains(px, py, anf.ri, anf.ro, anf.arc)) { - float px = dot(p, anf.x_axis); - float py = dot(p, anf.y_axis); - float theta = atan2f(py, px); - if (fabsf(theta) <= 0.5f * anf.arc) - { - optixReportIntersection(t, - 0, - __float_as_uint(n.x), - __float_as_uint(n.y), - __float_as_uint(n.z)); - } + optixReportIntersection(t, + 0, + __float_as_uint(n.x), + __float_as_uint(n.y), + __float_as_uint(n.z)); } } } @@ -834,28 +853,19 @@ extern "C" __global__ void __intersection__annulus_parabolic() ray_tmin, ray_tmax, ts, lxs, lys); - const float risq = anap.ri * anap.ri; - const float rosq = anap.ro * anap.ro; - const float theta_max = 0.5f * anap.arc; - for (int i = 0; i < nc; ++i) { - const float rsq = lxs[i] * lxs[i] + lys[i] * lys[i]; - if (risq <= rsq && rsq <= rosq) + if (annulus_contains(lxs[i], lys[i], anap.ri, anap.ro, anap.arc)) { - float theta = atan2f(lys[i], lxs[i]); - if (fabsf(theta) <= theta_max) - { - const float3 wn = parabolic_world_normal(lxs[i], lys[i], - anap.cx, anap.cy, - anap.x_axis, anap.y_axis, - n); - optixReportIntersection(ts[i], 0, - __float_as_uint(wn.x), - __float_as_uint(wn.y), - __float_as_uint(wn.z)); - return; - } + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + anap.cx, anap.cy, + anap.x_axis, anap.y_axis, + n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; } } } From a5d5b7f71fe128aaef385f40824f0167489ad49e Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Wed, 10 Jun 2026 14:41:20 -0600 Subject: [PATCH 22/30] Initial implementation of parabolic triangle --- .../OptixCSP/src/core/CspElement.cpp | 44 ++++++++++---- .../OptixCSP/src/core/geometry_manager.cpp | 4 ++ .../OptixCSP/src/shaders/GeometryDataST.h | 25 ++++++++ .../OptixCSP/src/shaders/intersection.cu | 60 ++++++++++++++++++- 4 files changed, 122 insertions(+), 11 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index eb1ebe28..9d93a55d 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -265,16 +265,40 @@ GeometryDataST CspElement::toDeviceGeometryData() const v2 = tri.get_v1(); v3 = tri.get_v2(); - // given the origin and rotation, compute global coordinates of the triangle vertices - Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d v1_global = rotation_matrix * v1 + m_origin; - Vec3d v2_global = rotation_matrix * v2 + m_origin; - Vec3d v3_global = rotation_matrix * v3 + m_origin; - - GeometryDataST::Triangle_Flat heliostat(OptixCSP::toFloat3(v1_global), - OptixCSP::toFloat3(v2_global), - OptixCSP::toFloat3(v3_global)); - geometry_data.setTriangle_Flat(heliostat); + if (surface_type == SurfaceType::FLAT) + { + // given the origin and rotation, compute global coordinates of the triangle vertices + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d v1_global = rotation_matrix * v1 + m_origin; + Vec3d v2_global = rotation_matrix * v2 + m_origin; + Vec3d v3_global = rotation_matrix * v3 + m_origin; + + GeometryDataST::Triangle_Flat heliostat(OptixCSP::toFloat3(v1_global), + OptixCSP::toFloat3(v2_global), + OptixCSP::toFloat3(v3_global)); + geometry_data.setTriangle_Flat(heliostat); + } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + float3 o = OptixCSP::toFloat3(m_origin); + + // Triangle vertices are in local element XY frame (z=0). + // The 2D local x,y coords for the aperture test equal the + // Vec3d x and y components directly (since rotation is orthonormal). + float2 lv0 = make_float2((float)v1.x, (float)v1.y); + float2 lv1 = make_float2((float)v2.x, (float)v2.y); + float2 lv2 = make_float2((float)v3.x, (float)v3.y); + + GeometryDataST::Triangle_Parabolic trip(o, vx, vy, cx, cy, lv0, lv1, lv2); + geometry_data.setTriangle_Parabolic(trip); + } } if (aperture_type == ApertureType::QUADRILATERAL) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index ca95ef29..fc5d7809 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -134,6 +134,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::TRIANGLE_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::TRIANGLE_PARABOLIC); + } else { std::stringstream ss; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index d9c090fd..c74ecb99 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -231,6 +231,31 @@ namespace OptixCSP struct Triangle_Parabolic { + Triangle_Parabolic() = default; + // Vertices v0, v1, v2 are in the local XY aperture frame. + // The constructor precomputes the barycentric inverse transform so + // the aperture test reduces to two dot products with no per-ray division. + Triangle_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float2 &v0, const float2 &v1, const float2 &v2) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y) + { + const float2 e1 = make_float2(v1.x - v0.x, v1.y - v0.y); + const float2 e2 = make_float2(v2.x - v0.x, v2.y - v0.y); + const float inv_det = 1.0f / (e1.x * e2.y - e1.y * e2.x); + // u = dot(utest, float3(px, py, 1.0f)) + utest = make_float3( e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; + // v = dot(vtest, float3(px, py, 1.0f)) + vtest = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; + } + float3 center; // element origin in global coordinates + float3 x_axis; // local x axis unit vector + float3 y_axis; // local y axis unit vector + float cx; // curvature along local x axis + float cy; // curvature along local y axis + float3 utest; // precomputed row: u = dot(utest, float3(px, py, 1.0f)) + float3 vtest; // precomputed row: v = dot(vtest, float3(px, py, 1.0f)) }; struct Annulus_Parabolic diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index df65ffe8..be045ff8 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -75,7 +75,7 @@ extern "C" __device__ __inline__ int parabolic_solve( const float B = cx * ox * dx + cy * oy * dy - dz; const float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; - const float eps = 1e-6f; + const float eps = 1e-12f; int count = 0; if (fabsf(A) < eps) @@ -210,6 +210,29 @@ extern "C" __device__ __inline__ bool annulus_contains(float px, float py, float return fabsf(theta) <= 0.5f * arc; } +// ----------------------------------------------------------------------- +// Shared helpers for triangle apertures +// +// The triangle aperture is encoded as two precomputed row vectors utest and +// vtest such that the barycentric coordinates of a point (px, py) are: +// u = dot(utest, float3(px, py, 1.0f)) +// v = dot(vtest, float3(px, py, 1.0f)) +// The point is inside when u >= 0, v >= 0, and u + v <= 1. +// These rows are computed once at setup time from the three triangle vertices. +// ----------------------------------------------------------------------- + +// Return true if the 2-D point (px, py) lies inside the triangle whose +// barycentric inverse transform is encoded in utest and vtest. +extern "C" __device__ __inline__ bool triangle_contains(float px, float py, + float3 utest, float3 vtest) +{ + const float3 p = make_float3(px, py, 1.0f); + const float u = dot(utest, p); + if (u < 0.0f || u > 1.0f) return false; + const float v = dot(vtest, p); + return v >= 0.0f && (u + v) <= 1.0f; +} + // ----------------------------------------------------------------------- /**************** Optix Intersection Functions ****************/ @@ -829,6 +852,41 @@ extern "C" __global__ void __intersection__hexagon_parabolic() extern "C" __global__ void __intersection__triangle_parabolic() { + const OptixCSP::GeometryDataST::Triangle_Parabolic &trip = + params.geometry_data_array[optixGetPrimitiveIndex()].getTriangle_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + trip.center, trip.x_axis, trip.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + trip.cx, trip.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + if (triangle_contains(lxs[i], lys[i], trip.utest, trip.vtest)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + trip.cx, trip.cy, + trip.x_axis, trip.y_axis, + n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } } extern "C" __global__ void __intersection__annulus_parabolic() From b3e1dc7ceb97d436ecd43b4605dcd0d45d2979b9 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 09:28:43 -0600 Subject: [PATCH 23/30] Initial implementation of parabolic quadrilateral for optix runner --- .../OptixCSP/src/shaders/GeometryDataST.h | 37 ++++++++++++++++-- .../OptixCSP/src/shaders/intersection.cu | 38 ++++++++++++++++++- 2 files changed, 70 insertions(+), 5 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index c74ecb99..97608211 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -245,9 +245,9 @@ namespace OptixCSP const float2 e2 = make_float2(v2.x - v0.x, v2.y - v0.y); const float inv_det = 1.0f / (e1.x * e2.y - e1.y * e2.x); // u = dot(utest, float3(px, py, 1.0f)) - utest = make_float3( e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; + utest = make_float3(e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; // v = dot(vtest, float3(px, py, 1.0f)) - vtest = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; + vtest = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; } float3 center; // element origin in global coordinates float3 x_axis; // local x axis unit vector @@ -257,7 +257,7 @@ namespace OptixCSP float3 utest; // precomputed row: u = dot(utest, float3(px, py, 1.0f)) float3 vtest; // precomputed row: v = dot(vtest, float3(px, py, 1.0f)) }; - + struct Annulus_Parabolic { Annulus_Parabolic() = default; @@ -278,9 +278,38 @@ namespace OptixCSP float ro; float arc; // in radians }; - + struct Quadrilateral_Parabolic { + Quadrilateral_Parabolic() = default; + Quadrilateral_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float2 &v0, const float2 &v1, + const float2 &v2, const float2 &v3) + : center(origin), x_axis(x_ax), y_axis(y_ax), cx(curv_x), cy(curv_y) + { + const float2 e1 = make_float2(v1.x - v0.x, v1.y - v0.y); + const float2 e2 = make_float2(v2.x - v0.x, v2.y - v0.y); + const float inv_det = 1.0f / (e1.x * e2.y - e1.y * e2.x); + // u = dot(utest, float3(px, py, 1.0f)) + u1test = make_float3(e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; + // v = dot(vtest, float3(px, py, 1.0f)) + v1test = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; + + const float2 e3 = make_float2(v3.x - v0.x, v3.y - v0.y); + const float inv_det2 = 1.0f / (e2.x * e3.y - e2.y * e3.x); + u2test = make_float3(e3.y, -e3.x, v0.y * e3.x - v0.x * e3.y) * inv_det2; + v2test = make_float3(-e2.y, e2.x, v0.x * e2.y - v0.y * e2.x) * inv_det2; + } + float3 center; // element origin in global coordinates + float3 x_axis; // local x axis unit vector + float3 y_axis; // local y axis unit vector + float cx; // curvature along local x axis + float cy; // curvature along local y axis + float3 u1test; // precomputed row: u = dot(utest, float3(px, py, 1.0f)) + float3 v1test; // precomputed row: v = dot(vtest, float3(px, py, 1.0f)) + float3 u2test; + float3 v2test; }; GeometryDataST() = default; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index be045ff8..069a6777 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -228,7 +228,8 @@ extern "C" __device__ __inline__ bool triangle_contains(float px, float py, { const float3 p = make_float3(px, py, 1.0f); const float u = dot(utest, p); - if (u < 0.0f || u > 1.0f) return false; + if (u < 0.0f || u > 1.0f) + return false; const float v = dot(vtest, p); return v >= 0.0f && (u + v) <= 1.0f; } @@ -930,4 +931,39 @@ extern "C" __global__ void __intersection__annulus_parabolic() extern "C" __global__ void __intersection__quadrilateral_parabolic() { + const OptixCSP::GeometryDataST::Quadrilateral_Parabolic &quap = + params.geometry_data_array[optixGetPrimitiveIndex()].getQuadrilateral_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + quap.center, quap.x_axis, quap.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + quap.cx, quap.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + if (triangle_contains(lxs[i], lys[i], quap.u1test, quap.v1test) || + triangle_contains(lxs[i], lys[i], quap.u2test, quap.v2test)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + quap.cx, quap.cy, + quap.x_axis, quap.y_axis, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } } From 50e1260b30be7655097535f2657d61dda28d5821 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 09:46:45 -0600 Subject: [PATCH 24/30] Add missed file change to parabolic quadrilateral implementation --- .../optix_runner/OptixCSP/src/core/geometry_manager.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index fc5d7809..3d49bf15 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -155,6 +155,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::QUADRILATERAL_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::QUADRILATERAL_PARABOLIC); + } else { std::stringstream ss; From cfef853d19a45dd4d033183bbe2b6dab4e889480 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 09:56:32 -0600 Subject: [PATCH 25/30] Change quadrilateral_flat triangle decomposition to be consistent with parabolic version --- .../optix_runner/OptixCSP/src/shaders/intersection.cu | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index 069a6777..092619e7 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -592,18 +592,17 @@ extern "C" __global__ void __intersection__quadrilateral_flat() const float3 ro = optixGetObjectRayOrigin(); const float3 rd = optixGetObjectRayDirection(); + // Decompose as (p0,p1,p2) ∪ (p0,p2,p3) — same p0-p2 diagonal as the parabolic kernel. const float3 p0 = quad.p0; - const float3 p2 = quad.p2; + const float3 e02 = quad.p2 - p0; float3 e1 = quad.p1 - p0; - float3 e2 = quad.p3 - p0; - float t = _triangle_intersect(p0, e1, e2, ro, rd); + float t = _triangle_intersect(p0, e1, e02, ro, rd); if (t < optixGetRayTmin() || t > optixGetRayTmax()) { - e1 = quad.p1 - p2; - e2 = quad.p3 - p2; - t = _triangle_intersect(p2, e1, e2, ro, rd); + const float3 e3 = quad.p3 - p0; + t = _triangle_intersect(p0, e02, e3, ro, rd); } if (t < optixGetRayTmin() || t > optixGetRayTmax()) From d0e15fb32bd5fd45189792279b38e33f56761be2 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 09:56:57 -0600 Subject: [PATCH 26/30] Enforce CCW ordering of points within optix runner --- .../optix_runner/optix_runner.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/coretrace/simulation_runner/optix_runner/optix_runner.cpp b/coretrace/simulation_runner/optix_runner/optix_runner.cpp index 532638a6..83d00e89 100644 --- a/coretrace/simulation_runner/optix_runner/optix_runner.cpp +++ b/coretrace/simulation_runner/optix_runner/optix_runner.cpp @@ -329,10 +329,19 @@ RunnerStatus OptixRunner::setup_elements(const SimulationData *data) auto el_aperture = std::dynamic_pointer_cast(el->get_aperture()); assert(el_aperture != nullptr); + OptixCSP::Vec3d p0(el_aperture->x1, el_aperture->y1, 0.0); OptixCSP::Vec3d p1(el_aperture->x2, el_aperture->y2, 0.0); OptixCSP::Vec3d p2(el_aperture->x3, el_aperture->y3, 0.0); + // Ensure CCW winding (right-hand rule) required by ApertureTriangle + // Aperture always lies in x-y-plane with positive z-axis corresponding to + // the front of the geometric element so CCW winding always gives the front + // of the element. + const double signed_area = (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x); + if (signed_area < 0.0) + std::swap(p1, p2); + auto aperture = std::make_shared(p0, p1, p2); optix_el->set_aperture(aperture); @@ -348,6 +357,16 @@ RunnerStatus OptixRunner::setup_elements(const SimulationData *data) OptixCSP::Vec3d p2(el_aperture->x3, el_aperture->y3, 0.0); OptixCSP::Vec3d p3(el_aperture->x4, el_aperture->y4, 0.0); + // Ensure CCW winding using the shoelace signed-area formula. + // For a simple (non-self-intersecting) quad the sign of 2*A tells + // the winding without decomposing into triangles. + const double signed_area2 = (p0.x * p1.y - p1.x * p0.y) + + (p1.x * p2.y - p2.x * p1.y) + + (p2.x * p3.y - p3.x * p2.y) + + (p3.x * p0.y - p0.x * p3.y); + if (signed_area2 < 0.0) + std::swap(p1, p3); // reverse winding, keeping p0 and p2 fixed + auto aperture = std::make_shared(p0, p1, p2, p3); optix_el->set_aperture(aperture); From 2edbf6fb20fb572c7335954190cdd4ad290de105 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 10:50:25 -0600 Subject: [PATCH 27/30] Build fixes --- .../OptixCSP/src/core/CspElement.cpp | 50 +++++++++++++------ 1 file changed, 35 insertions(+), 15 deletions(-) diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 9d93a55d..9f4dea55 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -292,9 +292,9 @@ GeometryDataST CspElement::toDeviceGeometryData() const // Triangle vertices are in local element XY frame (z=0). // The 2D local x,y coords for the aperture test equal the // Vec3d x and y components directly (since rotation is orthonormal). - float2 lv0 = make_float2((float)v1.x, (float)v1.y); - float2 lv1 = make_float2((float)v2.x, (float)v2.y); - float2 lv2 = make_float2((float)v3.x, (float)v3.y); + float2 lv0 = make_float2((float)v1[0], (float)v1[1]); + float2 lv1 = make_float2((float)v2[0], (float)v2[1]); + float2 lv2 = make_float2((float)v3[0], (float)v3[1]); GeometryDataST::Triangle_Parabolic trip(o, vx, vy, cx, cy, lv0, lv1, lv2); geometry_data.setTriangle_Parabolic(trip); @@ -312,18 +312,38 @@ GeometryDataST CspElement::toDeviceGeometryData() const p3 = quad.get_p2(); p4 = quad.get_p3(); - // given the origin and rotation, compute global coordinates of the triangle vertices - Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d p1_global = rotation_matrix * p1 + m_origin; - Vec3d p2_global = rotation_matrix * p2 + m_origin; - Vec3d p3_global = rotation_matrix * p3 + m_origin; - Vec3d p4_global = rotation_matrix * p4 + m_origin; - - GeometryDataST::Quadrilateral_Flat heliostat(OptixCSP::toFloat3(p1_global), - OptixCSP::toFloat3(p2_global), - OptixCSP::toFloat3(p3_global), - OptixCSP::toFloat3(p4_global)); - geometry_data.setQuadrilateral_Flat(heliostat); + if (surface_type == SurfaceType::FLAT) + { + // given the origin and rotation, compute global coordinates of the triangle vertices + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d p1_global = rotation_matrix * p1 + m_origin; + Vec3d p2_global = rotation_matrix * p2 + m_origin; + Vec3d p3_global = rotation_matrix * p3 + m_origin; + Vec3d p4_global = rotation_matrix * p4 + m_origin; + + GeometryDataST::Quadrilateral_Flat heliostat(OptixCSP::toFloat3(p1_global), + OptixCSP::toFloat3(p2_global), + OptixCSP::toFloat3(p3_global), + OptixCSP::toFloat3(p4_global)); + geometry_data.setQuadrilateral_Flat(heliostat); + } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + float3 o = OptixCSP::toFloat3(m_origin); + + GeometryDataST::Quadrilateral_Parabolic heliostat( + o, vx, vy, cx, cy, + OptixCSP::toFloat3(p1), OptixCSP::toFloat3(p2), + OptixCSP::toFloat3(p3), OptixCSP::toFloat3(p4)); + geometry_data.setQuadrilateral_Parabolic(heliostat); + } } if (aperture_type == ApertureType::CIRCLE) From ad6207ab7fe865a946b140b0fbf8d1cb82e814e6 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 11:01:38 -0600 Subject: [PATCH 28/30] Add tests for quadrilateral and triangle winding order in optix runner --- .../geometry_intersection_test.cpp | 145 ++++++++++++++++++ 1 file changed, 145 insertions(+) diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index 489da846..9403eceb 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -251,6 +251,78 @@ TEST(OptixRunner, FlatEquilateralTriangle) TEST(OptixRunner, FlatTriangle) { + // Same triangle as FlatTriangle but vertices supplied in CCW order: + // (0,0)->(2,0)->(1,2). The optix_runner's cross-product sign check sees + // a positive area and does NOT swap, so the CCW path is taken directly. + // Results should be identical to FlatTriangle. + const double x1 = 0.0, x2 = 2.0, x3 = 1.0; + const double y1 = 0.0, y2 = 0.0, y3 = 2.0; + const double ROT_DEG = 110.0; + auto surf = make_surface(); + auto aper = make_aperture(x1, y1, x2, y2, x3, y3); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, FlatTriangle_CW) +{ + // Vertices are in CW order (signed area = -2). The optix_runner's + // cross-product sign check detects this and swaps p1<->p2 to produce + // CCW winding before building the geometry. const double x1 = 0.0, x2 = 1.0, x3 = 2.0 * x2; const double y1 = 0.0, y2 = 2.0, y3 = y1; const double ROT_DEG = 110.0; @@ -381,6 +453,79 @@ TEST(OptixRunner, FlatQuadrilateral) << ", false negatives: " << fneg << std::endl; } +TEST(OptixRunner, FlatQuadrilateral_CW) +{ + // Same parallelogram as FlatQuadrilateral but vertices supplied in + // clockwise winding order: (0,0)->(1,2)->(4,2)->(3,0). + // The x1-x3 diagonal is already interior for this ordering, so + // ensure_valid_diagonal() leaves the vertices unchanged (CW). + // The optix_runner's shoelace sign check detects the negative area and + // swaps p1<->p3 to restore CCW winding before building the geometry. + // Results should be identical to FlatQuadrilateral. + const double x1 = 0.0, x2 = 1.0, x3 = 4.0, x4 = 3.0; + const double y1 = 0.0, y2 = 2.0, y3 = 2.0, y4 = 0.0; + const double ROT_DEG = -45.0; + auto surf = make_surface(); + auto aper = make_aperture( + x1, y1, x2, y2, x3, y3, x4, y4); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + TEST(OptixRunner, ParabolaRectangle) { constexpr double CX = 0.5; From 0f9fa96e580d6a80335c8df220373ccc429ae591 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 15:29:03 -0600 Subject: [PATCH 29/30] Add correction to irregular quadrilateral aperture to enforce it being a simple polygon --- coretrace/simulation_data/aperture.cpp | 49 ++++++++++++++++++++++++++ coretrace/simulation_data/aperture.hpp | 20 +++++++++-- 2 files changed, 67 insertions(+), 2 deletions(-) diff --git a/coretrace/simulation_data/aperture.cpp b/coretrace/simulation_data/aperture.cpp index 5982dd51..820d88b0 100644 --- a/coretrace/simulation_data/aperture.cpp +++ b/coretrace/simulation_data/aperture.cpp @@ -670,6 +670,7 @@ namespace SolTrace::Data x3(x3), y3(y3), x4(x4), y4(y4) { + ensure_valid_diagonal(); } IrregularQuadrilateral::IrregularQuadrilateral(const nlohmann::ordered_json &jnode) @@ -683,6 +684,7 @@ namespace SolTrace::Data this->y3 = jnode.at("y3"); this->x4 = jnode.at("x4"); this->y4 = jnode.at("y4"); + ensure_valid_diagonal(); } double IrregularQuadrilateral::aperture_area() const @@ -750,6 +752,53 @@ namespace SolTrace::Data return inquad(x1, y1, x2, y2, x3, y3, x4, y4, x, y); } + void IrregularQuadrilateral::ensure_valid_diagonal() + { + // Returns true when x2 and x4 are on opposite sides of the x1-x3 line, + // i.e. the x1-x3 diagonal is interior to the polygon. + const auto diagonal_ok = [this]() + { + const double ex = x3 - x1, ey = y3 - y1; + const double d2 = ex * (y2 - y1) - ey * (x2 - x1); + const double d4 = ex * (y4 - y1) - ey * (x4 - x1); + return d2 * d4 < 0.0; + }; + + if (!diagonal_ok()) + { + // Vertices are not in simple-polygon order (self-intersecting / + // bowtie input). Re-sort by angle from the centroid to recover a + // valid CCW traversal. + const double cx = 0.25 * (x1 + x2 + x3 + x4); + const double cy = 0.25 * (y1 + y2 + y3 + y4); + double pts[4][2] = {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}; + std::sort(std::begin(pts), std::end(pts), + [cx, cy](const double *a, const double *b) + { + return std::atan2(a[1] - cy, a[0] - cx) < + std::atan2(b[1] - cy, b[0] - cx); + }); + x1 = pts[0][0]; y1 = pts[0][1]; + x2 = pts[1][0]; y2 = pts[1][1]; + x3 = pts[2][0]; y3 = pts[2][1]; + x4 = pts[3][0]; y4 = pts[3][1]; + + // The sort produces a valid CCW polygon but the starting vertex is + // arbitrary. If the reflex vertex landed at position x1 or x3, the + // x1-x3 diagonal is still exterior. One cyclic shift swaps which + // pair of opposite vertices forms the diagonal (x1-x3 becomes the + // old x2-x4), which is guaranteed to be interior for a simple polygon. + if (!diagonal_ok()) + { + const double tx = x1, ty = y1; + x1 = x2; y1 = y2; + x2 = x3; y2 = y3; + x3 = x4; y3 = y4; + x4 = tx; y4 = ty; + } + } + } + aperture_ptr IrregularQuadrilateral::make_copy() const { // Invokes the implicit copy constructor diff --git a/coretrace/simulation_data/aperture.hpp b/coretrace/simulation_data/aperture.hpp index c460234f..d3c8fbb7 100644 --- a/coretrace/simulation_data/aperture.hpp +++ b/coretrace/simulation_data/aperture.hpp @@ -709,7 +709,14 @@ namespace SolTrace::Data struct IrregularQuadrilateral : public Aperture { - // Locations of the 4 vertices + // Locations of the 4 vertices. + // + // The quadrilateral must be a SIMPLE (non-self-intersecting) polygon. + // The vertices may be supplied in any order; the constructors call + // ensure_valid_diagonal() which sorts them by angle from the centroid + // into a consistent CCW traversal and ensures the x1-x3 diagonal is + // interior. All triangle-decomposition logic (inquad, the flat and + // parabolic intersection kernels) depends on this diagonal being interior. double x1; double y1; double x2; @@ -786,7 +793,16 @@ namespace SolTrace::Data */ virtual std::tuple, std::vector> triangulation() const override; - }; + + /** + * @brief Ensure the x1-x3 diagonal is interior to the quadrilateral. + * All runners decompose the quad into triangles (x1,x2,x3) and + * (x1,x3,x4) along this diagonal. If x2 and x4 lie on the same side + * of the x1-x3 line the diagonal is exterior and the decomposition + * is incorrect; swapping x2<->x4 moves the reflex vertex onto the + * shared diagonal endpoints and makes x1-x3 interior. + */ + void ensure_valid_diagonal(); /** * @brief Test if point is inside triangle defined by three vertices From ece31321cd34bbbbfc4d02797c072e049bc2ae00 Mon Sep 17 00:00:00 2001 From: Jonathan Maack Date: Thu, 11 Jun 2026 15:55:39 -0600 Subject: [PATCH 30/30] Build fixes; add tests for irregular quadrilateral vertex ordering --- coretrace/simulation_data/aperture.cpp | 7 +- coretrace/simulation_data/aperture.hpp | 1 + .../simulation_data/aperture_test.cpp | 87 +++++++++++++++++++ 3 files changed, 92 insertions(+), 3 deletions(-) diff --git a/coretrace/simulation_data/aperture.cpp b/coretrace/simulation_data/aperture.cpp index 820d88b0..a107f1c5 100644 --- a/coretrace/simulation_data/aperture.cpp +++ b/coretrace/simulation_data/aperture.cpp @@ -771,9 +771,10 @@ namespace SolTrace::Data // valid CCW traversal. const double cx = 0.25 * (x1 + x2 + x3 + x4); const double cy = 0.25 * (y1 + y2 + y3 + y4); - double pts[4][2] = {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}; - std::sort(std::begin(pts), std::end(pts), - [cx, cy](const double *a, const double *b) + using P2 = std::array; + std::array pts = {P2{x1, y1}, P2{x2, y2}, P2{x3, y3}, P2{x4, y4}}; + std::sort(pts.begin(), pts.end(), + [cx, cy](const P2 &a, const P2 &b) { return std::atan2(a[1] - cy, a[0] - cx) < std::atan2(b[1] - cy, b[0] - cx); diff --git a/coretrace/simulation_data/aperture.hpp b/coretrace/simulation_data/aperture.hpp index d3c8fbb7..bef66f4c 100644 --- a/coretrace/simulation_data/aperture.hpp +++ b/coretrace/simulation_data/aperture.hpp @@ -803,6 +803,7 @@ namespace SolTrace::Data * shared diagonal endpoints and makes x1-x3 interior. */ void ensure_valid_diagonal(); + }; /** * @brief Test if point is inside triangle defined by three vertices diff --git a/google-tests/unit-tests/simulation_data/aperture_test.cpp b/google-tests/unit-tests/simulation_data/aperture_test.cpp index ddf0f642..50944ef1 100644 --- a/google-tests/unit-tests/simulation_data/aperture_test.cpp +++ b/google-tests/unit-tests/simulation_data/aperture_test.cpp @@ -600,6 +600,93 @@ TEST(Aperture, IrregularQuadrilateral) std::cout << "Expected fraction of hits: " << frac << std::endl; } +TEST(Aperture, IrregularQuadrilateral_VertexOrder) +{ + // Verify that the same parallelogram described with four different vertex + // orderings — CCW, CW (reverse), a cyclic rotation, and a self-intersecting + // "bowtie" ordering — all produce identical area and is_in results after + // ensure_valid_diagonal() normalises the representation. + // + // Canonical CCW vertices: (0,0) -> (3,0) -> (4,2) -> (1,2) + const double TOL = 1e-12; + const double AREA = 6.0; // parallelogram: base 3 * height 2 + // One interior and one exterior point, chosen away from edges. + const double IX = 2.5, IY = 1.0; // clearly inside + const double OX = 4.5, OY = 1.0; // clearly outside + + // Helper: extract the four normalised vertices as an array of (x,y) pairs. + using V4 = std::array, 4>; + auto get_verts = [](const aperture_ptr &ap) -> V4 { + const auto *q = dynamic_cast(ap.get()); + return {{ {q->x1,q->y1}, {q->x2,q->y2}, {q->x3,q->y3}, {q->x4,q->y4} }}; + }; + // Returns true when both arrays contain the same four (x,y) pairs (any order). + auto same_vertex_set = [&](const V4 &a, const V4 &b) -> bool { + for (const auto &va : a) { + bool found = false; + for (const auto &vb : b) + if (std::abs(va.first - vb.first) < TOL && + std::abs(va.second - vb.second) < TOL) + { found = true; break; } + if (!found) return false; + } + return true; + }; + // Returns true when the x1-x3 diagonal of q is interior to the quad, + // i.e. x2 and x4 are on strictly opposite sides of the x1-x3 line. + auto diagonal_interior = [&](const aperture_ptr &ap) -> bool { + const auto *q = dynamic_cast(ap.get()); + const double ex = q->x3 - q->x1, ey = q->y3 - q->y1; + const double d2 = ex * (q->y2 - q->y1) - ey * (q->x2 - q->x1); + const double d4 = ex * (q->y4 - q->y1) - ey * (q->x4 - q->x1); + return d2 * d4 < 0.0; + }; + + // 1) CCW (canonical) + auto q_ccw = make_aperture( + 0.0, 0.0, 3.0, 0.0, 4.0, 2.0, 1.0, 2.0); + EXPECT_NEAR(q_ccw->aperture_area(), AREA, TOL) << "CCW area"; + EXPECT_TRUE( q_ccw->is_in(IX, IY)) << "CCW inside"; + EXPECT_FALSE(q_ccw->is_in(OX, OY)) << "CCW outside"; + EXPECT_TRUE(diagonal_interior(q_ccw)) << "CCW diagonal interior"; + const V4 canonical = get_verts(q_ccw); + + // 2) CW (full reversal): (0,0) -> (1,2) -> (4,2) -> (3,0) + // The x1-x3 diagonal (0,0)-(4,2) is already interior for this ordering, + // so ensure_valid_diagonal() is a no-op and the vertices stay CW. + auto q_cw = make_aperture( + 0.0, 0.0, 1.0, 2.0, 4.0, 2.0, 3.0, 0.0); + EXPECT_NEAR(q_cw->aperture_area(), AREA, TOL) << "CW area"; + EXPECT_TRUE( q_cw->is_in(IX, IY)) << "CW inside"; + EXPECT_FALSE(q_cw->is_in(OX, OY)) << "CW outside"; + EXPECT_TRUE(diagonal_interior(q_cw)) << "CW diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_cw))) << "CW vertex set"; + + // 3) Cyclic rotation of CCW: start from (3,0) + // (3,0) -> (4,2) -> (1,2) -> (0,0) + // The x1-x3 diagonal (3,0)-(1,2) is interior, so ensure_valid_diagonal() + // is again a no-op. + auto q_rot = make_aperture( + 3.0, 0.0, 4.0, 2.0, 1.0, 2.0, 0.0, 0.0); + EXPECT_NEAR(q_rot->aperture_area(), AREA, TOL) << "rotated area"; + EXPECT_TRUE( q_rot->is_in(IX, IY)) << "rotated inside"; + EXPECT_FALSE(q_rot->is_in(OX, OY)) << "rotated outside"; + EXPECT_TRUE(diagonal_interior(q_rot)) << "rotated diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_rot))) << "rotated vertex set"; + + // 4) Self-intersecting "bowtie" order: (0,0) -> (4,2) -> (3,0) -> (1,2) + // The x1-x3 diagonal test fails (same-sign cross products), so + // ensure_valid_diagonal() angle-sorts and cyclic-shifts to recover a + // valid simple polygon with an interior x1-x3 diagonal. + auto q_bowtie = make_aperture( + 0.0, 0.0, 4.0, 2.0, 3.0, 0.0, 1.0, 2.0); + EXPECT_NEAR(q_bowtie->aperture_area(), AREA, TOL) << "bowtie area"; + EXPECT_TRUE( q_bowtie->is_in(IX, IY)) << "bowtie inside"; + EXPECT_FALSE(q_bowtie->is_in(OX, OY)) << "bowtie outside"; + EXPECT_TRUE(diagonal_interior(q_bowtie)) << "bowtie diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_bowtie))) << "bowtie vertex set"; +} + TEST(Aperture, MakeApertureFromType) { const double TOL = 1e-12;