SM_PlanarAnisotropyOperator.h

#ifndef SM_PlanarAnisotropyOperator_H
#define SM_PlanarAnisotropyOperator_H
 
//base classes
#include "SM_AnisotropyOperator.h"
#include "SM_LinearSliceOperator.h"
 
//network class
#include "SM_Network.h"
 
 
 
class SM_PlanarAnisotropyOperator : public SM_AnisotropyOperator<SM_PlanarAnisotropyOperator>,public SM_LinearSliceOperator {
    
    //attributes
private :
 
    
protected:
    // CONSTRUCTORS 
    SM_PlanarAnisotropyOperator(void) {
        setName("PlanarAnisotropy");
    }
    
    // DESTRUCTORS 
    virtual ~SM_PlanarAnisotropyOperator(void) {
    }
    
public : 
 
    //Instance building
    //=================
 
    inline static CORE_UniquePointer<SM_PlanarAnisotropyOperator> New() {
        return CORE_UniquePointer<SM_PlanarAnisotropyOperator>(new SM_PlanarAnisotropyOperator(),
                                                               SM_Operator::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_PlanarAnisotropyOperator>::getContentsMemorySize();
        mem+=SM_LinearSliceOperator::getContentsMemorySize();
        return mem;
    }
    
    inline void computePhi(const tReal *U,const tReal *eU,
                           const tReal *S,const tReal * eS,
                           tReal& phi,tReal& W) const {
        //phi=<S,U>
        phi=0;
        while (S!=eS) {
            phi+=(*S)*(*U);
            S++;
            U++;
        }
      
        
        //phi=<S,U>^2
        phi*=phi;
    }
 
     
    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 {
        //phi=<S,U>
        phi=0;
        while (S!=eS) {
            phi+=(*S)*(*U);
            S++;
            U++;
        }
       
        //S++;
        //U++;
        //restore U
        U-=SM_Constants::DIM;
        
       
        //gradPhi=2<S,U>.U
        for(auto& gradPhi_k:gradPhi) {
            gradPhi_k=2*(*U)*phi;
            U++;
        }
        
        
        //phi=<S,U>^2
        phi*=phi;
        
        
    }
     
    
   
    
};
 
 
#endif