SM_ItoSystem.h

Go to the documentation of this file.

#ifndef SM_ItoSystem_H
#define SM_ItoSystem_H
 
#include "SM_TemplatedSystem.h"
 
class SM_ItoSystem : public SM_TemplatedSystem<SM_ItoSystem> {
   
    //attributes
private :
 
  
    
public:
    // CONSTRUCTORS 
    SM_ItoSystem(void) {
        
    }
    
    // DESTRUCTORS 
    virtual ~SM_ItoSystem(void) {
    }
 
 
public : 
   
 
 
public:
  
  
      
    inline static CORE_UniquePointer<SM_ItoSystem> New() {
        return CORE_UniquePointer<SM_ItoSystem>(new SM_ItoSystem(),SM_ItoSystem::Delete());
    }
 
     
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
 
    
public:
   
    inline tBoolean computeTemplatedMuAtNextTimeStep(const tReal& dt,const tReal& epsilon_t,
                                                     const tIndex nParticles,const tDimension& dim,
                                                     const tReal *mu_t,
                                                     tReal *mu_tpdt)  const {
        
        
        //inverse of norm of Mu
        tReal invNormMu=0;
 
        //iterator on point
        tIndex p;
        
        //iterator on dimension
        tDimension k;
        
        //iterator on mu(t+dt)
        tReal *mu_tpdt_p,*mu_tpdt_i=mu_tpdt;
        
        //iterator on mu(t)
        const tReal *mu_t_i=mu_t;
 
     
        for (p=0;p<nParticles;p++) {
            //mu(t+dt) at particule p
            mu_tpdt_p=mu_tpdt_i;
            //inverse of norm of mu
            invNormMu=0;
            
            for (k=0;k<dim;k++) {
                //mu(t+dt) at element i=(p,k)
                (*mu_tpdt_i)+=(*mu_t_i);
                //norm of mu(t+dt) at p
                invNormMu+=(*mu_tpdt_i)*(*mu_tpdt_i);
                
                //next element
                mu_tpdt_i++;
                mu_t_i++;
                    
            }
            if (invNormMu>EPSILON<tReal>) {
                //normalize mu(t+dt) at p
                invNormMu=1./sqrt(invNormMu);
                mu_tpdt_i=mu_tpdt_p;
                for (k=0;k<3;k++) {
                    (*mu_tpdt_i)*=invNormMu;
                    mu_tpdt_i++;
                }
            }
            
        }
        
        //because of nomalization the norm of mu is never greater than 1
        return true;
    }
    
};
 
 
#endif