SM_UniaxialAnisotropyOperator.h

#ifndef SM_UniaxialAnisotropyOperator_H
#define SM_UniaxialAnisotropyOperator_H
 
//base classes
#include "SM_AnisotropyOperator.h"
#include "SM_LinearSliceOperator.h"
 
//network class
#include "SM_Network.h"
 
 
 
 
 
class SM_UniaxialAnisotropyOperator : public SM_AnisotropyOperator<SM_UniaxialAnisotropyOperator>,public SM_LinearSliceOperator {
    
    //attributes
private :
 
    
protected:
    // CONSTRUCTORS 
    SM_UniaxialAnisotropyOperator(void) {
        setName("UniaxialAnisotropy");
    }
    
    // DESTRUCTORS 
    virtual ~SM_UniaxialAnisotropyOperator(void) {
    }
    
public : 
 
    //Instance building
    //=================
 
    inline static CORE_UniquePointer<SM_UniaxialAnisotropyOperator> New() {
        return CORE_UniquePointer<SM_UniaxialAnisotropyOperator>(new SM_UniaxialAnisotropyOperator(),
                                                                 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_UniaxialAnisotropyOperator>::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;
        phi*=-1;
    }
    
    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++;
        }
      
        //restore U pointer
        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;
        phi*=-1;
    }
     
    
   
    
};
 
 
#endif