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) const
        {
        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 _idx ,
               std::initializer_list<int> _i )
        : element<N,NPI>(_p_node,_idx,_i),surf(0),dMs(0)
        {
        if (existNodes())
            {
            surf = calc_surf();
            n = calc_norm();
            for(int i=0;i<NPI;i++)
                { weight[i] = 2.0 * surf * Triangle::pds[i]; }
            }
        else
            { std::cerr<<"Error: Tri constructor has out of bound index\n"; }
        }
 
    double surf;
 
    double dMs;
 
    Eigen::Vector3d n;
 
    void interpolation(const 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(const 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 ,
            const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,NPI>> u ) const;
 
    Eigen::Matrix<double,NPI,1>  charges(const double &dMs ,
            const std::function<Eigen::Vector3d(const Nodes::Node&)> &getter ) const override;
 
    void correctionCharges(const std::function<Eigen::Vector3d(Nodes::Node)>& getter ,
                           Eigen::Matrix<double,Triangle::NPI,1> &localCharges ,
                           std::vector<double> &corr ) const;
 
    double demagEnergy(Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u ,
                       const Eigen::Ref<const Eigen::Matrix<double,NPI,1>> phi ) const;
 
    double potential(const std::function<Eigen::Vector3d(Nodes::Node)>& getter, int i) const;
 
    bool operator==(const Tri &f) const
        {
        /* As we rule out degenerate triangles, it is sufficient to test whether all the vertices
         * of f are also vertices of this. */
        return std::find(ind.begin(), ind.end(), f.ind[0]) != ind.end()
            && std::find(ind.begin(), ind.end(), f.ind[1]) != ind.end()
            && std::find(ind.begin(), ind.end(), f.ind[2]) != ind.end();
        }
 
    inline Eigen::Vector3d calc_norm(void) const
        {
        Eigen::Vector3d _n = normal_vect();
        _n.normalize();
        return _n;
        }
 
    Eigen::Matrix<double,Nodes::DIM,NPI> getPtGauss(void) const override
        {
        Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
        for (int i = 0; i < N; i++)
            {
            vec_nod.col(i) << getNode(i).p;
            }
        return vec_nod*eigen_a;
        }
 
    inline double calc_surf(void) const { return 0.5 * normal_vect().norm(); }
 
private:
    void orientate(void) override {}
 
    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 */