AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

odeSolver.inl

Go to the documentation of this file.

#ifndef animal_odesolver_inl___________________________
#define animal_odesolver_inl___________________________


#include "odeSolver.h"
#include <iostream>



namespace animal
{

template<class P,class V,class R>
OdeSolver<P,V,R>::OdeSolver()
    : _method( EULER )
    , _MMIDsteps(3)
    , d1(dv1,da1)
    , d2(dv2,da2)
    , d3(dv3,da3)
    , d4(dv4,da4)
    , s1(sp1,sv1)
    , s2(sp2,sv2)
    , s3(sp3,sv3)
{}

template<class P,class V,class R>
void OdeSolver<P,V,R>::solveODE(
    Positions& pos,
    Vector& vel,
    Real dt
)
{
    //std::cerr<<"apply method " << _method << std::endl;
    if( _method == EULER ){
        //std::cerr<<"OdeSolver<P,V,R>::solveODE apply euler" << std::endl;
        euler(pos,vel,dt);
    }
    else if( _method == RK2 ){
        //std::cerr<<"apply rk2" << std::endl;
        rk2(pos,vel,dt);
    }
    else if( _method == RK4 ){
        //std::cerr<<"apply rk4" << std::endl;
        rk4(pos,vel,dt);
    }
    else if( _method == MODMID ){
        //std::cerr<<"apply modmid" << std::endl;
        modmid(pos,vel,dt,_MMIDsteps);
    }
    else if( _method == VVERLET ){
        //std::cerr<<"apply velocity verlet" << std::endl;
        velocityVerlet(pos,vel,dt);
    }
    else std::cerr <<"unknown integration method "<<_method<< std::endl;
}



template<class P,class V,class R>
void OdeSolver<P,V,R>::setMethod( int m )
{
    _method = m;
}

template<class P,class V,class R>
int OdeSolver<P,V,R>::method() const
{
    return _method;
}

//template<class P,class V,class R>
//void OdeSolver<P,V,R>::setMMIDsteps( int n )
//{
//  _MMIDsteps = n;
//}

//template<class P,class V,class R>
//int OdeSolver<P,V,R>::MMIDsteps() const
//{
//  return _MMIDsteps;
//}

template<class P,class V,class R>
void OdeSolver<P,V,R>::computeDerivative(
    Der& d,
    const Sta& s,
    Real
)
{
            computeAccelerations( d.accelerations, s.positions, s.velocities );
            v_eq( d.velocities, s.velocities );

}

template<class P,class V,class R>
void OdeSolver<P,V,R>::euler(Positions& pos, Vector& vel, Real h, Real t)
{
    d1.resize( size(pos) );
    Sta states(pos,vel);
//  cerr<<"OdeSolver<P,V,R>::euler: states.pos= "<< v_output(states.positions) <<endl;
//  cerr<<"OdeSolver<P,V,R>::euler: states.vel= "<< v_output(states.velocities) <<endl;
    integrate_euler(states,h,t,d1);
}

template<class P,class V,class R>
void OdeSolver<P,V,R>::rk2(Positions& pos, Vector& vel, Real h, Real t)
{
    d1.resize( size(pos) );
    s1.resize( size(pos) );
    States<P,V> states(pos,vel);
    integrate_rk2(states,h,t,d1,s1);
}

template<class P,class V,class R>
void OdeSolver<P,V,R>::rk4(Positions& pos, Vector& vel, Real h, Real t)
{
    d1.resize( size(pos) );
    d2.resize( size(pos) );
    d3.resize( size(pos) );
    d4.resize( size(pos) );
    s1.resize( size(pos) );
    States<P,V> states(pos,vel);
    integrate_rk4(states,h,t,d1,d2,d3,d4,s1);
}

template<class P,class V,class R>
void OdeSolver<P,V,R>::modmid(Positions& pos, Vector& vel, Real h, int n, Real t)
{
    d1.resize( size(pos) );
    s1.resize( size(pos) );
    s2.resize( size(pos) );
    s3.resize( size(pos) );
    States<P,V> states(pos,vel);
    integrate_modmid(states,h,t,n,d1,s1,s2,s3);
}

template<class P, class V, class R>
void OdeSolver<P, V, R>::velocityVerlet(Positions& pos, Vector& vel, Real h)
{
    // auxiliary variable for half step velocities - need to make sure its the right size
    velHalf.resize(size(vel));

    // derivatives and states storage
    d1.resize(size(pos));
    States<P,V> states(pos, vel);

    integrate_VVerlet(states,h,d1);

}
//=======================================================================
template<class P,class V,class R>
void OdeSolver<P,V,R>::integrate_euler
(
    Sta& s,
    Real h,
    Real t,
    Der& d
)
{
    computeDerivative(d,s,t);
    v_eq_euler_step(s,s,h,d);
}

template<class P,class V,class R>
void OdeSolver<P,V,R>::integrate_rk2
( 
    Sta& s, 
    Real h, 
    Real t, 
    Der& d, 
    Sta& s1 
)
{
    computeDerivative(d,s,t);
    v_eq_euler_step(s1,s,h/2,d);
    computeDerivative(d,s1,t+(h/2));
    v_eq_euler_step(s,s,h,d);
}


template<class P,class V,class R>
void OdeSolver<P,V,R>::integrate_rk4
(
    Sta& s, 
    Real h,
    Real t, 
    Der& d1, 
    Der& d2, 
    Der& d3, 
    Der& d4, 
    Sta& s1 
)
{
    computeDerivative(d1,s,t);
    v_eq_euler_step(s1,s,h/2,d1);
    computeDerivative(d2,s1,t+(h/2));
    v_eq_euler_step(s1,s,h/2,d2);
    computeDerivative(d3,s1,t+(h/2));
    v_eq_euler_step(s1,s,h,d3);
    computeDerivative(d4,s1,t+h);
    
    // accumulate composed derivative in d2
    v_peq(d2,d3); // d2+d3
    v_teq(d2,2.);  // 2(d2+d3)
    v_peq(d2,d1); // d1+2(d2+d3)
    v_peq(d2,d4); // d1+2(d2+d3)+d4
    
    v_eq_euler_step(s,s,h/6,d2);
}


template<class P,class V,class R>
void OdeSolver<P,V,R>::integrate_modmid
( 
    Sta& s, 
    Real h, 
    Real t, 
    int n, 
    Der& d1, 
    Sta& s1, 
    Sta& s2, 
    Sta& s3 
)
{
    Real hsub = h/n;       // h: total step to be made
    Real hhsub = 2.0*hsub; // hsub: substep value

    v_eq(s1,s);

    computeDerivative(d1,s,t);
    v_eq_euler_step(s2,s,hsub,d1); // First substep taken
    
    // take remaining substeps minus one
    int i;
    for( i=0; i<n-1; ++i )
    {
        computeDerivative(d1,s2,t+i*hsub);
        v_eq_euler_step(s3,s1,hhsub,d1);
        v_eq(s1,s2);
        v_eq(s2,s3);
    }
    
    //  Take last substep
    computeDerivative(d1,s2,t+i*hsub);
    v_eq_euler_step(s3,s1,hhsub,d1);
    
    // Take another step
    computeDerivative(d1,s3,t+h);
    v_eq_euler_step(s3,s3,hsub,d1);
    
    // result = average(s2,s3)
    v_peq(s2,s3);
    Real unDemi = 0.5;
    v_eq_ab( s,unDemi,s2 );
}

template<class P, class V, class R>
void OdeSolver<P, V, R>::integrate_VVerlet(
    Sta& s,
    Real h,
    Der& d)
{
    // get accelerations for this time step
    computeAccelerations(d.accelerations, s.positions, s.velocities);

    // velocity at half timestep is...
    v_eq(velHalf, s.velocities, h/2, d.accelerations);

    // the updated position at time t+Delta_t is...
    v_eq(s.positions, s.positions, h, velHalf);

    // acceleration at time t+Delta_t based on updated positions is...
    // note: we are using halfstep velocities to predict acceleration at update time here
    computeAccelerations(d.accelerations, s.positions, velHalf);

    // and finally, the velocity at time t+Delta_t is...
    v_eq(s.velocities, velHalf, h/2, d.accelerations);

}

}

#endif

---