wetted_surface_solver_openmp.cc

#include <vtestbed/config/openmp.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/bisection.hh>
#include <vtestbed/core/grid_interpolation.hh>
#include <vtestbed/core/policy.hh>
#include <vtestbed/core/ship.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/core/wetted_surface_solver.hh>
#include <vtestbed/geometry/basis.hh>
#include <vtestbed/geometry/tetrahedron.hh>
#include <vtestbed/geometry/triangle.hh>
#include <vtestbed/profile/profile.hh>

//#define VTB_DEBUG_TRIANGLES

#if defined(VTB_DEBUG_TRIANGLES)
#include <iostream>
#include <vtestbed/core/debug.hh>
#endif

namespace vtb {

    namespace core {

        template <class T, int N>
        class Wetted_surface_solver_openmp: public Wetted_surface_solver<T,N> {

        private:
            using base_type = Wetted_surface_solver<T,N>;
            using typename base_type::ship_type;
            using typename base_type::grid_type;
            using typename base_type::panel_type;
            using typename base_type::vertex_type;

        private:
            std::vector<vtb::geometry::Vertex<T,N>> _efvertices;

        public:

            void
            solve(
                const ship_type& ship,
                const grid_type& grid_txyz,
                Array3<T> wavy_surface
            ) override;

            void
            solve(
                const ship_type& ship,
                const grid_type& grid_txyz,
                Array3<vertex_type> wavy_surface
            ) override;

        private:

            template <class Interpolation> void
            do_solve(
                const ship_type& ship,
                const grid_type& grid_txyz,
                Interpolation interpolate
            );

        };

    }

}

template <class T, int N>
void
vtb::core::Wetted_surface_solver_openmp<T,N>::solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    Array3<T> wavy_surface
) {
    using function_type = Linear_function<T,2>;
    using interpolation_type = Grid_interpolation<T,2,function_type>;
    const auto _ = blitz::Range::all();
    Grid<T,2> grid_xy{grid_txyz.select(1,2)};
    function_type function(wavy_surface(wavy_surface.extent(0)-1,_,_), grid_xy);
    this->do_solve(ship, grid_txyz, interpolation_type(grid_xy,std::move(function)));
}

template <class T, int N>
void
vtb::core::Wetted_surface_solver_openmp<T,N>::solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    Array3<vertex_type> wavy_surface
) {
    using function_type = Barycentric_function<T>;
    using interpolation_type = Grid_interpolation<T,2,function_type>;
    const auto _ = blitz::Range::all();
    Grid<T,2> grid_xy{grid_txyz.select(1,2)};
    function_type function(wavy_surface(wavy_surface.extent(0)-1,_,_));
    this->do_solve(ship, grid_txyz, interpolation_type(grid_xy,std::move(function)));
}

template <class T, int N>
template <class Interpolation>
void
vtb::core::Wetted_surface_solver_openmp<T,N>::do_solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    Interpolation interpolate
) {
    using vec2 = vtb::core::Vector2<T>;
    using vec3 = vtb::core::Vector3<T>;
    using bool3 = blitz::TinyVector<bool,3>;
    using vtb::geometry::Line_segment;
    using vtb::geometry::Tetrahedron;
    using vtb::geometry::Triangle;
    using vtb::geometry::unit;
    const auto& hull = ship.hull();
    const auto& bf_vertices = hull.vertices();
    const int nvertices = bf_vertices.size();
    auto& vertices = this->_efvertices;
    vertices.resize(nvertices);
    // Transform all hull vertices to Earth-fixed coordinate system.
    #if defined(VTB_OPENMP)
    #pragma omp parallel for
    #endif
    for (int i=0; i<nvertices; ++i) {
        vertices[i] = ship.rotation_matrix()*bf_vertices[i] + ship.position();
    }
    Grid<T,2> grid_xy{grid_txyz.select(1,2)};
    const auto& faces = hull.faces();
    const int npanels = faces.size();
    auto& mask = this->_wetted_vertices_mask;
    this->_all_panels.resize(npanels);
    auto& wetted_panels = this->_wetted_panels;
    wetted_panels.clear();
    wetted_panels.reserve(npanels/2);
    auto& dry_panels = this->_dry_panels;
    dry_panels.clear();
    dry_panels.reserve(npanels/2);
    mask.resize(nvertices);
    // set mask for visualisation
    mask_array valid_vertex(mask.size());
    for (int i=0; i<nvertices; ++i) {
        const auto& v = vertices[i];
        valid_vertex[i] = grid_xy.contains(vec2{v(0),v(1)});
    }
    mask.assign(nvertices, false);
    // clear waterline
    auto& waterline = this->_waterline;
    waterline.clear();
    Bisection<T,1> bisection;
    bisection.max_iterations(100);
    bisection.argument_precision(T{1e-3});
    bisection.function_precision(T{1e-3});
    auto intersection_point =
        [&] (vec3 from, vec3 to) -> vec3 {
            const vec3 delta = to - from;
            const T alpha = bisection(
                T{0}, T{1},
                [&] (T alpha) -> T {
                    vec3 x = from + alpha*delta;
                    return x(2) - interpolate({x(0),x(1)});
                }
            );
            return from + alpha*delta;
        };
    #if defined(VTB_OPENMP)
    #pragma omp parallel for schedule(dynamic,1)
    #endif
    for (int i=0; i<npanels; ++i) {
        const auto& face = faces[i];
        const int i0 = face[0];
        const int i1 = face[1];
        const int i2 = face[2];
        if (!valid_vertex[i0] || !valid_vertex[i1] || !valid_vertex[i2]) {
            continue;
        }
        Triangle<T,3> tr{vertices[i0], vertices[i1], vertices[i2]};
        vec3 centre = centroid(tr);
        T z_water[3] = {
            interpolate({tr[0][0], tr[0][1]}),
            interpolate({tr[1][0], tr[1][1]}),
            interpolate({tr[2][0], tr[2][1]})
        };
        bool3 above{
            tr[0][2] > z_water[0],
            tr[1][2] > z_water[1],
            tr[2][2] > z_water[2]
        };
        vec3 n = cross(vec3(tr[0]-tr[1]), vec3(tr[2]-tr[0]));
        T area = T{0.5}*length(n);
        // ignore degenerate panels
        if (!(area > 0)) { continue; }
        vec3 origin = ship.position();
        origin(2) = 0; // TODO ???
        T volume = signed_volume(Tetrahedron<T,3>(origin,tr));
        Line_segment<T,3> waterline_segment;
        bool has_waterline_segment = false;
        bool wet = !all(above);
        // if we have at least one vertex above the water level
        // and at least one vertex below the water level
        if (any(above) && wet) {
            vec3 a, b;
            bool x01 = above(0) ^ above(1);
            bool x02 = above(0) ^ above(2);
            bool x12 = above(1) ^ above(2);
            int outstanding = 0;
            if (x01 && x02) { outstanding = 0, a = tr[1], b = tr[2]; }
            else if (x01 && x12) { outstanding = 1, a = tr[0], b = tr[2]; }
            else { outstanding = 2, a = tr[0], b = tr[1]; }
            const vec3& outstanding_vertex = tr[outstanding];
            Triangle<T,3> subtr{
                outstanding_vertex,
                intersection_point(outstanding_vertex, a),
                intersection_point(outstanding_vertex, b)
            };
            T subarea = vtb::geometry::area(subtr);
            if (!(subarea > T{0})) { continue; }
            waterline_segment[0] = subtr[1];
            waterline_segment[1] = subtr[2];
            has_waterline_segment = true;
            #if defined(VTB_DEBUG_TRIANGLES)
            std::string tag;
            #endif
            // TODO why ship.position()
            T new_volume = signed_volume(Tetrahedron<T,3>(ship.position(),subtr));
            if (outstanding_vertex(2) > z_water[outstanding]) {
                area -= subarea;
                volume -= new_volume;
                centre = T{0.25} * (subtr[1] + subtr[2] + a + b);
                #if defined(VTB_DEBUG_TRIANGLES)
                tag = "polygon";
                #endif
            } else {
                area = subarea;
                volume = new_volume;
                centre = centroid(subtr);
                #if defined(VTB_DEBUG_TRIANGLES)
                tag = "triangle";
                #endif
            }
            #if defined(VTB_DEBUG_TRIANGLES)
            const auto& grid_xyz = grid_txyz.select(1,2,3);
            if (!grid_xyz.contains(centre)) {
                std::clog << "grid_xyz=" << grid_xyz << std::endl;
                std::clog << "centre=" << centre << std::endl;
                std::clog << "tr=" << tr << std::endl;
                std::clog << "subtr=" << subtr << std::endl;
                std::clog << "tag=" << tag << std::endl;
                std::clog << "area=" << area << std::endl;
                std::clog << "subarea=" << subarea << std::endl;
                dbg::write_triangle("tr", tr);
                dbg::write_triangle("subtr", subtr);
            }
            #endif
        }
        auto& panel = this->_all_panels[i];
        panel.centre(centre);
        panel.normal(-unit(n));
        panel.area(area);
        panel.volume(volume);
        panel.wetted(wet);
        panel.waterline(has_waterline_segment);
        if (wet) {
            #if defined(VTB_WITH_OPENMP)
            #pragma omp critical
            #endif
            {
                wetted_panels.emplace_back(panel);
                mask[i0] = true, mask[i1] = true, mask[i2] = true;
                if (has_waterline_segment) {
                    waterline.emplace_back(waterline_segment);
                }
            }
        } else {
            #if defined(VTB_OPENMP)
            #pragma omp critical
            #endif
            {
                dry_panels.emplace_back(panel);
            }
        }
    }
}

template <>
auto
vtb::core::make_wetted_surface_solver<VTB_REAL_TYPE,3,vtb::core::Policy::OpenMP>() ->
std::unique_ptr<Wetted_surface_solver<VTB_REAL_TYPE,3>> {
    return std::unique_ptr<Wetted_surface_solver<VTB_REAL_TYPE,3>>(
        new Wetted_surface_solver_openmp<VTB_REAL_TYPE,3>
    );
}

template class vtb::core::Wetted_surface_solver_openmp<VTB_REAL_TYPE,3>;