00001 #ifndef OPENTISSUE_DYNAMICS_SPH_SOLVERS_SPH_SURFACE_H 00002 #define OPENTISSUE_DYNAMICS_SPH_SOLVERS_SPH_SURFACE_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_solver.h> 00013 00014 namespace OpenTissue 00015 { 00016 namespace sph 00017 { 00018 //#define SPH_TENSION_THRESHOLD 0.000001 00019 //#define SPH_TENSION_THRESHOLD 0.001 00020 //#define SPH_TENSION_THRESHOLD 4.5 00021 //#define SPH_TENSION_THRESHOLD 56.25 // 2250/40 00022 //#define SPH_TENSION_THRESHOLD 12.5 // 500/40 00023 00027 template< typename Types, typename KernelPolicy > 00028 class SurfaceForce : Solver< Types, typename Types::vector > 00029 { 00030 public: 00031 typedef Solver<Types, typename Types::vector> base_type; 00032 typedef KernelPolicy smoothing_kernel; 00033 typedef typename base_type::value value; 00034 typedef typename Types::vector vector; 00035 typedef typename Types::real_type real_type; 00036 typedef typename Types::particle particle; 00037 typedef typename Types::particle_cptr_container::const_iterator particle_cptr_container_citerator; 00038 00039 public: 00043 SurfaceForce(const real_type& tension = 0, const real_type& threshold = 1) : base_type() 00044 , m_s(tension) 00045 , m_l(threshold) 00046 { 00047 } 00048 00052 ~SurfaceForce() 00053 { 00054 } 00055 00059 SurfaceForce& operator=(const SurfaceForce& rhs) 00060 { 00061 m_s = rhs.m_s; 00062 m_l = rhs.m_l; 00063 return *this; 00064 } 00065 00066 public: 00067 real_type& tension() 00068 { 00069 return m_s; 00070 } 00071 00072 const real_type& tension() const 00073 { 00074 return m_s; 00075 } 00076 00077 real_type& threshold() 00078 { 00079 return m_l; 00080 } 00081 00082 const real_type& threshold() const 00083 { 00084 return m_l; 00085 } 00086 00092 virtual value apply(const particle& par, particle_cptr_container_citerator begin, particle_cptr_container_citerator end) const 00093 { 00094 const vector& n = par.normal(); 00095 const real_type chk = n*n; 00096 if (chk < m_l) 00097 return value(0); 00098 real_type tmp = 0; 00099 for (particle_cptr_container_citerator p_ = begin; p_ != end; ++p_) { 00100 const particle* p = *p_; 00101 // if (&par == p) continue; // all sph forces apply i != j 00102 tmp += (p->mass()/p->density())*m_W.laplacian(par.position()-p->position()); 00103 } 00104 return value((-m_s*tmp)*normalize(n)); 00105 } 00106 00107 private: 00108 const smoothing_kernel m_W; 00109 real_type m_s; 00110 real_type m_l; 00111 00112 }; // End class SurfaceForce 00113 00114 } // namespace sph 00115 } // namespace OpenTissue 00116 00117 // OPENTISSUE_DYNAMICS_SPH_SOLVERS_SPH_SURFACE_H 00118 #endif
1.7.1