tetra.h

Go to the documentation of this file.

#ifndef tetra_h
#define tetra_h
 
#include <set>
 
#include "spinTransferTorque.h"
#include "facette.h"
#include "node.h"
#include "time_integration.h"
#include "element.h"
 
namespace Tetra
    {
const int N = 4;   
const int NPI = 5; 
const double epsilon =
        EPSILON; 
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]}};
 
static const Eigen::Matrix<double,N,NPI> eigen_a = [] {
        Eigen::Matrix<double,N,NPI> tmp; tmp << 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];
                                            return tmp; }();
 
struct prm
    {
    std::string regName; 
    double alpha_LLG;    
    double A;            
    double J;            
    double K;            
    Eigen::Vector3d uk;         
    double K3;   
    Eigen::Vector3d ex; 
    Eigen::Vector3d ey; 
    Eigen::Vector3d ez; 
    STT p_STT;   
    inline void infos()
        {
        std::cout << "volume region name = " << regName << std::endl;
        std::cout << "alpha_LLG = " << alpha_LLG << std::endl;
        std::cout << "A = " << A << std::endl;
        std::cout << "J = " << J << std::endl;
 
        if (K != 0)
            {
            std::cout << "K*uk =" << K << "*[ " << uk << "]" << std::endl;
            }
 
        if (K3 != 0)
            {
            std::cout << "K3 = " << K3 << "; ex=[ " << ex << "], ey=[" << ey << "], ez=[" << ez
                      << "]\n";
            }
        };
    };
 
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;
            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]; }
            }
        // do nothing lambda's (usefull for spin transfer torque)
        extraField = [] ( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> ) {};
        extraCoeffs_BE = [](double, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
                            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
                            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
                            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>>,
                            Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,N>>) {};
        }
 
    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;
        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;
 
    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<void( Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> Hext )> 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 anisotropyEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,NPI>> u) const;
 
    void charges(Tetra::prm const &param,
                 std::function<Eigen::Vector3d(Nodes::Node)> getter,
                 std::vector<double> &srcDen,
                 int &nsrc) 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,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;
 
    std::function<void(double Js, 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)> extraCoeffs_BE;
 
    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);
    }   // end namespace Tetra
 
#endif /* tetra_h */