triangle.h

Go to the documentation of this file.

#ifndef triangle_h
#define triangle_h
 
#include "element.h"
 
namespace Triangle
    {
const int N = 3;
 
#if ONE_GAUSS_POINT  // Single Gauss point at barycenter of triangle
    const int NPI = 1;
 
    constexpr double u[NPI] = {1. / 3.};
 
    constexpr double v[NPI] = {1. / 3.};
 
    constexpr double pds[NPI] = {1. / 2.};
 
    constexpr double a[N][NPI] = { {1. - u[0] - v[0]}, {u[0]}, {v[0]}};
 
    const Eigen::Matrix<double,N,NPI> eigen_a =
        (Eigen::MatrixXd(N,NPI) << a[0][0], a[1][0], a[2][0] ).finished();
#else
    const int NPI = 4;
 
    constexpr double u[NPI] = {1 / 3., 1 / 5., 3 / 5., 1 / 5.};
 
    constexpr double v[NPI] = {1 / 3., 1 / 5., 1 / 5., 3 / 5.};
 
    constexpr double pds[NPI] = {-27 / 96., 25 / 96., 25 / 96., 25 / 96.};
 
    constexpr double a[N][NPI] = {
        {1. - u[0] - v[0], 1. - u[1] - v[1], 1. - u[2] - v[2], 1. - u[3] - v[3]},
        {u[0], u[1], u[2], u[3]},
        {v[0], v[1], v[2], v[3]}};
 
    const Eigen::Matrix<double,N,NPI> eigen_a =
            (Eigen::MatrixXd(N,NPI) << a[0][0], a[0][1], a[0][2], a[0][3],
                                       a[1][0], a[1][1], a[1][2], a[1][3],
                                       a[2][0], a[2][1], a[2][2], a[2][3] ).finished();
#endif
struct prm
    {
    std::string regName;   
    bool suppress_charges; 
    double Ks;             
    Eigen::Vector3d uk;    
    double V;              
    double jn;             
    Eigen::Vector3d uP;    
    Eigen::Vector3d s;     
    inline void infos(const bool spinAcc = false)
        {
        std::cout << "surface region name = " << regName
                  << " ; suppress charges = " << suppress_charges << std::endl;
 
        if (Ks != 0)
            {
            std::cout << "Ks*a = " << Ks << "*[ " << uk << "]" << std::endl;
            }
        else
            std::cout << "no surface anisotropy" << std::endl;
 
        if (spinAcc)
            {
            std::cout<< "V= " << V << "; jn= " << jn << "; uP= " << uP << "; s= " << s << std::endl;
            }
        };
    };
 
class Tri : public element<N,NPI>
    {
public:
    inline Tri(const std::vector<Nodes::Node> &_p_node ,
               const int _NOD ,
               const int _idx ,
               std::initializer_list<int> _i )
        : element<N,NPI>(_p_node,_idx,_i),dMs(0)
        {
        if (_NOD > 0)
            {
            zeroBasing();
            surf = calc_surf();
            n = calc_norm();
            }
        else
            {
            surf = 0.0;
            n = Eigen::Vector3d(0,0,0);
            }  // no index shift here if NOD == 0 : usefull while reordering triangle indices
 
        for(int i=0;i<NPI;i++)
            { weight[i] = 2.0 * surf * Triangle::pds[i]; }
        }
 
    double surf;
 
    double dMs;
 
    Eigen::Vector3d n;
 
    void interpolation(std::function<Eigen::Vector3d(Nodes::Node)> getter ,
                       Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result ) const
        {
        Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
        for (int i = 0; i < N; i++) vec_nod.col(i) = getter(getNode(i));
        
        result = vec_nod * eigen_a;
        }
 
    void interpolation(std::function<double(Nodes::Node)> getter ,
                       Eigen::Ref<Eigen::Matrix<double,NPI,1>> result ) const
        {
        Eigen::Matrix<double,N,1> scalar_nod;
        for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
        
        result = scalar_nod.transpose() * eigen_a;
        }
 
    void integrales(Triangle::prm const &params );
 
    double anisotropyEnergy(Triangle::prm const &param ,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> const u ) const;
 
    Eigen::Matrix<double,NPI,1>  charges(const double &dMs ,
            std::function<Eigen::Vector3d(const Nodes::Node&)> getter ) const;
 
    void correctionCharges(std::function<Eigen::Vector3d(Nodes::Node)> getter ,
                           Eigen::Matrix<double,Triangle::NPI,1> &localCharges ,
                           std::vector<double> &corr );
 
    double demagEnergy(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u ,
                       Eigen::Ref<Eigen::Matrix<double,NPI,1>> phi ) const;
 
    double potential(std::function<Eigen::Vector3d(Nodes::Node)> getter, int i) const;
 
    inline bool operator<(const Tri &f) const
        {
        return (this->ind[0] < f.ind[0])
               || ((this->ind[0] == f.ind[0])
                   && ((this->ind[1] < f.ind[1])
                       || ((this->ind[1] == f.ind[1]) && (this->ind[2] < f.ind[2]))));
        }
 
    bool operator==(const Tri &f) const
        {
        const int a= this->ind[0];
        const int b= this->ind[1];
        const int c= this->ind[2];
        const int i= f.ind[0];
        const int j= f.ind[1];
        const int k= f.ind[2];
        return ((a==i)&&(b==j)&&(c==k))||((a==j)&&(b==i)&&(c==k))
             ||((a==i)&&(b==k)&&(c==j))||((a==k)&&(b==i)&&(c==j))
             ||((a==j)&&(b==k)&&(c==i))||((a==k)&&(b==j)&&(c==i));
        }
 
    inline Eigen::Vector3d calc_norm(void) const
        {
        Eigen::Vector3d _n = normal_vect();
        _n.normalize();
        return _n;
        }
 
    void getPtGauss(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result) const
        {
        Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
        for (int i = 0; i < N; i++)
            {
            vec_nod.col(i) << getNode(i).p;
            }
        result = vec_nod * eigen_a;
        }
 
    inline double calc_surf(void) const { return 0.5 * normal_vect().norm(); }
 
private:
    void orientate(void) {}
 
    inline Eigen::Vector3d normal_vect() const
        {
        Eigen::Vector3d p0p1 = getNode(1).p - getNode(0).p;
        Eigen::Vector3d p0p2 = getNode(2).p - getNode(0).p;
 
        return p0p1.cross(p0p2);
        }
    };  // end class Tri
 
    }  // namespace Triangle
 
#endif /* triangle_h */