tetra.h

Go to the documentation of this file.

#ifndef tetra_h
#define tetra_h
 
#include <set>
 
#include "absorbing_ends.h"
#include "triangle.h"
#include "node.h"
#include "time_integration.h"
#include "element.h"
 
namespace Tetra
    {
constexpr double epsilon = EPSILON;
 
constexpr int N = 4;   
#if ONE_GAUSS_POINT
    constexpr 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
    constexpr 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 DMI;          
    double Ms;           
    double K;            
    Eigen::Vector3d uk;  
    double K3;           
    Eigen::Matrix3d e;   
    double P;            
    double N0;           
    double sigma;        
    double lsd;          
    double lsf;          
    double volume = 0; 
    bool absorbing = false;
 
    void infos(void);  
    };
 
class Tet : public element<N,NPI>
    {
public:
    inline Tet(const std::vector<Nodes::Node> &_p_node ,
               const int _idx ,
               const std::initializer_list<int> _i )
        : element<N,NPI>(_p_node,_idx,_i), idx(0)
        {
        if (existNodes())
            {
            orientate(); // enforce the correct orientation
 
            Eigen::Matrix3d J;
            double detJ = Jacobian(J);
            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 = [] ( const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
            }
        else
            { std::cerr<<"Error: Tet constructor has out of bound index\n"; }
        }
 
    Eigen::Matrix<double,N,Nodes::DIM> da;
 
    inline 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;
        }
 
    inline void interpolation(const 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(const 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(const std::function<double(Nodes::Node, Nodes::index)>& getter,
            const 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(const Eigen::Ref<const 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,
                                              const Eigen::Matrix<double,N,1> &x) const;
 
    Eigen::Matrix<double,N,1> calcOffDiagBlock(const Nodes::index idx) const;
 
 
    void DMI_bulk_OffDiag(timing const &prm_t, const double DMIbulk_bis,
                          Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE) const;
 
    void DMI_bulk_BE(const double cst,
                     Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> U,
                     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,Nodes::DIM,N>> BE) const;
 
    void add_drift_BE(const Eigen::Ref<const Eigen::Matrix<double,NPI,1>> 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(const 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(const Eigen::Ref<const Eigen::Matrix3d> e ,
            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;
 
    Eigen::Matrix<double,Nodes::DIM,NPI> Hst_order2_contribution(const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,NPI>> Hst,
            const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,NPI>> U ,
            const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,NPI>> V) const;
 
    void integrales( const Tetra::prm &param, const timing &prm_t,
                    const std::function< Eigen::Matrix<double,Nodes::DIM,NPI> (void)>& calc_Hext,
                    Nodes::index idx_dir, double Vdrift, const AbsorbingEnds &absorb);
 
    double exchangeEnergy(const Tetra::prm &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 DMIBulkEnergy(const Tetra::prm &param,
                         Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u,
                         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(const Tetra::prm &param,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
 
    double cubicAnisotropyEnergy(const Tetra::prm &param,
            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
 
    Eigen::Matrix<double,Tetra::NPI,1> charges(const double &Ms,
            const std::function<Eigen::Vector3d(const Nodes::Node&)> &getter) const override;
 
    double demagEnergy(const Tetra::prm &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(const Tetra::prm &param  ,
            const Eigen::Ref<const Eigen::Vector3d> Hext ,
            const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> u ) const;
 
    double zeemanEnergy(const Tetra::prm &param ,
                        double fieldAmp ,
                        const std::vector<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> &spaceField,
                        const Eigen::Ref<const Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> u) const;
 
    double Jacobian(Eigen::Ref<Eigen::Matrix3d> J) const;
 
    double calc_vol(void) const;
 
    int idx;
 
    std::function<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> H)> extraField;
 
    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;
        }
 
    Eigen::Matrix<double,Nodes::DIM,NPI> gradV(const std::vector<double> &V) const;
 
private:
    void orientate(void) override;
    };  // end class Tetra
 
    Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> calc_Hst(const Tetra::Tet &tet,
                                                         const double prefactor,
                                                         const std::vector<Eigen::Vector3d> &s);
}   // end namespace Tetra
 
#endif /* tetra_h */