SM_DemagnetizedOperator.h

Go to the documentation of this file.

#ifndef SM_DemagnetizedOperator_H
#define SM_DemagnetizedOperator_H
 
//base classes
#include "SM_TemplatedOperator.h"
 
//network class
#include "SM_Network.h"
 
 
class SM_DemagnetizedOperator : public SM_TemplatedOperator<SM_DemagnetizedOperator> {
    
    //attributes
private :
    
    
    //association
    CORE_BooleanArray mAreMagnetizedCells;
   
protected:
    // CONSTRUCTORS 
    SM_DemagnetizedOperator(void) {
        setName("Demagnetized");
    }
    
    // DESTRUCTORS 
    virtual ~SM_DemagnetizedOperator(void) {
    }
    
public : 
 
    //Instance building
    //=================
 
   
    inline static CORE_UniquePointer<SM_DemagnetizedOperator> New() {
        return CORE_UniquePointer<SM_DemagnetizedOperator>(new SM_DemagnetizedOperator(),
                                                           SM_DemagnetizedOperator::Delete());
    }
     
    virtual CORE_UniquePointer<SM_Operator> NewInstance() const override {
        return New();
    }
 
   
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
 
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SM_TemplatedOperator<SM_DemagnetizedOperator>::getContentsMemorySize();
        mem+=mAreMagnetizedCells.getSize()*sizeof(tBoolean);
        return mem;
    }
public:
 
    virtual void discretize(const SM_System& system);
 
    inline void updateMagnetizedCells(const SM_RealField& Mu) {
      
        tDimension k,dim=Mu.getDimension();
        mAreMagnetizedCells.setSize(Mu.getElementsNumber());
        auto Mu_jk=Mu.cbegin();
        tReal n2Muj;
        std::for_each(mAreMagnetizedCells.begin(),mAreMagnetizedCells.end(),
                      [&](auto& isM) {
                          n2Muj=0;
                          for (k=0;k<dim;k++) {
                              n2Muj+=(*Mu_jk)*(*Mu_jk);
                              Mu_jk++;
                          }
                          isM=(n2Muj>1.e-12);
                          
                      });
    }
    
    
    virtual void computeMagneticField(const tIndex& t,
                                      const SM_Network& network,
                                      const tReal *mu,
                                      tReal *B) const override {
        return computeTemplatedMagneticField(t,network,mu,B);
    }
    virtual void computeMagneticFieldAndEnergy(const tIndex& t,
                                               const SM_Network& network,
                                               const tReal *mu,
                                               tReal *B,
                                               tReal& E) const override{
        return computeTemplatedMagneticFieldAndEnergy(t,network,mu,B,E);
    }
    inline void computeTemplatedMagneticField(const tIndex& t,
                                              const SM_Network& network,
                                              const tReal *mu,
                                              tReal *B) const {
        tReal E=0;
        computeTemplatedMagneticFieldAndEnergy(t,network,mu,B,E);
    }
    inline void computeTemplatedMagneticFieldAndEnergy(const tIndex& t,
                                                       const SM_Network& network,
                                                       const tReal *mu,
                                                       tReal *B,
                                                       tReal& E) const {
        
        //get the number of particles of the network
        const tIndex& nParticles=network.getParticlesNumber();
 
        //get the dimension of the netwok
        const tDimension& dim=network.getDimension();
 
        //strength of the netweok
        const tReal& lambda=network.getLambda();
        //coordinates of the network
        const tReal *particles=network.getParticlesCoordinates().getValues();
        
        //iterator on Mu at point j
        const tReal *Mu_jk,*Mu_j=null;
        tReal lambda3=lambda*lambda*lambda;
 
        //iterators on dimension
        tDimension k;
 
        //iterators on particles
        tIndex i,j;
        
        //distance d
        tReal d3PiPjL3,d2PiPj,dd;
        
        //iterator on particle i
        const tReal *Pik,*Pi;
        const tReal *Mu_ik,*Mu_i=mu;
        
        //iterators on particle j
        const tReal *Pjk,*Pj=null;
        
        //<P^i-P^j,mu^j>
        tReal sRijMuj;
 
        //initializer
        tReal *Bik,*Bi=B;
 
        tReal Bdem_ik;
        
        //energy
        E=0;
        
        //coordinates of the particle i
        Pi=particles;
 
        //magnetized
       
        
        for (i=0;i<nParticles;i++) {//loop on i particle Pi
 
 
            //bdem_i=0
            memset(Bi,0,dim*sizeof(tReal));
           
            
           
            //magnetic moment of the particle j
            Mu_j=mu;
            //coordinates of the particle j
            Pj=particles;//first particle
            
            auto isMj=mAreMagnetizedCells.cbegin();
            for (j=0;j<nParticles;j++) { //loop on j particle
                if (j!=i) {
                    //test if the particule j is magnetized
                    if (*isMj) {
                        
                        //d2PiPj=|P_j-Pj|^2
                        d2PiPj=0;
                        //sXMu=<P^i-P^j,\mu^j>
                        sRijMuj=0;
                        //iterator
                        Pik=Pi;
                        Pjk=Pj;
                        Mu_jk=Mu_j;
                        for (k=0;k<dim;k++) {
                            //compute |Pi_Pj|^2
                            dd=(*Pik)-(*Pjk);
                            d2PiPj+=dd*dd;
                            //compute <Mu^j,Pi-Pj>
                            sRijMuj+=dd*(*Mu_jk);
                            
                            //next iterator
                            Pik++;
                            Pjk++;
                            Mu_jk++;
                        }
                        
                        //d3PiPj=|P_j-Pj|^3.\lambda^3
                        d3PiPjL3=(d2PiPj*sqrt(d2PiPj)*lambda3);
                        
                        Bik=Bi;
                        Mu_jk=Mu_j;
                        Pik=Pi;
                        Pjk=Pj;
                        Mu_ik=Mu_i;
                        for (k=0;k<dim;k++) {
                            
                            //Bd
                            Bdem_ik=-2*((*Mu_jk)-3.*(sRijMuj/d2PiPj)*((*Pik)-(*Pjk)) )/d3PiPjL3;
                            (*Bik)+=Bdem_ik;
                            
                            //energy
                            E-=0.5*Bdem_ik*(*Mu_ik);
                            
                            //next coordinate
                            Bik++;
                            Pik++;
                            Pjk++;
                            Mu_jk++;
                            Mu_ik++;
                        }
                    }//|mu_j|!=0
                } //i!=j
                
                
                //next iterator on j
                Mu_j+=dim;//magnetic moment at next particle j
                Pj+=dim;//coordinates of next particle j
                isMj++;//next particle is magnetized
            }//end loop on point j
                
 
            
            //coordinates of next particle i
            Pi+=dim;
            //mu at next particle i
            Mu_i+=dim;
            //Bd at next particle i
            Bi+=dim;
        }//end loop on particle i
        
        
    }//end method
    
  
    
    virtual  tReal computeEnergy(const tIndex& t,
                                 const SM_Network& network,
                                 const tReal *mu) const override{
        return computeTemplatedEnergy(t,network,mu);
    }
  
    inline tReal computeTemplatedEnergy(const tIndex& t,
                                        const SM_Network& network,
                                        const tReal *mu) const {
        
         //get the number of particles of the network
        const tIndex& nParticles=network.getParticlesNumber();
 
        //get the dimension of the netwok
        const tDimension& dim=network.getDimension();
 
        //strength of the network
        const tReal& lambda=network.getLambda();
        tReal iLambda3=1./(lambda*lambda*lambda);
 
        //coordinates of the network
        const tReal *particles=network.getParticlesCoordinates().getValues();
 
        
        //iterator on particles
        tIndex i,j;
        
        //iterator on dimension
        tDimension d;
 
        //Mu at particle i
        const tReal *Mui=mu,*Mui_d;
 
        //Mu at particle j
        const tReal *Muj,*Muj_d;
        
        //iterator on particle i
        const tReal *Pi_d,*Pi=particles;
 
        //iterators on particle j
        const tReal *Pj_d,*Pj=null;
 
        //rij=P_j-P_i
        tReal r[3];
        tReal *rij;
        //dij2=|P_j-P_i|^2 
        tReal dij2;
        //dij3=|P_j-P_i|^3 
        tReal dij3;
 
        //sMu=<Mu_i,Mu_j>
        //rMui=<rij,Mu_i>
        //rMuj=<rij,Mu_j>
        tReal sMu,rMui,rMuj;
            
        //energy
        tReal E=0;
        
        for (i=0;i<nParticles;i++) {//loop on particles i
 
            Pj=particles;
            Muj=mu;
            for (j=0;j<nParticles;j++) {//loop on particles j \=i
                if (j!=i) {
                    //rij=P_j-Pj
                    //dij=|rij|^2
                    rij=r;
                    dij2=0;
                    sMu=0;
                    rMui=0;
                    rMuj=0;
                    Pj_d=Pj;
                    Pi_d=Pi;
                    Mui_d=Mui;
                    Muj_d=Muj;
                    for (d=0;d<dim;d++) {
                        (*rij)=(*Pj_d)-(*Pi_d);
                        dij2+=(*rij)*(*rij);
                        sMu+= (*Mui_d)*(*Muj_d);
                        rMui+=(*Mui_d)*(*rij);
                        rMuj+=(*Muj_d)*(*rij);
                        Pj_d++;
                        Pi_d++;
                        rij++;
                        Mui_d++;
                        Muj_d++;
                    }
 
                    dij3=dij2*sqrt(dij2);
                
                    E+=iLambda3*(sMu-(3*rMui*rMuj/dij2))/dij3;
                }//end condition  i!=j
                      
                //next j-particles
                Muj+=dim;
                Pj+=dim;
                
            }//end loop on particle j
            
 
            //next i-particle
            Mui+=dim;
            Pi+=dim;
         
            
        }//end loop on particle i
 
 
        return E;
    }//end method
};
 
 
#endif