00001 #ifndef OPENTISSUE_DYNAMICS_SPH_INTEGRATORS_SPH_LEAP_FROG_H 00002 #define OPENTISSUE_DYNAMICS_SPH_INTEGRATORS_SPH_LEAP_FROG_H 00003 // 00004 // OpenTissue Template Library 00005 // - A generic toolbox for physics-based modeling and simulation. 00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen. 00007 // 00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php 00009 // 00010 #include <OpenTissue/configuration.h> 00011 00012 #include <OpenTissue/dynamics/sph/sph_integrator.h> 00013 00014 namespace OpenTissue 00015 { 00016 namespace sph 00017 { 00018 00033 template< typename Types > 00034 class LeapFrog : public Integrator<Types> 00035 { 00036 public: 00037 typedef Integrator<Types> base_type; 00038 typedef typename Types::real_type real_type; 00039 typedef typename Types::vector vector; 00040 typedef typename Types::particle particle; 00041 typedef typename Types::particle_container particle_container; 00042 typedef typename base_type::collision_detection collision_detection; 00043 public: 00047 LeapFrog(const real_type& timestep, const real_type& restitution) 00048 : base_type(timestep, restitution) 00049 {} 00050 00054 ~LeapFrog() 00055 {} 00056 00057 public: 00058 00062 virtual void initialize(particle& par) const 00063 { 00064 vector &a0 = par.acceleration(); 00065 vector &v_apx = par.velocity(); 00066 vector &v0 = par.position_old(); 00067 00068 // compute a at time 0 00069 a0 = par.force()/par.density(); 00070 00071 // compute v at time - 00072 v0 = v_apx - .5*base_type::m_dt*a0; 00073 } 00074 00080 void integrate(particle& par) const 00081 { 00082 vector &a = par.acceleration(); // a at time t 00083 vector &v_apx = par.velocity(); // v at time t 00084 vector &v = par.position_old(); // v at time t - 00085 vector v_old = v; // v at time t - 00086 vector &x = par.position(); // x at time t 00087 00088 // compute a at time t 00089 a = par.force()/par.density(); 00090 00091 // compute v at time n + 00092 v += base_type::m_dt*a; 00093 00094 // compute x at time n + 1, using v at time n + 00095 x += base_type::m_dt*v; 00096 00097 // COLLISION 00098 if (base_type::m_colisys) { 00099 typename collision_detection::collision_type coli; 00100 if (base_type::m_colisys->collision(coli, par)) { 00101 x = coli.contact(); 00102 const real_type r = base_type::m_r>0?base_type::m_r*coli.penetration()/(base_type::m_dt*(length(v))):0; 00103 v -= (1+r)*(v*coli.normal())*coli.normal(); 00104 // v_apx += coli.normal()*(v_apx*coli.normal()); 00105 // v -= (1+base_type::m_r)*(v*coli.normal())*coli.normal(); 00106 // v -= (v*coli.normal())*coli.normal(); 00107 // real_type l = sqr_length(base_type::m_dt*v); 00108 // if (l>0) l = 1./l; 00109 // v -= (1.+coli.penetration()*coli.penetration()*l)*(v*coli.normal())*coli.normal(); 00110 //v -= ((1+m_r)*(v*n))*n; 00111 //v_apx -= ((1+m_r)*(v_apx*n))*n; 00112 } 00113 } 00114 00115 // estimate velocity at time t 00116 v_apx = .5*(v_old+v); 00117 00118 } 00119 }; // End class LeapFrog 00120 00121 } // namespace sph 00122 } // namespace OpenTissue 00123 00124 // OPENTISSUE_DYNAMICS_SPH_INTEGRATORS_SPH_LEAP_FROG_H 00125 #endif
1.7.1