SM_CubicAnisotropyOperator.h

#ifndef SM_CubicAnisotropyOperator_H
#define SM_CubicAnisotropyOperator_H
 
//base classes
#include "SM_AnisotropyOperator.h"
 
//network class
#include "SM_Network.h"
 
 
 
class SM_CubicAnisotropyOperator : public SM_AnisotropyOperator<SM_CubicAnisotropyOperator> {
    
    //attributes
private :
 
   
    
protected:
    // CONSTRUCTORS 
    SM_CubicAnisotropyOperator(void);
    
    // DESTRUCTORS 
    virtual ~SM_CubicAnisotropyOperator(void);
public : 
 
    //Instance building
    //=================
 
    inline static CORE_UniquePointer<SM_CubicAnisotropyOperator> New() {
        return CORE_UniquePointer<SM_CubicAnisotropyOperator>(new SM_CubicAnisotropyOperator(),
                                                              SM_CubicAnisotropyOperator::Delete());
    }
    
    virtual CORE_UniquePointer<SM_SliceOperator> NewInstance() const override {
        return New();
    }
public:
 
     //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
  
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SM_AnisotropyOperator<SM_CubicAnisotropyOperator>::getContentsMemorySize();
        return mem;
    }
    inline void computePhi(const tReal *U,const tReal* eU,
                           const tReal*S,const tReal *eS,
                           tReal& phi,tReal& W) const {
        // OMP_CRITICAL(cout) {
        //     std::cout<<"U size:"<<(eU-U)<<" U={";
        //     const tReal *iU=U;
        //     while (iU!=eU) {
        //         std::cout<<(*iU)<<" ";
        //         iU++;
        //     }
        //     std::cout<<"}\n";
        // }
       
        phi=0;
        while (U!=eU) {
            
            //W=<S,Ud>
            W=0;
            while (S!=eS) {
 
                W+=(*S)*(*U);
                
                //next coordinate
                S++;
                U++;
            }
          
            
            //W=<S,Ud>^2
            W*=W;
            //W=<S,Ud>^4
            W*=W;
          
            // \Phi_U(S)= \sum_{d=0}^{d=2} .<S,U_d>^4 
            phi+=W;
 
            //restore S pointer
            S-=SM_Constants::DIM;
            
           
            
        }//end loop on direction d 
 
        // \Phi_U(S)= 0.5 \sum_{d=0}^{d=2} .<S,U_d>^4 
        phi*=0.5;
       
       
    }
    
    inline void  computePhiDerivatives(const tReal *U,const tReal* eU,
                                       const tReal*S,const tReal *eS,
                                       tReal& phi,std::array<tReal,SM_Constants::DIM>& gradPhi,
                                       tReal& W,tReal& W3) const {
 
        //gradPhi=0;
        memset(gradPhi.data(),0,sizeof(tReal)*SM_Constants::DIM);
        
        phi=0;
        while (U!=eU) {//loop on anistropy direction
            //W=<S,Ud>
            W=0;
            while (S!=eS) {
                W+=(*S)*(*U);
                //next coordinate
                S++;
                U++;
            }
           
            //restore S pointer
            S-=SM_Constants::DIM;
            
            //restore U pointer
            U-=SM_Constants::DIM;
            
            //W3=<S,Ud>^3
            W3=W;
            W3*=W3;
            W3*=W;
         
            //gradPhi+=2*Ud.<S,U_d>^3
            for(auto& gradPhi_k:gradPhi) {   
                gradPhi_k+=2*W3*(*U);
                U++;
            }
           
            //W^4
            W3*=W;
            
            //\Phi_U(S)= sum_{d=0}^{d=2} .<S,U_d>^4 
            phi+=W3;
            
            //next anisotropy direction
            
        }//end loop on direction
        //\Phi_U(S)= - \frac{1}{2} \sum_{d=0}^{d=2} .<S,U_d>^4 
        phi*=0.5;
        
        
        
    }
     
    
   
    
};
 
 
#endif