TrianglePN.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_TrianglePN_h__
#define __geometry_TrianglePN_h__

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Evaluate barycenter on triangle using Curved PN Triangles

    "Curved PN Triangles"
    2001 Alex Vlachos, Jorg Peters, Chas Boyd, Jason L. Mitchell

    @ingroup geometry
*//***************************************************************************/
template <typename T>
class TrianglePN
{
    public:
                /// @brief halfway between two of the vertices
        enum            Edge
                        {
                            e_v1v2,
                            e_v2v3,
                            e_v3v1
                        };

                        TrianglePN(void)                                    {}

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

        Vector<3,T>     midpoint(Edge a_edge);
        Vector<3,T>     midnormal(Edge a_edge);

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

    private:
        Vector<3,T> m_b300;
        Vector<3,T> m_b030;
        Vector<3,T> m_b003;
        Vector<3,T> m_b210;
        Vector<3,T> m_b120;
        Vector<3,T> m_b021;
        Vector<3,T> m_b012;
        Vector<3,T> m_b102;
        Vector<3,T> m_b201;
        Vector<3,T> m_b111;

        Vector<3,T> m_n200;
        Vector<3,T> m_n020;
        Vector<3,T> m_n002;
        Vector<3,T> m_n110;
        Vector<3,T> m_n011;
        Vector<3,T> m_n101;

};

template <typename T>
inline void TrianglePN<T>::configure(
    const Vector<3,T>& v1,
    const Vector<3,T>& v2,
    const Vector<3,T>& v3,
    const Vector<3,T>& n1,
    const Vector<3,T>& n2,
    const Vector<3,T>& n3)
{
    if(isZero(n1) || isZero(n2) || isZero(n3))
    {
#if FE_CODEGEN<=FE_DEBUG
        feLog("v1 %s\n",c_print(v1));
        feLog("v2 %s\n",c_print(v2));
        feLog("v3 %s\n",c_print(v3));
        feLog("n1 %s\n",c_print(n1));
        feLog("n2 %s\n",c_print(n2));
        feLog("n3 %s\n",c_print(n3));
#endif
        feX(e_unsolvable,"TrianglePN<T>::configure","zero length normal(s)");
    }

    m_b300=v1;
    m_b030=v2;
    m_b003=v3;

    const Vector<3,T> v2_v1=v2-v1;
    const Vector<3,T> v3_v2=v3-v2;
    const Vector<3,T> v1_v3=v1-v3;

    const T w12=dot(v2_v1,n1);
    const T w21=dot(-v2_v1,n2);
    const T w23=dot(v3_v2,n2);
    const T w32=dot(-v3_v2,n3);
    const T w31=dot(v1_v3,n3);
    const T w13=dot(-v1_v3,n1);

    m_b210=(T(2)*v1+v2-w12*n1)*T(1.0/3.0);
    m_b120=(T(2)*v2+v1-w21*n2)*T(1.0/3.0);
    m_b021=(T(2)*v2+v3-w23*n2)*T(1.0/3.0);
    m_b012=(T(2)*v3+v2-w32*n3)*T(1.0/3.0);
    m_b102=(T(2)*v3+v1-w31*n3)*T(1.0/3.0);
    m_b201=(T(2)*v1+v3-w13*n1)*T(1.0/3.0);

    const Vector<3,T> E=(m_b210+m_b120+m_b021+m_b012+m_b102+m_b201)*T(1.0/6.0);
    const Vector<3,T> V=(v1+v2+v3)*T(1.0/3.0);
    m_b111=E+(E-V)*T(1.0/2.0);

    m_n200=n1;
    m_n020=n2;
    m_n002=n3;

    const Vector<3,T> n1_n2=n1+n2;
    const Vector<3,T> n2_n3=n2+n3;
    const Vector<3,T> n3_n1=n3+n1;

    const T v12=T(2)*dot(v2_v1,n1_n2)/dot(v2_v1,v2_v1);
    const T v23=T(2)*dot(v3_v2,n2_n3)/dot(v3_v2,v3_v2);
    const T v31=T(2)*dot(v1_v3,n3_n1)/dot(v1_v3,v1_v3);

    const Vector<3,T> h110=n1_n2-v12*(v2_v1);
    const Vector<3,T> h011=n2_n3-v23*(v3_v2);
    const Vector<3,T> h101=n3_n1-v31*(v1_v3);

    m_n110=unit(h110);
    m_n011=unit(h011);
    m_n101=unit(h101);
}

template <typename T>
inline void TrianglePN<T>::solve(
    const Barycenter<T>& barycenter,
    Vector<3,T>& vertex,Vector<3,T>& normal) const
{
    const T u=barycenter[0];
    const T v=barycenter[1];
    const T w=1.0-u-v;

    const T uu=u*u;
    const T vv=v*v;
    const T ww=w*w;

    vertex=m_b300*uu*u+m_b030*vv*v+m_b003*ww*w+
            T(3)*(m_b210*uu*v+m_b120*u*vv+m_b201*uu*w+
            m_b021*vv*w+m_b102*u*ww+m_b012*v*ww)+
            T(6)*m_b111*u*v*w;

    normal=unit(m_n200*uu+m_n020*vv+m_n002*ww+
            m_n110*u*v+m_n011*v*w+m_n101*u*w);
}

template <typename T>
inline Vector<3,T> TrianglePN<T>::midpoint(Edge a_edge)
{
    switch(a_edge)
    {
        case e_v1v2:
            return 0.125*(m_b300+m_b030+T(3)*(m_b210+m_b120));

        case e_v2v3:
            return 0.125*(m_b030+m_b003+T(3)*(m_b021+m_b012));

        case e_v3v1:
            return 0.125*(m_b300+m_b003+T(3)*(m_b201+m_b102));
    }
    return Vector<3,T>(0.0,0.0,0.0);
}

template <typename T>
inline Vector<3,T> TrianglePN<T>::midnormal(Edge a_edge)
{
    switch(a_edge)
    {
        case e_v1v2:
            return unit(m_n200+m_n020+m_n110);

        case e_v2v3:
            return unit(m_n020+m_n002+m_n011);

        case e_v3v1:
            return unit(m_n200+m_n002+m_n101);
    }
    return Vector<3,T>(0.0,0.0,1.0);
}

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

#endif /* __geometry_TrianglePN_h__ */