SM_AnisotropyOperator.hpp

#ifndef SM_AnisotropyOperator_HPP
#define SM_AnisotropyOperator_HPP
 
template<class I>
void SM_AnisotropyOperator<I>::computeMagneticFieldSlice(const tIndex& timeIndex,
                                                         const SM_Network& network,
                                                         const SM_Material& material,
                                                         const tIndex& startIndex,
                                                         const tIndex& endIndex,
                                                         const tReal *S,
                                                         const tBoolean& alpha,
                                                         const tIndex& nH,
                                                         tReal *H) const {
 
 
  
    //number of particles of the slice 
    tIndex s=endIndex-startIndex;
    
    //iterator on H at start index
    tReal *iH=H;
    if (nH>s) iH+=startIndex*SM_Constants::DIM;
    
    
    if (alpha==0) {   
        //H=0
        memset(iH,0,sizeof(tReal)*s*SM_Constants::DIM);
    }
   
 
    
      
    //anisotropy variables
    tReal K=material.getAdimensionizedAnisotropyEnergyFactor();
    //dividide by \mu_s depending on adimenonizized variable or not
    //mu_s=\tilde mu_s . mu_B for dimensionized variable
    //mu_s=\tilde mu_s  for adimenisonized variable
    K/=material.getAtomicSpinMoment();
    K*=material.getAdimensionizedDerivativeEnergyFactor();
    
    //begin anisotropy direction
    const tReal*U=&material.getAnisotropyDirections()[0];
 
    //end anisotropy direction
    const tReal *eU=U;
    eU+=material.getAnisotropyDirections().size();
        
        
    
    //phi value
    tReal phi;
    
    //gradPhi value
    std::array<tReal,SM_Constants::DIM> gradPhi;
    const tReal *iGradPhi=gradPhi.data();
    const tReal *eGradPhi=iGradPhi;eGradPhi+=SM_Constants::DIM;
        
            
 
    //begin iterator on S
    const tReal *iS=S;iS+=startIndex*SM_Constants::DIM;
    
    //end iterator on S
    const tReal *eS=S;eS+=endIndex*SM_Constants::DIM;
    
    //end iterator on coodinate of k
    const tReal *eSk=iS;eSk+=SM_Constants::DIM;
    
    //work values
    tReal W,W3;
    while (iS!=eS) {
        
        //compute phi & GradPhi
        computePhiDerivatives(U,eU,iS,eSk,phi,gradPhi,W,W3);
        
        //H+=K.GradPhi
        iGradPhi=gradPhi.data();
        while (iGradPhi!=eGradPhi) {
            (*iH)-=K*(*iGradPhi);
            iGradPhi++;
            iH++;
        }
        
        //iterator at next particle
        iS=eSk;
        eSk+=SM_Constants::DIM;
        
    }//end loop on particles
            
}
 
template<class I>
tReal SM_AnisotropyOperator<I>::computeEnergySlice(const tIndex& timeIndex,
                                                   const SM_Network& network,
                                                   const SM_Material& material,
                                                   const tIndex& startIndex,
                                                   const tIndex& endIndex,
                                                   const tReal *S) const   {
    
        
    //energy of the operator
    tReal E=0;
 
    //anisotropy constant
    const tReal &K=material.getAdimensionizedAnisotropyEnergyFactor();
  
    //begin anisotropy direction
    const tReal *U=&material.getAnisotropyDirections()[0];
    
    //end anisotropy direction
    const tReal *eU=U;
    eU+=material.getAnisotropyDirections().size();
        
    
    
    //begin iterator on S
    const tReal *iS=S;iS+=startIndex*SM_Constants::DIM;
    
    //end iterator on S
    const tReal *eS=S;eS+=endIndex*SM_Constants::DIM;
    
    //end iterator on coodinate of k
    const tReal *eSk=iS;eSk+=SM_Constants::DIM;
    
    
    tReal W,Ei;
    
    
    
    while (iS!=eS) {
        
        //compute the spin energy
        computePhi(U,eU,iS,eSk,Ei,W);
        
        //add to total energy
        E+=Ei;
        
        
        //iterator at next particles
        iS=eSk;
        eSk+=SM_Constants::DIM;
        
    }//end particles loop
    
    //  E= \sum_{i=0}^{i=P-1}  K.  \Phi
    E*=K;
    
    return E;
}//end method
 
 
#endif