EMM_CubicAnisotropyOperator.h

Go to the documentation of this file.

#ifndef EMM_CubicAnisotropyOperator_H
#define EMM_CubicAnisotropyOperator_H


#include "EMM_AnisotropyOperator.h"


DEFINE_SPTR(EMM_CubicAnisotropyOperator);
class EMM_CubicAnisotropyOperator : public virtual EMM_AnisotropyOperator {
    
    SP_OBJECT(EMM_CubicAnisotropyOperator);
    // ATTRIBUTES
    
    
private:

    
public:
   

    // ASSOCIATION
 
private:
    
   
protected: 

    // METHODS
  
    // CONSTRUCTORS 
  
    EMM_CubicAnisotropyOperator(void);
  
    // DESTRUCTORS 
    
    virtual ~EMM_CubicAnisotropyOperator(void);   
  
  
public:
   
   
    // SET methods
    // ===========
   
    
    // GET methods
    // ===========
    virtual tULLInt getMemorySize() const {
        return EMM_AnisotropyOperator::getMemorySize();
    }

    //DATA Methods
    //===========
    
public:
    /* \brief build the symmetric matrix of size 6
     * @param nU: memory size of directions U
     * @param U: anisotropy direction of size 3
     * @param A: anisotropy matrix of size 6
     *
     * copy the nU elements for U into the morse matrix A
     *
     * @return the size of matrix A (nU)
     */
    static inline tUSInt BuildDiscretizedMatrix(const tUCInt& nU,
                                                const tReal* U,
                                                tReal* A) {
        //copy only the directions
        const tReal *Uk=U;
        tReal *Ak=A;
        for (tDimension k=0;k<nU;k++) {
            (*Ak)=(*Uk);
            Uk++;
            Ak++;
        }
        return nU;
    }
 
    //MAGNETIC Field Methods
    //======================
    
protected:
 
    virtual void computeMagneticExcitationField(const tUIndex& nCells,
                                                const tDimension& dim,
                                                const EMM_RealArray& sigma,
                                                const tReal* M,
                                                tReal *H) const;
    
    
    virtual tBoolean computeMagneticExcitationFieldGradient(const tUIndex& nCells,
                                                            const tDimension& dim,
                                                            const EMM_RealArray& sigma,
                                                            const tReal *M,
                                                            const tReal *D,
                                                            tReal  *gradH)const;

    
public:
    static inline tUSInt ComputeMagneticExcitation(const tReal& S,
                                                   const tReal& K1,
                                                   const tReal& K2,
                                                   const tUCInt& nA,
                                                   const tReal* A,
                                                   const tDimension& dim,
                                                   const tReal* M,
                                                   tReal* H,
                                                   const tUCInt& ldW,
                                                   tReal *work) {


       
     
        
       
        //get the number of directions
        tUCInt d,nDirs=nA/dim;
        tDimension k;

        
        //compute Mud=Mu[d]=\sigma_i*<M,Ud> = Mu[d+3]=Mut
        memset(work,0,ldW*sizeof(tReal));
        
        const tReal* Ud=A;
        tReal *Mud=work;
        tReal *Mut=work+dim;
        const tReal *Mk;
        for(d=0;d<nDirs;d++) {
            Mk=M;
            for (k=0;k<dim;k++) {
                (*Mud)+=(*Ud)*(*Mk);
                Ud++;
                Mk++;
            }
            (*Mud)*=S;
            (*Mut)=(*Mud);
            Mud++;
            Mut++;
        }
        
        //dH=Mu[d] * [ K1* (Mu[d+1]*Mu[d+1]+Mu[d+2]*Mu[d+2]) + K2* (Mu[d+1]*Mu[d+1]*Mu[d+2]*Mu[d+2]) ]
        //H-=dH  U[d]
        
        const tReal *Mu0=work;
        const tReal *Mu1=(Mu0+1);
        const tReal *Mu2=(Mu1+1);
        tReal *Hk=H;
        tReal dH;
        const tReal *Adk=A;
        for (d=0;d<nDirs;d++) {
            //compute dH
            dH=(*Mu0)*( K1*((*Mu1)*(*Mu1)+(*Mu2)*(*Mu2))+
                        K2*((*Mu1)*(*Mu1)*(*Mu2)*(*Mu2)) );
            Mu0++;
            Mu1++;
            Mu2++;
            
            //add Ud.dH to H
            Hk=H;
            for (k=0;k<dim;k++) {
                (*Hk)-=dH*(*Adk);
                Adk++;
                Hk++;
            }
        }
        return nDirs*dim;
    };
    


public:
    static inline tUSInt ComputeMagneticExcitationGradient(const tReal& S,
                                                           const tReal& K1,
                                                           const tReal& K2,
                                                           const tUCInt& nA,
                                                           const tReal* A,
                                                           const tDimension& dim,
                                                           const tReal* M,
                                                           const tReal* D,
                                                           tReal* gradH,
                                                           const tUCInt& ldW,
                                                           tReal *work) {
       
        
        //get the number of directions
        tUCInt d,nDirs=nA/dim;
        tDimension k;
        
        
        //compute Mud=Mu[d]=\sigma_i*<M,Ud> = Mu[d+3]
        //compute Dud=Du[d]=\sigma_i*<D,Ud> = Mu[d+6]=Mu[d+9]
        const tReal *Ud=A;
        tReal *Mud=work;
        tReal *Mut=work+3;

        const tReal *Mk;
        const tReal *Dk;
        tReal *Dud=work+6;
        tReal *Dut=work+9;
        tReal *gradHk;
        tReal dG;
        
        //initialize the work array to 0
        memset(work,0,ldW*sizeof(tReal));
        
        for(d=0;d<nDirs;d++) {
            Mk=M;
            Dk=D;
            for (k=0;k<dim;k++) {
                (*Mud)+=(*Ud)*(*Mk);
                (*Dud)+=(*Ud)*(*Dk);
                Ud++;
                Mk++;
                Dk++;
            }
            (*Mud)*=S;
            (*Dud)*=S;
            (*Mut)=(*Mud);
            (*Dut)=(*Dud);
            Mud++;
            Mut++;
            Dud++;
            Dut++;
            
        }
        
        //dG=Mu[d] * [ 2K1* (Mu[d+1]*Du[d+1]+Mu[d+2]*Du[d+2]) 
        //            +2K2* (Mu[d+1]*Du[d+1]*Mu[d+2]*Mu[d+2] + Mu[d+1]*Mu[d+1]*Mu[d+2]*Du[d+2]) ] +
        //   Du[d] * [ K1* (Mu[d+1]*Mu[d+1]+Mu[d+2]*Mu[d+2]) +
        //             K2* (Mu[d+1]*Mu[d+1]*Mu[d+2]*Mu[d+2]) ]
        //=>
        //dG=2.Mu[d] * [ K1* (Mu[d+1]*Du[d+1]+Mu[d+2]*Du[d+2]) 
        //              +K2* Mu[d+1]*Mu[d+2]*(Du[d+1]*Mu[d+2] + Mu[d+1]*Du[d+2]) ] +
        //     Du[d] * [ K1* (Mu[d+1]*Mu[d+1]+Mu[d+2]*Mu[d+2]) +
        //               K2* (Mu[d+1]*Mu[d+1]*Mu[d+2]*Mu[d+2]) ]
        
        //gradH-=dG  U[d]
        
        const tReal *Mu0=work;
        const tReal *Mu1=(Mu0+1);
        const tReal *Mu2=(Mu1+1);
        const tReal *Du0=work+6;
        const tReal *Du1=(Du0+1);
        const tReal *Du2=(Du1+1);
        const tReal *Adk=A;
        
        for (d=0;d<nDirs;d++) {
            //compute dG
            dG=2*(*Mu0)*( K1*((*Mu1)*(*Du1)+(*Mu2)*(*Du2)) +
                          K2*(*Mu1)*(*Mu2)*((*Du1)*(*Mu2)+(*Mu1)*(*Du2)))+
                (*Du0)*( K1*((*Mu1)*(*Mu1)+(*Mu2)*(*Mu2))+
                         K2*((*Mu1)*(*Mu1)*(*Mu2)*(*Mu2)) );
            
            Mu0++;
            Mu1++;
            Mu2++;
            Du0++;
            Du1++;
            Du2++;
            
            //add Ud.dH to H
            gradHk=gradH;
            for (k=0;k<dim;k++) {
                (*gradHk)-=dG*(*Adk);
                Adk++;
                gradHk++;
            }
            
        }
        return dim*nDirs;
        
    };




    //ENERGY methods
    //==============
    
    tReal computeEnergy(const tUIndex& nCells,const tDimension& dim,
                        const EMM_RealArray& sigma,const tReal* M) const;
    

  
public:
    static inline void  ComputeEnergy(const tReal& S,
                                      const tReal& K1,
                                      const tReal& K2,
                                      const tUCInt& nA,
                                      const tReal* A,
                                      const tDimension& dim,
                                      const tReal* M,
                                      const tUCInt& ldW,
                                      tReal* work,
                                      tReal& E) {
        tReal dE;
        E=0;
        //compute Mud=Mu[d]=\sigma_i*<M,Ud> = Mu[d+3]=Mut
        memset(work,0,ldW*sizeof(tReal));
        tReal *Muk,*Mud=work;
        const tReal *Mk;
        tUCInt d,nDirs=nA/dim;//number of anisotropy directions in [0,3[
        tDimension k;
        
        const tReal *Adk=A;
        for(d=0;d<nDirs;d++) {
            Mk=M;
            for (k=0;k<dim;k++) {
                (*Mud)+=(*Adk)*(*Mk);
                Adk++;
                Mk++;
            }
            (*Mud)*=S;
            Mud++;
        }

        //K2 contribution
        if (K2!=0) {
            dE=K2;
            Mud=work;
            for (d=0;d<nDirs;d++) {
                dE*=(*Mud)*(*Mud);
                Mud++;
            }
            E+=dE;
        }
        
        //K1 contribution
        if (K1!=0) {
            Mud=work;
            for (d=0;d<nDirs;d++) {
                dE=0;
                Muk=work;
                for (k=0;k<d;k++) {
                    dE+=(*Muk)*(*Muk);
                    Muk++;
                }
                E+=K1*(*Mud)*(*Mud)*dE;
                Mud++;
            }
            
        }

   
        
    }




    
        
};

#endif