SM_MacroCellsDemagnetizedField.h

#ifndef SM_MacroCellsDemagnetizedField_H
#define SM_MacroCellsDemagnetizedField_H
 
//base classes
#include "SM_Object.h"
 
//macro cells network header
#include "SM_MacroCellsNetwork.h"
 
//magnetization field defined on macro cells network header
#include "SM_MacroCellsMagnetizationField.h"
 
//storage matrix header
#include "SM_PackedBlockSymmetricMatrix.h"
 
 
class SM_MacroCellsDemagnetizedField : public virtual SM_Object  {
    
 
public:
    static constexpr tFlag PACKED_STORAGE=1;
 
    static constexpr tFlag NO_STORAGE=0;
private :
 
    //type class
    typedef SM_Object SuperClass;
    typedef SM_MacroCellsDemagnetizedField SelfClass;
    
    //single attributes
    //-----------------
  
    //constants for computing H
    tReal mLambdaHs;  
    tReal mLambdaH;
 
    
    //type of satorage of the demagnetized matrix in {NO_STORAGE, PACKED_STORAGE}
    tFlag mStorage;
    
   
   
    //demagnetized field Hdem  at each macro cell
    SM_RealField mHdem;
 
    //demagnetized matrix
    CORE_UniquePointer<SM_PackedBlockMatrix<tReal,SM_Constants::DIM>> mMatrix;
 
 
    //multiple associations
    //---------------------
 
    
   
protected:
    // CONSTRUCTORS 
    SM_MacroCellsDemagnetizedField(void) {
        mLambdaHs=0;
        mLambdaH=0;
       
        mStorage=NO_STORAGE;
    }
    
    // DESTRUCTORS 
    virtual ~SM_MacroCellsDemagnetizedField(void) {
    }
   
    
public : 
 
    //Instance building
    //=================
 
   
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+this->getContentsMemorySize();
    }
 
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mHdem.getContentsMemorySize();
        mem+=(mMatrix.get()==null)?0:mMatrix->getContentsMemorySize();
        return mem;
    }
 
    virtual CORE_UniquePointer<SM_MacroCellsDemagnetizedField> newInstance() const=0;
    
public:
    inline void setStorage(const tFlag& f) {
        mStorage=f;
    }
    inline const tFlag& getStorage() const {
        return mStorage;
    }
    
   
    inline const SM_RealField& getField() const {
        return mHdem;
    }
 
    inline  SM_RealField& getField()  {
        return mHdem;
    }
 
 
 
public:
  
    inline tBoolean hasMatrix() const {
        return (mMatrix.get()!=null);
    }
    inline void  setMatrix(CORE_UniquePointer<SM_PackedBlockMatrix<tReal,SM_Constants::DIM>> matrix) {
        mMatrix=std::move(matrix);
    }
    inline SM_PackedBlockMatrix<tReal,SM_Constants::DIM>& getMatrix() {
        return *mMatrix.get();
    }
    inline const SM_PackedBlockMatrix<tReal,SM_Constants::DIM>& getMatrix() const {
        return *mMatrix.get();
    }
 
    
 
public:
    inline const tReal& getLambdaH() const {
        return mLambdaH;
    }
    inline const tReal& getLambdaHs() const {
        return mLambdaHs;
    }
 
     //discretize the field
    //======================
public:
    virtual  void discretize(const SM_Material& material,const SM_MacroCellsNetwork& network);
    
 
    //demagnetized matrix build by macro cells
    //=========================================
public:
    virtual void computeMatrix(const SM_Material& material,
                               const SM_MacroCellsNetwork& network)=0;
    
protected:
    inline static void BuildPackedBlockMatrix(tReal lambdaH,const tReal& lambdaHs,const tReal& Wi,
                                              tInteger i,tInteger j,
                                              const tReal *Xi,const tReal* Xj,
                                              tReal *vBij,
                                              tReal& rij,std::array<tReal,SM_Constants::DIM>& U) {
 
        //zero tolerance
        tReal eps=1.e-9;
        
        //\f$ r_{ij}=|X_j-X_i| f$
        //\f$ U_{ij}=\frac{1}{r_ij}. (X_j-X)\f$
        //\f$ H_i=\lambda_H \sum_{j=0, j \neq i }^{j=P-1} (3*<U_ij,S_j>*U_ij}-S_j))/r_{ij}^3  - \frac{\lambda_{H_s}}{N_i} S_i  \f$ 
        //\f$ H_i= \sum_{ij} M_{ij}. S_j \f$
        // M_{ii}^{pq}= - \frac{\lambda_{H_s}}{N_i}  . \delta_{pq}
        // M_{ij]^{pq}=\frac{\lambdaH}{r_{ij}^3}  (3.Up.U.q -\delta_{pq})
        
        
        //I_{ij]=j.(j+1)/2+i
        //next row : i+1
        //next column C_{j+1}-C_{j]=(j+1).(j+2)/2 - j.(j+1)/2=0.5.(j+1).(j+2-j)=j+1
        //next diagonal D_{j+1}-D_j=0.5.(j+1).(j+2)+j+1-0.5.j.(j+1)-j=0.5.(j+1).(j+2-j)+1=j+2
        
        //compute U,rij
 
        //rij^2
        rij=0;
        for(auto& Uk:U) {
            //Uij
            Uk=(*Xj)-(*Xi);
            //rij
            rij+=Uk*Uk;
            //ietartors at next coorinates
            Xj++;
            Xi++;
        }
 
        if (rij<eps) {
            BuildPackedDiagonalBlockMatrix(Wi,lambdaHs,vBij);
        } else {
            
            //lambdaH=LambdaH/rij^2
            lambdaH/=rij;
            //rij
            rij=sqrt(rij);
            //LambdaH=LambdaH/rij^3
            lambdaH/=rij;
            
            //U is normalized
            for(auto& Ui:U) Ui/=rij;
            
            const tReal *Ui=null,*Uj=U.data();
            for(j=0;j<SM_Constants::DIM;j++) {//loop on column j in [0,dim[
                //H_ij=3.Ui.Uj.lambdaH - lambdaH.delta_ij
                Ui=U.data();
                for(i=0;i<j;i++) {//loop on row i in[0,j[ of column j  
                    (*vBij)=3;
                    (*vBij)*=(*Ui);
                    (*vBij)*=(*Uj);
                    (*vBij)*=lambdaH;
                    //iterator at next row i
                vBij++;
                Ui++;
                }
                //i=j diagonal element
                (*vBij)=3;
                (*vBij)*=(*Ui);
                (*vBij)*=(*Uj);
                (*vBij)-=1.0L;
                (*vBij)*=lambdaH;
                
                
                //iterator at next colum j
                vBij++;
                Uj++;
            }//end packed  values
        }//end null rij condition
    }
        
    
    inline static void BuildPackedSupBlockMatrix(tReal lambdaH,
                                                 tInteger i,tInteger j,
                                                 const tReal *Xi,const tReal* Xj,
                                                 tReal *vBij,
                                                 tReal& rij,std::array<tReal,SM_Constants::DIM>& U) {
        //\f$ r_{ij}=|X_j-X_i| f$
        //\f$ U_{ij}=\frac{1}{r_ij}. (X_j-X)\f$
        //\f$ H_i=\lambda_H \sum_{j=0, j \neq i }^{j=P-1} (3*<U_ij,S_j>*U_ij}-S_j))/r_{ij}^3  - \frac{\lambda_{H_s}}{N_i} S_i  \f$ 
        //\f$ H_i= \sum_{ij} M_{ij}. S_j \f$
        // M_{ii}^{pq}= - \frac{\lambda_{H_s}}{N_i}  . \delta_{pq}
        // M_{ij]^{pq}=\frac{\lambdaH}{r_{ij}^3}  (3.Up.U.q -\delta_{pq})
        
        
        //I_{ij]=j.(j+1)/2+i
        //next row : i+1
        //next column C_{j+1}-C_{j]=(j+1).(j+2)/2 - j.(j+1)/2=0.5.(j+1).(j+2-j)=j+1
        //next diagonal D_{j+1}-D_j=0.5.(j+1).(j+2)+j+1-0.5.j.(j+1)-j=0.5.(j+1).(j+2-j)+1=j+2
        
        //compute U,rij
 
        //rij^2
        rij=0;
        for(auto& Uk:U) {
            //Uij
            Uk=(*Xj)-(*Xi);
            //rij
            rij+=Uk*Uk;
            //iterators at next coorinates
            Xj++;
            Xi++;
        }
        
        //lambdaH=LambdaH/rij^2
        lambdaH/=rij;
        //rij
        rij=sqrt(rij);
        //LambdaH=LambdaH/rij^3
        lambdaH/=rij;
    
        //U is normalized
        for(auto& Ui:U) Ui/=rij;
        
        const tReal *Ui=null,*Uj=U.data();
        for(j=0;j<SM_Constants::DIM;j++) {//loop on column j in [0,dim[
            //H_ij=3.Ui.Uj.lambdaH - lambdaH.delta_ij
            Ui=U.data();
            for(i=0;i<j;i++) {//loop on row i in[0,j[ of column j  
                (*vBij)=3;
                (*vBij)*=(*Ui);
                (*vBij)*=(*Uj);
                (*vBij)*=lambdaH;
                
#ifdef DEBUG
                if (std::isnan(*vBij)) {
                    throw CORE_Exception("stochmagnet/operators/fields",
                                         "SM_MacroCellsDemagnetizedField::BuildPackedSupBlockMatrix(Xi="+functions_array::toString(3,Xi-3)+
                                         " Xj="+functions_array::toString(3,Xj-3)+",LH="+std::to_string(lambdaH)+" rij="+std::to_string(rij)+" U:"+functions_array::toString(U),
                                         "nan values");
                }
#endif
        
                //iterator at next row i
                vBij++;
                Ui++;
            }
            //i=j diagonal element
            (*vBij)=3;
            (*vBij)*=(*Ui);
            (*vBij)*=(*Uj);
            (*vBij)-=1.0L;
            (*vBij)*=lambdaH;
            
            
            //iterator at next colum j
            vBij++;
            Uj++;
        }
    }
    inline static void BuildDiagonalBlockMatrix(const tReal& V,
                                                const tReal& lambdaHs,
                                                tReal& diag) {
 
        diag=lambdaHs;
        diag*=-1;
        diag/=V;
#ifdef DEBUG
        if (std::isnan(diag)) {
            throw CORE_Exception("stochmagnet/operators/fields",
                                 "SM_MacroCellsDemagnetizedField::BuildDiagonalBlockMatrix("+std::to_string(V)+","+std::to_string(lambdaHs)+")",
                                 "nan values");
        }
#endif
        
    }
    inline static void BuildPackedDiagonalBlockMatrix(const tReal& V,
                                                      const tReal& lambdaHs,
                                                      tReal *iDiag) {
 
        memset(iDiag,0,sizeof(tReal)*SM_PackedBlockMatrix<tReal,SM_Constants::DIM>::BLOCK_SIZE);
        BuildDiagonalBlockMatrix(V,lambdaHs,(*iDiag));
        //std::cout<<"diag:"<<(*iDiag)<<"\n";
        const tReal &d=(*iDiag);
        for(tDimension j=2;j<=SM_Constants::DIM;j++) {//loop on column j in [0,dim[
            iDiag+=j;
            (*iDiag)=d;
        }
         
        
    }
  
 
 
protected:
   
    void computeBlocksSymmetricMatrix(const tReal& lambdaH,const tReal& lambdaHs,
                                      const tReal* vMCsVolume,
                                      const tReal *vMCsX,
                                      const tInteger& start,const tInteger& end,
                                      tReal  *supBlocksValues,
                                      tReal *diagBlocksValues) const;
 
    void computeBlocksGeneralMatrix(const tReal& lambdaH,const tReal& lambdaHs,
                                    const tReal* vMCsVolumes,
                                    const tInteger& nX,const tReal *vX,
                                    const tInteger& nP,const tReal *vP,
                                    tReal *supBlocksValues) const;
 
 
public:
    virtual void computeField(const SM_MacroCellsNetwork& network,
                              const SM_MacroCellsMagnetizationField& M)=0;
    
 
protected:
    void computeFieldSlice(const tReal& lambdaH,const tReal& lambdaHs,
                           const tReal *P,const tReal *eP,
                           const tReal* M_P,
                           const tInteger& nMacroCells,
                           const tReal *X,
                           const tReal *M_X,
                           const tReal* iMCsVolume,
                           tReal *H) const;
private:
    void computeFieldSliceInSeparatedNetworkStorage(const tReal& lambdaH,const tReal& lambdaHs,
                                                    const tReal *P,const tReal *eP,
                                                    const tReal* M_P,
                                                    const tInteger& nMacroCells,
                                                    const tReal *X,
                                                    const tReal *M_X,
                                                    const tReal* iMCsVolume,
                                                    tReal *H) const;
 
public:
    void computeSeparatedFieldSlice(const tReal& lambdaH,const tReal& lambdaHs,
                                    const tReal *P,const tReal *eP,
                                    const tReal* M_P,
                                    const tInteger& nMacroCells,
                                    const tReal *X,
                                    tReal *H) const;
        
 
    // String representation
    // =======================
public:
 
      
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString()<<"\n";
        ret<<"\t storage of demagnetized field: "<<((int) mStorage)<<"\n";
        ret<<"\t Hdem:\n";
        tIndex s=mHdem.getElementsNumber()/10;
        s=(s==0)?1:s;
        for(tIndex i=0;i<mHdem.getElementsNumber();i+=s) {
            ret<<"\t\t ["<<i<<"]="<<functions_array::toString(SM_Constants::DIM,mHdem(i))<<"\n";
        }
        return ret.str();
    }
 
public:
    
   
  
};
 
 
#endif