flooding_solver.cc

#include <bitset>

#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/bisection.hh>
#include <vtestbed/core/flooding_solver.hh>
#include <vtestbed/geometry/line_segment.hh>
#include <vtestbed/geometry/polyhedron.hh>
#include <vtestbed/geometry/tetrahedron.hh>

template <class T>
vtb::core::Flooding_solver<T>::Flooding_solver(const polyhedron_type& polyhedron):
_polyhedron(polyhedron) {
    this->_z_bounds = this->_polyhedron.bounds(2);
}

template <class T>
void
vtb::core::Flooding_solver<T>::volume(T rhs) {
    this->_volume = rhs;
    this->_level = calc_level();
    const auto& result = calc_volume(level());
    this->_volume = result.volume();
    this->_centre = result.point();
}

template <class T>
void
vtb::core::Flooding_solver<T>::level(T rhs) {
    this->_level = rhs;
    const auto& result = calc_volume(level());
    this->_volume = result.volume();
    this->_centre = result.point();
}

template <class T>
auto
vtb::core::Flooding_solver<T>::calc_level() -> T {
    Bisection<T,1> bisection;
    // fluid level precision (1 cm)
    bisection.argument_precision(T{1e-2});
    // fluid volume precision (1 cm^3)
    bisection.function_precision(T{1e-6});
    bisection.max_iterations(20);
    return bisection(
        z_bounds().min(), z_bounds().max(),
        [this] (T level) { return calc_volume(level).volume() - this->_volume; }
    );
}

template <class T>
auto
vtb::core::Flooding_solver<T>::calc_volume(T level) -> Volume<T> {
    Volume<T> volume;
    if (level <= z_bounds().min()) { return volume; }
    if (level > z_bounds().max()) { level = z_bounds().max(); }
    using vtb::geometry::Tetrahedron;
    using vec3 = Vector<T,3>;
    using line_segment = vtb::geometry::Line_segment<T,3>;
    using vtb::geometry::intersection_point;
    const auto& vertices = polyhedron().vertices();
    const auto& faces = polyhedron().faces();
    vec3 origin(0,0,level);
    #if defined(DEBUG_VOLUME)
    std::ofstream("origin.tmp") << origin;
    std::ofstream intersectionPoints("intersection.tmp");
    std::ofstream belowPoints("below.tmp");
    std::ofstream borderPoints("border.tmp");
    #endif
    for (const auto& p : faces) {
        vec3 v[3] = {vertices[p[0]], vertices[p[1]], vertices[p[2]]};
        std::bitset<3> above;
        above[0] = v[0](2) > level;
        above[1] = v[1](2) > level;
        above[2] = v[2](2) > level;
        if (above.none()) {
            // All three vertices are under water.
            Tetrahedron<T,3> t(origin, v[0], v[1], v[2]);
            volume.accumulate(signed_volume(t), centroid(t));
            #if defined(DEBUG_VOLUME)
            belowPoints << t;
            #endif
        } else if (above.count() < 3) {
            // One or two vertices are under water.
            int outstanding = 0;
            int i1 = 0;
            int i2 = 0;
            const bool x01 = above[0] ^ above[1];
            const bool x02 = above[0] ^ above[2];
            const bool x12 = above[1] ^ above[2];
            if (x01 && x02) { outstanding = 0, i1 = 1, i2 = 2; }
            else if (x01 && x12) { outstanding = 1, i1 = 2, i2 = 0; }
            else { outstanding = 2, i1 = 0, i2 = 1; }
            const vec3 outstandingVertex = v[outstanding];
            vec3 v_new[3];
            v_new[outstanding] = outstandingVertex;
            v_new[i1] = intersection_point<2>(line_segment{outstandingVertex, v[i1]}, level);
            v_new[i2] = intersection_point<2>(line_segment{outstandingVertex, v[i2]}, level);
            #if defined(DEBUG_VOLUME)
            intersectionPoints << v_new[i1] << '\n';
            intersectionPoints << v_new[i2] << '\n';
            #endif
            // Correct vertices order in tetrahedron constructor is needed
            // for correct volume sign.
            if (outstandingVertex(2) > level) {
                Tetrahedron<T,3> t1(origin, v_new[i1], v[i1], v[i2]);
                volume.accumulate(signed_volume(t1), centroid(t1));
                #if defined(DEBUG_VOLUME)
                borderPoints << t1;
                #endif
                Tetrahedron<T,3> t2(origin, v_new[i1], v[i2], v_new[i2]);
                volume.accumulate(signed_volume(t2), centroid(t2));
                #if defined(DEBUG_VOLUME)
                borderPoints << t2;
                #endif
            } else {
                Tetrahedron<T,3> t(origin, outstandingVertex, v_new[i1], v_new[i2]);
                volume.accumulate(signed_volume(t), centroid(t));
                #if defined(DEBUG_VOLUME)
                borderPoints << t;
                #endif
            }
        }
    }
    volume.normalise();
    return volume;
}

template class vtb::core::Flooding_solver<VTB_REAL_TYPE>;