#include <calculate_integral.hpp>

void expand_integral_on_u(uint8_t                                           q,
                          const pcurve_relation_graph_t&                    pcurve_relation_graph,
                          uint32_t                                          subface_index,
                          double                                            v,
                          double                                            u_start,
                          double                                            u_end,
                          bool                                              with_start,
                          bool                                              with_end,
                          double                                            w_u,
                          std::function<double(int, int, double, double)>&& integrator_x,
                          std::function<double(int, int, double, double)>&& integrator_w,
                          std::vector<Eigen::Vector3d>&                     integral_points)
{
    auto results = v_line_subface_intersection(pcurve_relation_graph, subface_index, v, u_start, u_end, with_start, with_end);
    assert(results.size() % 2 == 0);

    integral_points.reserve(integral_points.size() + q * results.size() / 2);
    for (auto iter = results.begin(), next = iter + 1; iter != results.end(); iter += 2, next += 2) {
        // iff degenerate, ignore it
        if (sorted_double_equal(*iter, *(iter + 1))) continue;

        for (auto i = 0; i < q; ++i)
            integral_points.emplace_back(integrator_x(q, i, *iter, *(iter + 1)),
                                         v,
                                         w_u * integrator_w(q, i, *iter, *(iter + 1)));
    }
}

void expand_integral_on_v(uint8_t                                           q,
                          const pcurve_relation_graph_t&                    pcurve_relation_graph,
                          uint32_t                                          subface_index,
                          double                                            u,
                          double                                            v_start,
                          double                                            v_end,
                          bool                                              with_start,
                          bool                                              with_end,
                          double                                            w_v,
                          std::function<double(int, int, double, double)>&& integrator_x,
                          std::function<double(int, int, double, double)>&& integrator_w,
                          std::vector<Eigen::Vector3d>&                     integral_points)
{
    auto results = u_line_subface_intersection(pcurve_relation_graph, subface_index, u, v_start, v_end, with_start, with_end);
    assert(results.size() % 2 == 0);

    integral_points.reserve(integral_points.size() + q * results.size() / 2);
    for (auto iter = results.begin(), next = iter + 1; iter != results.end(); iter += 2, next += 2) {
        // iff degenerate, ignore it
        if (sorted_double_equal(*iter, *(iter + 1))) continue;

        for (auto i = 0; i < q; ++i)
            integral_points.emplace_back(u,
                                         integrator_x(q, i, *iter, *(iter + 1)),
                                         w_v * integrator_w(q, i, *iter, *(iter + 1)));
    }
}