116 support all geometries gpu runner

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,54 @@ 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);
+            using P2 = std::array<double, 2>;
+            std::array<P2, 4> 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);
+                      });
+            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

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,6 +793,16 @@ namespace SolTrace::Data
          */
         virtual std::tuple<std::vector<double>, std::vector<int>>
         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();
     };
 
     /**

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;
 
@@ -264,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[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);
+        }
     }
 
     if (aperture_type == ApertureType::QUADRILATERAL)
@@ -287,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)
@@ -307,39 +352,79 @@ 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);
+            geometry_data.setCircle_Parabolic(heliostat);
+        }
     }
 
     if (aperture_type == ApertureType::HEXAGON)
     {
+        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());
+
         ApertureHexagon hex = static_cast<ApertureHexagon &>(*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 (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)
     {
+        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());
+
         ApertureAnnulus anf = static_cast<ApertureAnnulus &>(*m_aperture);
         float radius_in = anf.get_radius_inner();
         float radius_out = anf.get_radius_outer();
         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, radius_in, radius_out, arc);
+            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;

coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp

@@ -46,6 +46,10 @@ void GeometryManager::collect_geometry_info(const std::vector<std::shared_ptr<Cs
             {
                 sbt_offset = static_cast<uint32_t>(OpticalEntityType::ANNULUS_FLAT);
             }
+            else if (element->get_surface_type() == SurfaceType::PARABOLIC)
+            {
+                sbt_offset = static_cast<uint32_t>(OpticalEntityType::ANNULUS_PARABOLIC);
+            }
             else
             {
                 std::stringstream ss;
@@ -63,6 +67,10 @@ void GeometryManager::collect_geometry_info(const std::vector<std::shared_ptr<Cs
             {
                 sbt_offset = static_cast<uint32_t>(OpticalEntityType::CIRCLE_FLAT);
             }
+            else if (element->get_surface_type() == SurfaceType::PARABOLIC)
+            {
+                sbt_offset = static_cast<uint32_t>(OpticalEntityType::CIRCLE_PARABOLIC);
+            }
             else
             {
                 std::stringstream ss;
@@ -80,6 +88,10 @@ void GeometryManager::collect_geometry_info(const std::vector<std::shared_ptr<Cs
             {
                 sbt_offset = static_cast<uint32_t>(OpticalEntityType::HEXAGON_FLAT);
             }
+            else if (element->get_surface_type() == SurfaceType::PARABOLIC)
+            {
+                sbt_offset = static_cast<uint32_t>(OpticalEntityType::HEXAGON_PARABOLIC);
+            }
             else
             {
                 std::stringstream ss;
@@ -122,6 +134,10 @@ void GeometryManager::collect_geometry_info(const std::vector<std::shared_ptr<Cs
             {
                 sbt_offset = static_cast<uint32_t>(OpticalEntityType::TRIANGLE_FLAT);
             }
+            else if (element->get_surface_type() == SurfaceType::PARABOLIC)
+            {
+                sbt_offset = static_cast<uint32_t>(OpticalEntityType::TRIANGLE_PARABOLIC);
+            }
             else
             {
                 std::stringstream ss;
@@ -139,6 +155,10 @@ void GeometryManager::collect_geometry_info(const std::vector<std::shared_ptr<Cs
             {
                 sbt_offset = static_cast<uint32_t>(OpticalEntityType::QUADRILATERAL_FLAT);
             }
+            else if (element->get_surface_type() == SurfaceType::PARABOLIC)
+            {
+                sbt_offset = static_cast<uint32_t>(OpticalEntityType::QUADRILATERAL_PARABOLIC);
+            }
             else
             {
                 std::stringstream ss;

coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp

@@ -45,7 +45,12 @@ const std::map<OpticalEntityType, std::string> 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) {}