TrianglePower.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __geometry_TrianglePower_h__
#define __geometry_TrianglePower_h__

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Evaluate barycenter on triangle of surface using MatrixPower

    @ingroup geometry
*//***************************************************************************/
template <typename T>
class TrianglePower
{
    public:
                TrianglePower(void)                                         {}

        void    configure(const Vector<3,T>& v0,
                        const Vector<3,T>& v1,
                        const Vector<3,T>& v2,
                        const Vector<3,T>& n0,
                        const Vector<3,T>& n1,
                        const Vector<3,T>& n2);

        void    solve(const Barycenter<T>& barycenter,
                        Vector<3,T>& v,Vector<3,T>& n) const;

        void    setIterations(U32 iterations)
                {   m_matrixPower.setIterations(iterations); }

    private:
        void    calcFrames(
                        Matrix<3,4,T>& a_rTransform1,
                        Matrix<3,4,T>& a_rTransform2,
                        const Vector<3,T>& a_rDirection,
                        const Vector<3,T>& a_rVertex1,
                        const Vector<3,T>& a_rNormal1,
                        const Vector<3,T>& a_rVertex2,
                        const Vector<3,T>& a_rNormal2) const;

        MatrixPower< Matrix<3,4,T> >    m_matrixPower;

        Matrix<3,4,T>                   m_frame0;
        Matrix<3,4,T>                   m_delta01;
        Vector<3,T>                     m_direction;
        Vector<3,T>                     m_v2;
        Vector<3,T>                     m_n2;
};

template <typename T>
inline void TrianglePower<T>::configure(
    const Vector<3,T>& v0,
    const Vector<3,T>& v1,
    const Vector<3,T>& v2,
    const Vector<3,T>& n0,
    const Vector<3,T>& n1,
    const Vector<3,T>& n2)
{
    m_direction=unit(v1-v0);

    Matrix<3,4,T> frame1;
    calcFrames(m_frame0,frame1,m_direction,v0,n0,v1,n1);
    Matrix<3,4,T> inv0;
    invert(inv0,m_frame0);
    m_delta01=inv0*frame1;

    m_v2=v2;
    m_n2=n2;
}

template <typename T>
inline void TrianglePower<T>::solve(
    const Barycenter<T>& barycenter,
    Vector<3,T>& vertex,Vector<3,T>& normal) const
{
    const T f01=(barycenter[0]+barycenter[1] > T(0))?
            barycenter[1]/(barycenter[0]+barycenter[1]): T(0);
    const T f2=T(1)-barycenter[0]-barycenter[1];

    Matrix<3,4,T> partial;
    m_matrixPower.solve(partial,m_delta01,f01);
    Matrix<3,4,T> mid=m_frame0*partial;
    const Vector<3,T> mid01=mid.translation();
    const Vector<3,T> norm01=mid.up();

    FEASSERT(magnitude(norm01)>0.99);
    FEASSERT(magnitude(norm01)<1.01);

    Matrix<3,4,T> frame2;
    Matrix<3,4,T> frame3;
    calcFrames(frame2,frame3,m_direction,mid01,norm01,m_v2,m_n2);
    Matrix<3,4,T> inv2;
    invert(inv2,frame2);
    Matrix<3,4,T> delta2=inv2*frame3;
    m_matrixPower.solve(partial,delta2,f2);
    mid=frame2*partial;
    vertex=mid.translation();
    normal=mid.up();

    FEASSERT(magnitude(normal)>0.99);
    FEASSERT(magnitude(normal)<1.01);
}

template <typename T>
inline void TrianglePower<T>::calcFrames(
    Matrix<3,4,T>& a_rTransform1,Matrix<3,4,T>& a_rTransform2,
    const Vector<3,T>& a_rDirection,
    const Vector<3,T>& a_rVertex1,const Vector<3,T>& a_rNormal1,
    const Vector<3,T>& a_rVertex2,const Vector<3,T>& a_rNormal2) const
{
    a_rTransform1.left()=unit(cross(a_rNormal1,a_rDirection));
    a_rTransform2.left()=unit(cross(a_rNormal2,a_rDirection));
    a_rTransform1.direction()=cross(a_rTransform1.left(),a_rNormal1);
    a_rTransform2.direction()=cross(a_rTransform2.left(),a_rNormal2);
    a_rTransform1.up()=a_rNormal1;
    a_rTransform2.up()=a_rNormal2;
    a_rTransform1.translation()=a_rVertex1;
    a_rTransform2.translation()=a_rVertex2;

#if FALSE
    feLog("a_rDirection\n%s\n",c_print(a_rDirection));
    feLog("frame1\n%s\n",c_print(a_rTransform1));
    feLog("frame2\n%s\n",c_print(a_rTransform2));
#endif
}

} /* namespace ext */
} /* namespace fe */

#endif /* __geometry_TrianglePower_h__ */