Geometric.hpp

Go to the documentation of this file.

#pragma once
#include "DNDS/Defines.hpp"
 
namespace DNDS::Geom
{
    using t_index = int32_t;
    const t_index invalid_index = INT32_MAX;
    using t_real = double;
    using tPoint = Eigen::Vector3d;
    using tJacobi = Eigen::Matrix3d;
    using tGPoint = Eigen::Matrix3d;
    using tPointMap = Eigen::Map<tPoint, Eigen::Unaligned>;
    using tPointConstMap = Eigen::Map<const tPoint, Eigen::Unaligned>;
    using tGPointMap = Eigen::Map<tGPoint, Eigen::Unaligned>;
    using tGPointConstMap = Eigen::Map<const tGPoint, Eigen::Unaligned>;
 
    struct tPointPortable
    {
        std::array<real, 3> d;
        DNDS_DEVICE_TRIVIAL_COPY_DEFINE(tPointPortable, tPointPortable)
        DNDS_DEVICE_CALLABLE auto map()
        {
            return tPointMap{d.data()};
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] auto map() const
        {
            return tPointConstMap{d.data()};
        }
    };
 
    struct tGPointPortable
    {
        std::array<real, 9> d;
        DNDS_DEVICE_TRIVIAL_COPY_DEFINE(tGPointPortable, tGPointPortable)
        DNDS_DEVICE_CALLABLE auto map()
        {
            return tGPointMap{d.data()};
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] auto map() const
        {
            return tGPointConstMap{d.data()};
        }
    };
 
    static_assert(std::is_trivially_copyable_v<tPointPortable>);
    static_assert(std::is_trivially_copyable_v<tGPointPortable>);
 
    using tSmallCoords = Eigen::Matrix<real, 3, Eigen::Dynamic>;
    struct SmallCoordsAsVector : public tSmallCoords
    {
        auto operator[](Eigen::Index i)
        {
            return tSmallCoords::operator()(EigenAll, i);
        }
 
        auto operator[](Eigen::Index i) const
        {
            return tSmallCoords::operator()(EigenAll, i);
        }
    };
 
    static_assert(std::is_signed_v<t_index>);
 
    using tLocalMatStruct = std::vector<std::vector<index>>;
}
 
namespace DNDS::Geom
{
 
    inline std::vector<real> JacobiToSTDVector(const tJacobi &j)
    {
        std::vector<real> ret(9);
        for (int i = 0; i < 9; i++)
            ret[i] = j(i);
        return ret; // TODO this is bullshit code
    }
 
    inline tJacobi STDVectorToJacobi(const std::vector<real> &v)
    {
        tJacobi ret;
        DNDS_assert(v.size() == 9);
        for (int i = 0; i < 9; i++)
            ret(i) = v[i];
        return ret; //  TODO this is bullshit code
    }
 
    inline std::vector<real> VectorToSTDVector(const Eigen::VectorXd &v)
    {
        return {v.begin(), v.end()}; //  TODO this is bullshit code
    }
 
    inline Eigen::VectorXd STDVectorToVector(const std::vector<real> &v)
    {
        auto ret = Eigen::VectorXd(v.size());
        for (size_t i = 0; i < v.size(); i++)
            ret[i] = v[i]; //  TODO this is bullshit code
        return ret;
    }
}
 
namespace DNDS::Geom
{
    /**
     * @brief get normal vector from facial jacobian
     * note that
     * $$
     * J_{ij} = \partial_{\xi_j} x_i
     * $$
     */
    template <int d>
    inline tPoint FacialJacobianToNormVec(const tJacobi &J)
    {
        static_assert(d == 2 || d == 3);
        if constexpr (d == 2)
        {
            DNDS_assert_info(J(2, 0) == 0, "Must be a line in x-y plane");
            DNDS_assert_info(J(EigenAll, 1).norm() == 0, "Must be a line in x-y plane");
            DNDS_assert_info(J(EigenAll, 2).norm() == 0, "Must be a line in x-y plane");
            return tPoint{J(1, 0), -J(0, 0), 0};
        }
        if constexpr (d == 3)
        {
            DNDS_assert_info(J(EigenAll, 2).norm() == 0, "Must be a face");
            return J(EigenAll, 0).cross(J(EigenAll, 1));
        }
    }
 
    /**
     * @brief input uNorm should already be normalized
     *
     */
    template <int dim>
    inline Eigen::Matrix<real, dim, dim> NormBuildLocalBaseV(const Eigen::Ref<const Eigen::Vector<real, dim>> &uNorm)
    {
        Eigen::Matrix<real, dim, dim> base;
        static_assert(dim == 2 || dim == 3, "dim is bad");
        if constexpr (dim == 3)
        {
            base({0, 1, 2}, 0) = uNorm;
            if (std::abs(uNorm(2)) < 0.9)
                base({0, 1, 2}, 1) = -uNorm.cross(tPoint{0, 0, 1}).normalized();
            else
                base({0, 1, 2}, 1) = -uNorm.cross(tPoint{0, 1, 0}).normalized();
 
            base({0, 1, 2}, 2) = base({0, 1, 2}, 0).cross(base({0, 1, 2}, 1)).normalized();
        }
        else
        {
            base({0, 1}, 0) = uNorm;
            base(0, 1) = -uNorm(1); // x2=-y1
            base(1, 1) = uNorm(0);  // y2=x1
        }
        return base;
    }
 
    inline tGPoint RotZ(real theta)
    {
        theta *= pi / 180.0;
        return tGPoint{{cos(theta), -sin(theta), 0},
                       {sin(theta), cos(theta), 0},
                       {0, 0, 1}};
    }
 
    inline tGPoint RotX(real theta)
    {
        theta *= pi / 180.0;
        return tGPoint{{1, 0, 0},
                       {0, cos(theta), -sin(theta)},
                       {0, sin(theta), cos(theta)}};
    }
 
    inline tGPoint RotY(real theta)
    {
        theta *= pi / 180.0;
        return tGPoint{
            {cos(theta), 0, sin(theta)},
            {0, 1, 0},
            {-sin(theta), 0, cos(theta)}};
    }
 
    inline tGPoint RotateAxis(const tPoint &axis)
    {
        real theta = axis.norm();
        if (theta < smallReal)
            return tGPoint::Identity();
        tPoint axn = axis / theta;
        tGPoint Omega{{0, -axn(2), axn(1)},
                      {axn(2), 0, -axn(0)},
                      {-axn(1), axn(0), 0}};
        return axn * axn.transpose() + std::sin(theta) * Omega - std::cos(theta) * (Omega * Omega);
    }
 
    inline tGPoint CrossVecToMat(const tPoint &axn)
    {
        return tGPoint{{0, -axn(2), axn(1)},
                       {axn(2), 0, -axn(0)},
                       {-axn(1), axn(0), 0}};
    }
 
    inline tPoint CrossMatToVec(const tGPoint &axnM)
    {
        return tPoint{axnM(2, 1), axnM(0, 2), axnM(1, 0)};
    }
 
    template <int dim>
    tPoint ToThreeDim(const Eigen::Vector<real, dim> &v)
    {
        tPoint r;
        r.setZero();
        r(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>)) = v;
        return r;
    }
}