triangle.hpp

Go to the documentation of this file.

// OpenBEM - Copyright (C) 2026 Shashwat Sharma
 
// This file is part of OpenBEM.
 
// OpenBEM is free software: you can redistribute it and/or modify it under the terms of the
// GNU General Public License as published by the Free Software Foundation, either version 3
// of the License, or (at your option) any later version.
 
// You should have received a copy of the GNU General Public License along with OpenBEM.
// If not, see <https://www.gnu.org/licenses/>.
 
 
#ifndef GEOM_TRIANGLE_H
#define GEOM_TRIANGLE_H
 
#include <stdexcept>
#include <array>
#include <cassert>
 
#include "types.hpp"
#include "constants.hpp"
 
 
namespace bem
{
 
const Float TRIANGLE_DEFAULT_TOL = 1.0e-6;
 
template <uint8_t dim>
class Triangle
{
public:
 
    Triangle()
    {
        set_v(reference_triangle().v());
        return;
    };
 
 
    Triangle(
        ConstEigRef<EigMatMN<Float, dim, 3>> v,
        EigRowVecN<Float, 3> edge_polarities = EigRowVecN<Float, 3>::Constant(1, 3, 1)
        )
    {
        set_v(v, edge_polarities);
        return;
    };
 
 
    Triangle(
        std::array<EigColVecN<Float, dim>, 3> v_array,
        EigRowVecN<Float, 3> edge_polarities = EigRowVecN<Float, 3>::Constant(1, 3, 1)
        )
    {
        EigMatMN<Float, dim, 3> v;
        for (uint8_t ii = 0; ii < 3; ++ii)
            v.col(ii) = v_array[ii];
        set_v(v, edge_polarities);
        return;
    };
 
 
    void set_v(
        ConstEigRef<EigMatMN<Float, dim, 3>> v,
        EigRowVecN<Float, 3> edge_polarities = EigRowVecN<Float, 3>::Constant(1, 3, 1)
        );
 
 
    const EigMatMN<Float, dim, 3>& v() const
    {
        return v_;
    };
 
 
    EigColVecN<Float, dim> v(const uint8_t idx) const
    {
        if (idx >= 3)
            throw std::out_of_range("Triangle::v(): vertex index out of bounds.");
        return v_.col(idx);
    };
 
 
    const EigRowVecN<Float, 3>& edge_polarities() const
    {
        return edge_polarities_;
    };
 
 
    Float edge_polarities(const uint8_t idx) const
    {
        if (idx >= 3)
            throw std::out_of_range("Triangle::edge_polarities(): edge index out of bounds.");
        return edge_polarities_[idx];
    };
 
 
    void reverse()
    {
        set_v(v_.rowwise().reverse(), edge_polarities_.reverse());
        return;
    };
 
 
    Float area() const
    { return area_; };
 
 
    EigColVecN<Float, 3> normal() const
    { return normal_; };
 
 
    EigColVecN<Float, dim> centroid() const
    { return centroid_; };
 
 
    Float mean_edge_length() const
    {
        Float len = 0.0;
        for (uint8_t ii = 0; ii < 3; ii++)
            len += (v_.col(ii) - v_.col((ii + 1) % 3)).norm();
        return len / (Float)3;
    };
 
 
    Triangle<2> to_2d() const;
 
 
    Triangle<3> to_3d(const Float z = 0.0) const;
 
 
    EigColVecN<Float, dim> local_origin() const
    {
        return v_.col(0);
    };
 
 
    EigMatMN<Float, dim, dim> local_coordinate_basis() const;
 
 
    void get_plane_projection(
        EigMatNX<Float, dim>& r_proj,
        EigRowVec<Float>& d,
        ConstEigRef<EigMatNX<Float, 3>> r,
        uint8_t ref_idx = 0
        ) const;
 
 
    static Triangle<dim> reference_triangle();
 
 
    EigColVecN<Float, 3> edge_lengths() const
    {
        EigColVecN<Float, 3> lengths;
        for (uint8_t ii = 0; ii < 3; ii++)
            lengths[ii] = (v_.col((ii + 1) % 3) - v_.col((ii) % 3)).norm();
        return lengths;
    };
 
 
    uint8_t longest_edge_index() const
    {
        uint8_t result;
        edge_lengths().maxCoeff(&result);
        assert(result < 3 && "Triangle::longest_edge(): Result out of range.");
        return result;
    };
 
 
    uint8_t shortest_edge_index() const
    {
        uint8_t result;
        edge_lengths().minCoeff(&result);
        assert(result < 3 && "Triangle::shortest_edge(): Result out of range.");
        return result;
    };
 
 
    Float longest_edge_length() const
    {
        return edge_lengths().maxCoeff();
    };
 
 
    Float shortest_edge_length() const
    {
        return edge_lengths().minCoeff();
    };
 
 
    EigMatNX<Float, 3> barycentric_coords(ConstEigRef<EigMatNX<Float, dim>> p) const;
 
 
    EigRowVec<uint8_t> projection_loc(ConstEigRef<EigMatNX<Float, 3>> r) const;
 
 
    bool point_in_triangle(ConstEigRef<EigColVecN<Float, dim>> r, const Float tol = TRIANGLE_DEFAULT_TOL) const;
 
 
    template <uint8_t dim_out>
    EigColVecN<Float, dim_out> map_points(
        ConstEigRef<EigColVecN<Float, dim>> points,
        const Triangle<dim_out>& tri_out
        ) const
    {
        EigMatNX<Float, 3> lambda = barycentric_coords(points);
        EigColVecN<Float, dim_out> points_out = tri_out.v(0) * lambda.row(0) +
                                      tri_out.v(1) * lambda.row(1) +
                                      tri_out.v(2) * lambda.row(2);
        return points_out;
    }
 
 
    static Float area(ConstEigRef<EigMatMN<Float, dim, 3>> v);
 
 
    static EigColVecN<Float, 3> normal(ConstEigRef<EigMatMN<Float, dim, 3>> v);
 
 
    static EigColVecN<Float, dim> centroid(ConstEigRef<EigMatMN<Float, dim, 3>> v)
    { return v.rowwise().sum() / (Float)3; };
 
 
protected:
 
    EigMatMN<Float, dim, 3> v_;
    Float area_;
    EigColVecN<Float, 3> normal_;
    EigColVecN<Float, dim> centroid_;
    EigRowVecN<Float, 3> edge_polarities_ = EigRowVecN<Float, 3>::Constant(1, 3, 1);
 
};
 
}
 
#ifndef BEM_LINKED
#include "geometry/primitives/triangle.cpp"
#endif
 
#endif