tetra.h

Go to the documentation of this file.

#ifndef tetra_h
#define tetra_h
 
#include <set>
 
#include "facette.h"
#include "node.h"
#include "time_integration.h"
#include "element.h"
 
namespace Tetra
    {
const double epsilon = EPSILON;
 
const int N = 4;   
#if ONE_GAUSS_POINT
    const int NPI = 1; 
    constexpr double A = 1. / 4.;           
    constexpr double u[NPI] = {A};          
    constexpr double v[NPI] = {A};          
    constexpr double w[NPI] = {A};          
    constexpr double pds[NPI] = {1./6.};    
    constexpr double a[N][NPI] = {{1. - u[0] - v[0] - w[0]},
                              {u[0]},
                              {v[0]},
                              {w[0]}};
 
    const Eigen::Matrix<double,N,NPI> eigen_a = (Eigen::MatrixXd(N,NPI) << a[0][0],
                                                                       a[1][0],
                                                                       a[2][0],
                                                                       a[3][0] ).finished();
#else
    const int NPI = 5; 
    constexpr double A = 1. / 4.;                
    constexpr double B = 1. / 6.;                
    constexpr double C = 1. / 2.;                
    constexpr double D = -2. / 15.;              
    constexpr double E = 3. / 40.;               
    constexpr double u[NPI] = {A, B, B, B, C};   
    constexpr double v[NPI] = {A, B, B, C, B};   
    constexpr double w[NPI] = {A, B, C, B, B};   
    constexpr double pds[NPI] = {D, E, E, E, E}; 
    constexpr double a[N][NPI] = {{1. - u[0] - v[0] - w[0], 1. - u[1] - v[1] - w[1],
                                   1. - u[2] - v[2] - w[2], 1. - u[3] - v[3] - w[3],
                                   1. - u[4] - v[4] - w[4]},
                                  {u[0], u[1], u[2], u[3], u[4]},
                                  {v[0], v[1], v[2], v[3], v[4]},
                                  {w[0], w[1], w[2], w[3], w[4]}};
 
    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[0][4],
                                       a[1][0], a[1][1], a[1][2], a[1][3], a[1][4],
                                       a[2][0], a[2][1], a[2][2], a[2][3], a[2][4],
                                       a[3][0], a[3][1], a[3][2], a[3][3], a[3][4] ).finished();
#endif
 
struct prm
    {
    std::string regName; 
    double alpha_LLG;    
    double A;            
    double Ms;           
    double K;            
    Eigen::Vector3d uk;  
    double K3;           
    Eigen::Vector3d ex;  
    Eigen::Vector3d ey;  
    Eigen::Vector3d ez;  
    double P;            
    double N0;           
    double sigma;        
    double lsd;          
    double lsf;          
    double spinHall;     
    double volume = 0; 
    void infos(void);  
    };
 
class Tet : public element<N,NPI>
    {
public:
    inline Tet(const std::vector<Nodes::Node> &_p_node ,
               const int _idx ,
               std::initializer_list<int> _i )
        : element<N,NPI>(_p_node,_idx,_i), idx(0)
        {
        zeroBasing();
        da.setZero();
 
        if (existNodes())
            {
            orientate();
 
            Eigen::Matrix3d J;
            // we have to rebuild the jacobian in case of ill oriented tetrahedron
            double detJ = Jacobian(J);
 
            if (fabs(detJ) < Tetra::epsilon)
                {
                std::cerr << "Singular jacobian in tetrahedron: |det(J)|= "
                        << fabs(detJ) << std::endl;
                element::infos();
                SYSTEM_ERROR;
                }
            Eigen::Matrix<double,N,Nodes::DIM> dadu; // Shape function derivatives
                                                     // (constant for linear tetrahedron)
            dadu << -1., -1., -1., 1., 0., 0., 0., 1., 0., 0., 0., 1.;
            da = dadu * J.inverse();
 
            for (int j = 0; j < NPI; j++)
                { weight[j] = detJ * Tetra::pds[j]; }// if NPI=1 weight[0] is the tetrahedron volume
            }
        // do nothing lambda's (usefull for spin transfer torque)
        extraField = [] ( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
        }
 
    Eigen::Matrix<double,N,Nodes::DIM> da;
 
    inline 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;
        }
 
    inline void interpolation(std::function<Eigen::Vector3d(Nodes::Node)> getter,
                              Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> result,
                              Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tx,
                              Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Ty,
                              Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Tz) 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;// interpolated value at Gauss point
 
        //Spatial derivatives: constant for linear elements, evaluated at any point including Gauss
        //point
        Tx = vec_nod * (da.col(Nodes::IDX_X)).replicate(1,NPI);
        Ty = vec_nod * (da.col(Nodes::IDX_Y)).replicate(1,NPI);
        Tz = vec_nod * (da.col(Nodes::IDX_Z)).replicate(1,NPI);
        }
 
    inline void interpolation_field(std::function<double(Nodes::Node)> getter,
                                    Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> X) const
        {
        Eigen::Matrix<double,N,1> scalar_nod;
        for (int i = 0; i < N; i++) scalar_nod(i) = getter(getNode(i));
        
        X.setZero();
        for (int j = 0; j < NPI; j++)
            {
            for (int i = 0; i < N; i++)
                {
                X.col(j) -= (scalar_nod[i] * da.row(i));
                }
            }
        }
 
    inline void interpolation(std::function<double(Nodes::Node, Nodes::index)> getter,
            Nodes::index idx, Eigen::Ref<Eigen::Matrix<double,Tetra::NPI,1>> result) const
        {
        Eigen::Matrix<double,N,1> scalar_nod;
 
        for (int i = 0; i < N; i++)
            {
            scalar_nod[i] = getter(getNode(i),idx);
            }
        result = scalar_nod.transpose() * eigen_a;
        }
    
    void lumping(Eigen::Ref<Eigen::Matrix<double,NPI,1>> alpha_eff, double prefactor,
                  Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE ) const;
 
    Eigen::Matrix<double,N,N> calcDiagBlock(const double c, Eigen::Matrix<double,N,1> &x) const;
 
    Eigen::Matrix<double,N,1> calcOffDiagBlock(const Nodes::index idx) const;
 
    void add_drift_BE(double alpha, double s_dt, double Vdrift,
                      Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
                      Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V, 
                      Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dUd_,
                      Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dVd_,
                      Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>> BE) const;
 
    Eigen::Matrix<double,NPI,1> calc_aniso_uniax(Eigen::Ref<const Eigen::Vector3d> uk,
            const double Kbis, const double s_dt,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
 
    Eigen::Matrix<double,NPI,1> calc_aniso_cub(Eigen::Ref<const Eigen::Vector3d> ex,
            Eigen::Ref<const Eigen::Vector3d> ey,
            Eigen::Ref<const Eigen::Vector3d> ez,
            const double K3bis, const double s_dt,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> V,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H_aniso) const;
 
    void integrales( Tetra::prm const &param, timing const &prm_t,
                    std::function< Eigen::Matrix<double,Nodes::DIM,NPI> (void)> calc_Hext,
                    Nodes::index idx_dir, double Vdrift);
 
    double exchangeEnergy(Tetra::prm const &param,
                          Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
                          Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
                          Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz) const;
 
    double uniaxialAnisotropyEnergy(Tetra::prm const &param,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
 
    double cubicAnisotropyEnergy(Tetra::prm const &param,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
 
    Eigen::Matrix<double,Tetra::NPI,1> charges(const double &Ms,
            std::function<Eigen::Vector3d(const Nodes::Node&)> getter) const;
 
    double demagEnergy(Tetra::prm const &param,
                       Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudx,
                       Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudy,
                       Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> dudz,
                       Eigen::Ref<Eigen::Matrix<double,NPI,1>> phi) const;
 
    double zeemanEnergy(Tetra::prm const &param  ,
            Eigen::Ref<Eigen::Vector3d> const Hext ,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> const u ) const;
 
    double zeemanEnergy(Tetra::prm const &param ,
                        double fieldAmp ,
                        std::vector<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> &spaceField,
                        Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> const u) const;
 
    double Jacobian(Eigen::Ref<Eigen::Matrix3d> J);
 
    double calc_vol(void) const;
 
    int idx;
 
    std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H)> extraField;
 
    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;
        }
 
private:
    void orientate(void)
        {
        if (calc_vol() < 0.0)
                std::swap(ind[2], ind[3]);
        }
 
    };  // end class Tetra
 
    Eigen::Matrix<double,NPI,1> calc_alpha_eff(const double dt, const double alpha,
                                               Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff);
 
    Eigen::Matrix<double,Nodes::DIM,NPI> calc_gradV(Tet const &tet,std::vector<double> &V);
 
    Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> calc_Hst(Tetra::Tet const &tet,
                                                         const double prefactor,
                                                         std::vector<Eigen::Vector3d> &s);
}   // end namespace Tetra
 
#endif /* tetra_h */