SM_LandauLifschitzFunction.h

#ifndef SM_LandauLifschitzFunction_H
#define SM_LandauLifschitzFunction_H
 
//inherit class header
#include "SM_Object.h"
 
//constants value
#include "SM_Constants.h"
//field heaer
#include "SM_Field.h"
 
class SM_LandauLifschitzFunction : public  SM_Object {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
    
    using LambdaF = std::function<void(const tReal&,const tReal& ,tReal&)>;
    LambdaF mAlphaF;
    
    //gamma function
    tReal mGamma;
 
 
    //lambda parameter
    tReal mLambda;
  
    //alpha  parameter  \f$ \alpha=\alpha(\gamma,\lambda)=-\frac{\gamma]{1+\lambda^2} \f$
    tReal mAlpha;
 
    
    //normalized mu
    tBoolean mIsSNormalized;
  
protected:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_LandauLifschitzFunction(void) {
    
        mGamma=SM_Constants::GAMMA;
        mLambda=0.1;
        
        mIsSNormalized=false;
        mAlphaF=([](const tReal& gamma,const tReal& lambda,tReal& alpha) {
                     //alpha=lambda^2
                     alpha=lambda;
                     alpha*=alpha;
                     alpha+=1;
                     //alpha=gamma/(1+lambda^2)
                     alpha=gamma/alpha;
                     //alpha=-gamma/(1+lambda^2)
                     alpha*=-1;
                 });
        //compute alpha
        mAlphaF(mGamma,mLambda,mAlpha);
                        
    }
  
  
    // DESTRUCTORS 
  
    
    virtual ~SM_LandauLifschitzFunction(void) {
    }
  
  
public:
 
 
      //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
 
    
    
    // SET methods 
 
    inline void copy(const SM_LandauLifschitzFunction& LL) {
        mGamma=LL.getGamma();
        mLambda=LL.getLambda();
        mAlphaF(mGamma,mLambda,mAlpha);
        mIsSNormalized=LL.isSNormalized();
    }
 
   
 
    inline void setAlphaFunction(const LambdaF& L) {
        mAlphaF=L;
    }
    inline void setGamma(const tReal& gamma) {
        mGamma=gamma;
        mAlphaF(mGamma,mLambda,mAlpha);
    }
    inline void setLambda(const tReal& v) {
        mLambda=v;
        mAlphaF(mGamma,mLambda,mAlpha);
    }
 
    inline void setAlpha(const tReal& v) {
        mAlpha=v;
    }
public:
   
    
    inline void setIsSNormalized(const tBoolean& v) {
        mIsSNormalized=v;
    }
    
    // GET methods
    // ============
    
    inline const tReal&  getGamma() const {
        return mGamma;
    }
    
    inline const tReal&  getAlpha() const {
        return mAlpha;
    }
    
    inline const tReal&  getLambda() const {
        return mLambda;
    }
  
    inline tBoolean  isSNormalized() const {
        return mIsSNormalized;
    }
    
  
public:
    // OTHERS methods
    inline void computeFunction(const tIndex& nParticles,const SM_RealField& S,const SM_RealField& B,SM_RealField& F) const {
 
        tIndex nB=B.getElementsNumber();
        if (nB>0) {
            //set the output field dimension F=LL(S,B) size of F = size of B
            if (F.getElementsNumber()!=B.getElementsNumber()) F.setElementsNumber(nB);
            
            computeFunction(nParticles,&S[0],B.getElementsNumber(),&B[0],&F[0]);
        }
    }
    inline void computeFunction(const tIndex& nSParticles,const tReal *S,
                                const tIndex& nBParticles,const tReal *B,
                                tReal *F) const {
        if (mIsSNormalized) computeNLLFunction(nSParticles,S,nBParticles,B,F);
        else computeLLFunction(nSParticles,S,nBParticles,B,F);
    }
 
protected:
 
    virtual void computeLLFunction(const tIndex& nSParticles,const tReal *S,
                                   const tIndex& nBParticles,const tReal *B,tReal *F) const=0;
    
    virtual void computeNLLFunction(const tIndex& nSParticles,const tReal *S,
                                    const tIndex& nBParticles,const tReal *B,tReal *F) const=0;
    
 
    
    inline void computeSliceLLFunction(tIndex start,tIndex end,
                                       const tIndex& nSParticles,const tReal *S,
                                       const tIndex& nBParticles,const tReal *B,tReal *F) const {
        
        //particle index -> memory index
        start*=SM_Constants::DIM;
        end*=SM_Constants::DIM;
 
        //sSB=<S,B>
        //sSS=<S,S>
        tReal sSB,sSS;
        //begin iterator on mu
        const tReal *iS=S;iS+=start;
        //end iterator on mu
        const tReal *eS=S;eS+=end;
        
        //begin iterator on B
        const tReal *iB=B;
        if (nBParticles==nSParticles) iB+=start;
        
        //iterator on F
        tReal *iF=F;
        if (nBParticles==nSParticles) iF+=start;
        
        tReal S0,S1,S2,B0,B1,B2;
        while (iS!=eS)  {//loop on point
            
            
            // A = alpha(  S \wedge B + lambda ( <\S,B> \S -  <\S,\S> B ) )
            S0=(*iS);iS++;
            S1=(*iS);iS++;
            S2=(*iS);iS++;
            
            B0=(*iB);iB++;
            B1=(*iB);iB++;
            B2=(*iB);iB++;
            
            sSB=S0*B0+S1*B1+S2*B2;
            sSS=S0*S0+S1*S1+S2*S2;
            (*iF)= mAlpha*((S1 * B2 - S2 * B1) + mLambda * (S0 * sSB - B0 * sSS));iF++;
            (*iF)= mAlpha*((S2 * B0 - S0 * B2) + mLambda * (S1 * sSB - B1 * sSS));iF++;
            (*iF)= mAlpha*((S0 * B1 - S1 * B0) + mLambda * (S2 * sSB - B2 * sSS));iF++;
            
#ifdef DEBUG
            if (std::isnan(iF[-1]) || std::isnan(iF[-2]) || std::isnan(iF[-3])) {
                throw CORE_Exception("stochmagnet/landauLifschitz",
                                     "SM_LandauLifschitzFunction::computeSliceLLFunction("+std::to_string(start)+","+std::to_string(end)+",..)",
                                     "nan detected : S=("+std::to_string(S0)+","+std::to_string(S1)+","+std::to_string(S2)+") B=("+std::to_string(B0)+","+std::to_string(B1)+","+std::to_string(B2)+") F=("+std::to_string(iF[-3])+","+std::to_string(iF[-2])+","+std::to_string(iF[-3])+")");
            }
#endif
            
        }//end loop on values
    }
    
    inline void computeSliceNLLFunction(tIndex start,tIndex end,
                                        const tIndex& nSParticles,const tReal *S,
                                        const tIndex& nBParticles,const tReal *B,tReal *F) const {
 
        //particle index -> memory index
        start*=SM_Constants::DIM;
        end*=SM_Constants::DIM;
       
        
        //sSB=<S,B>
        //sSS=<S,S>
        tReal sSB;
        //begin iterator on S
        const tReal *iS=S;iS+=start;
        //end iterator on S
        const tReal *eS=S;eS+=end;
        //begin iterator on B
        const tReal *iB=B;
        if (nBParticles==nSParticles) iB+=start;
        
        //iterator on F
        tReal *iF=F;
        if (nBParticles==nSParticles) iF+=start;
        
        tReal S0,S1,S2,B0,B1,B2;
        
        while (iS!=eS) {//loop on point
            
            // A = alpha(  S \wedge B + lambda ( <\S,B> \S -  <\S,\S> B ) )
            S0=(*iS);iS++;
            S1=(*iS);iS++;
            S2=(*iS);iS++;
            
            B0=(*iB);iB++;
            B1=(*iB);iB++;
            B2=(*iB);iB++;
            
            sSB=S0*B0+S1*B1+S2*B2;
            //sSS=S0*S0+S1*S1+S2*S2;
            //sSS=1
            (*iF)= mAlpha*((S1 * B2 - S2 * B1) + mLambda * (S0 * sSB - B0 ));iF++;
            (*iF)= mAlpha*((S2 * B0 - S0 * B2) + mLambda * (S1 * sSB - B1 ));iF++;
            (*iF)= mAlpha*((S0 * B1 - S1 * B0) + mLambda * (S2 * sSB - B2 ));iF++;
 
#ifdef DEBUG
            if (std::isnan(iF[-1]) || std::isnan(iF[-2]) || std::isnan(iF[-3])) {
                throw CORE_Exception("stochmagnet/landauLifschitz",
                                     "SM_LandauLifschitzFunction::computeSliceNLLFunction("+std::to_string(start)+","+std::to_string(end)+",..)",
                                     "nan detected : S=("+std::to_string(S0)+","+std::to_string(S1)+","+std::to_string(S2)+") B=("+std::to_string(B0)+","+std::to_string(B1)+","+std::to_string(B2)+") F=("+std::to_string(iF[-3])+","+std::to_string(iF[-2])+","+std::to_string(iF[-3])+")");
            }
#endif
            
        }//end loop on particle
   
    }
 
    
public:
    
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<"LL gamma:"<<mGamma<<"\n";
        ret<<"LL alpha:"<<mAlpha<<"\n";
        ret<<"LL lambda:"<<mLambda<<"\n";
        ret<<((mIsSNormalized)?
              "LL(S,H)=alpha.( S ^ H + lambda.(<S,H>S - H )\n":
              "LL(S,H)=alpha.( S ^ H + lambda.(<S,H>S - |S|^2 H )\n");
        return ret.str();
    }
 
 
};
 
#endif