river_sss.h

Go to the documentation of this file.

//
// Created by julian on 1/15/26.
//
 
#ifndef UG4_RIVER_SSS_H
#define UG4_RIVER_SSS_H
#include "common/common.h"
#include "common/math/misc/shapes.h"
#include "lib_disc/spatial_disc/disc_util/fv1_sss.h"
#include "../d3f_sss.h"
namespace ug{
    namespace d3f{
 
        template <int dim>
        bool boxes_overlap(const AABox<MathVector<dim>>& a,
                           const AABox<MathVector<dim>>& b)
        {
            for (int d = 0; d < dim; ++d)
            {
                if (a.max[d] < b.min[d] || a.min[d] > b.max[d])
                    return false; // no overlap in this dimension
            }
            return true; // overlaps in all dimensions
        }
 
template <int dim>
void expand_bbox(AABox<MathVector<dim>>& bbox,
                 const AABox<MathVector<dim>>& other)
{
    for (int d = 0; d < dim; ++d)
    {
        if (other.min[d] < bbox.min[d]) bbox.min[d] = other.min[d];
        if (other.max[d] > bbox.max[d]) bbox.max[d] = other.max[d];
    }
}
 
template <int dim>
int longest_axis(const AABox<MathVector<dim>>& box)
{
    int axis = 0;
    number max_len = box.max[0] - box.min[0];
 
    for (int d = 1; d < dim; ++d)
    {
        number len = box.max[d] - box.min[d];
        if (len > max_len)
        {
            max_len = len;
            axis = d;
        }
    }
    return axis;
}
 
    template <int dim>
AABox<MathVector<dim>>
LineBoundingBox(const MathVector<dim>& a,
                const MathVector<dim>& b)
    {
        AABox<MathVector<dim>> bb;
    // initialize min and max component-wise
    for (int d = 0; d < dim; ++d)
    {
        if (a[d] < b[d])
        {
            bb.min[d] = a[d];
            bb.max[d] = b[d];
        }
        else
        {
            bb.min[d] = b[d];
            bb.max[d] = a[d];
        }
    }
 
    return bb;
    }
 
    template <int dim>
struct LineBVHNode
    {
        AABox<MathVector<dim>> bbox; 
        std::vector<size_t> indices; 
        SmartPtr<LineBVHNode> left;  
        SmartPtr<LineBVHNode> right; 
 
        bool is_leaf() const { return !left && !right; }
    };
 
template <int dim, typename TLineData>
class LineBVHTree
{
public:
    using line_type = FVLineSourceOrSink<dim, TLineData>;
 
    void build(std::vector<SmartPtr<line_type>>& lines)
    {
        std::vector<size_t> ids(lines.size());
        for (size_t i = 0; i < ids.size(); ++i)
            ids[i] = i;
        m_root = build_node(lines, ids, 0);
    }
 
    void query(const AABox<MathVector<dim>>& box,
               std::vector<size_t>& out) const
    {
        query_node(m_root.get(), box, out);
    }
 
private:
    SmartPtr<LineBVHNode<dim>> m_root;
 
    SmartPtr<LineBVHNode<dim>>
    build_node(std::vector<SmartPtr<line_type>>& lines,
               std::vector<size_t>& ids,
               int depth)
    {
        auto node = make_sp(new LineBVHNode<dim>());
 
        // 1. Compute the combined bounding box for all lines in this node
        if (!ids.empty())
        {
            node->bbox = LineBoundingBox<dim>(
                lines[ids[0]]->from_position(),
                lines[ids[0]]->to_position());
 
            for (size_t i = 1; i < ids.size(); ++i)
            {
                AABox<MathVector<dim>> bb = LineBoundingBox<dim>(
                    lines[ids[i]]->from_position(),
                    lines[ids[i]]->to_position());
                expand_bbox<dim>(node->bbox, bb);
            }
        }
 
        // 2. Base case: If few lines remain, create a leaf node
        if (ids.size() <= 8)
        {
            node->indices = ids;
            return node;
        }
 
        // 3. Choose split axis and sort lines to partition them
        int axis = longest_axis<dim>(node->bbox);
 
        std::sort(ids.begin(), ids.end(),
            [&](size_t a, size_t b)
            {
                return lines[a]->from_position()[axis] <
                       lines[b]->from_position()[axis];
            });
 
        // 4. Split and recurse
        size_t mid = ids.size() / 2;
        std::vector<size_t> left(ids.begin(), ids.begin() + mid);
        std::vector<size_t> right(ids.begin() + mid, ids.end());
 
        node->left = build_node(lines, left, depth + 1);
        node->right = build_node(lines, right, depth + 1);
        return node;
    }
 
    void query_node(const LineBVHNode<dim>* node,
                    const AABox<MathVector<dim>>& box,
                    std::vector<size_t>& out) const
    {
        // Early exit if node is empty or query box doesn't overlap the node's volume
        if (!node || !boxes_overlap<dim>(node->bbox, box))
            return;
 
        // If leaf, add all contained indices to the output
        if (node->is_leaf())
        {
            out.insert(out.end(),
                       node->indices.begin(),
                       node->indices.end());
            return;
        }
 
        // Recurse into children
        query_node(node->left.get(), box, out);
        query_node(node->right.get(), box, out);
    }
};
 
    template <int dim, typename TPointData, typename TLineData = TPointData>
class FVSingularSourcesAndSinksKD
{
public:
    using point_sss_type = FVPointSourceOrSink<dim, TPointData>;
    using line_sss_type  = FVLineSourceOrSink<dim, TLineData>;
 
private:
    std::vector<SmartPtr<point_sss_type>> ListP; 
    std::vector<SmartPtr<line_sss_type>>  ListL; 
 
    LineBVHTree<dim, TLineData> m_bvh;           
    bool m_bvh_valid = false;                    
 
public:
    void add_point(SmartPtr<point_sss_type> p)
    {
        ListP.push_back(p);
    }
 
    void add_line(SmartPtr<line_sss_type> l)
    {
        ListL.push_back(l);
        m_bvh_valid = false; // Mark for rebuild on next query
    }
 
    void build_spatial_index()
    {
        m_bvh.build(ListL);
        m_bvh_valid = true;
    }
 
    size_t num_lines() const { return ListL.size(); }
    line_sss_type* line(size_t i) { return ListL[i].get(); }
 
    template <typename TElem, typename TAAPos, typename TFVGeom>
    class line_iterator
    {
        using master_type = FVSingularSourcesAndSinksKD;
 
        master_type* m_sss;
        std::vector<size_t> m_candidates; 
        size_t m_pos = 0;                 
 
        TElem* m_elem;
        Grid& m_grid;
        TAAPos& m_aaPos;
        const TFVGeom& m_geo;
        size_t m_co;
 
        MathVector<dim> m_ls, m_le;
 
    public:
        line_iterator(master_type* sss,
                      TElem* elem,
                      Grid& grid,
                      TAAPos& aaPos,
                      const TFVGeom& geo,
                      size_t co)
            : m_sss(sss), m_elem(elem), m_grid(grid),
              m_aaPos(aaPos), m_geo(geo), m_co(co)
        {
            if (!m_sss->m_bvh_valid)
                m_sss->build_spatial_index();
 
            // Find all line segments whose bounding box overlaps the element's bounding box
            auto elem_bb = CalculateBoundingBox(elem, aaPos);
            m_sss->m_bvh.query(elem_bb, m_candidates);
            
            // Advance to the first segment that actually intersects the element geometry
            advance();
        }
 
        bool is_over() const { return m_pos >= m_candidates.size(); }
 
        line_sss_type* operator*()
        {
            return m_sss->ListL[m_candidates[m_pos]].get();
        }
 
        const MathVector<dim>& seg_start() const { return m_ls; }
        const MathVector<dim>& seg_end()   const { return m_le; }
 
        line_iterator& operator++()
        {
            ++m_pos;
            advance();
            return *this;
        }
 
    private:
        void advance()
        {
            while (m_pos < m_candidates.size())
            {
                auto* l = m_sss->ListL[m_candidates[m_pos]].get();
                // BVH only tells us if boxes overlap. corresponds_to() performs the 
                // expensive exact intersection test between the line segment and the FV cell.
                if (l->corresponds_to(m_elem, m_grid, m_aaPos,
                                      m_geo, m_co, m_ls, m_le))
                    return;
                ++m_pos;
            }
        }
    };
};
 
template <int dim>
class RiverSingularSourcesAndSinks
    : public FVSingularSourcesAndSinksKD<dim, river_point_sss_data<dim> ,river_line_sss_data<dim> >
{};
 
} // namespace d3f
} // namespace ug
#endif //UG4_RIVER_SSS_H