00001 #ifndef OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_VISCOSITY_H 00002 #define OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_VISCOSITY_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_kernel.h> 00013 #include <OpenTissue/core/math/math_constants.h> 00014 00015 namespace OpenTissue 00016 { 00017 namespace sph 00018 { 00019 00025 template< typename Types, 00026 struct OpenTissue::utility::RuntimeType<typename Types::real_type>* Radius, 00027 bool CheckRange > 00028 class WViscosity : public FixedSmoothingKernel<Types, CheckRange> 00029 { 00030 public: 00031 typedef FixedSmoothingKernel<Types, CheckRange> base_type; 00032 typedef typename Types::real_type real_type; 00033 typedef typename Types::vector vector; 00034 public: 00038 WViscosity() : base_type(*Radius) 00039 { 00040 m_inv_hSqr = 1./base_type::m_radiusSqr; 00041 m_inv_2hTri = m_inv_hSqr/(2.*base_type::m_radius); 00042 m_k = (15./math::detail::pi<real_type>())*m_inv_2hTri; 00043 m_l = (15./math::detail::pi<real_type>())*m_inv_2hTri; 00044 m_m = 45./(math::detail::pi<real_type>()*pow(base_type::m_radius, 6)); 00045 } 00046 00047 public: 00048 00052 real_type evaluate(const vector& r) const 00053 { 00054 if (!checkRange(r)) 00055 return 0.; 00056 register real_type tmp = length(r); 00057 if (tmp <= 0.) 00058 return 1e38; // +oo 00059 tmp = -m_inv_2hTri*tmp*tmp*tmp + m_inv_hSqr*tmp*tmp + (base_type::m_radius/(2.*tmp)) - 1.; 00060 return real_type(m_k*tmp); 00061 } 00062 00066 vector gradient(const vector& r) const 00067 { 00068 if (!checkRange(r)) 00069 return vector(0); 00070 register real_type tmp = length(r); 00071 if (tmp <= 0.) { 00072 vector rnd; 00073 random(rnd,-0.0001, 0.0001); 00074 return vector(-1e38*rnd); // -oo 00075 } 00076 tmp = -3*tmp*m_inv_2hTri + 2/base_type::m_radiusSqr - base_type::m_radiusSqr/(2*tmp*tmp*tmp); 00077 return vector((m_l*tmp)*r); 00078 } 00079 00083 real_type laplacian(const vector& r) const 00084 { 00085 if (!checkRange(r)) 00086 return 0.; 00087 return real_type(m_m*(base_type::m_radius-length(r))); 00088 } 00089 00090 protected: 00091 real_type m_k; 00092 real_type m_l; 00093 real_type m_m; 00094 real_type m_inv_hSqr; 00095 real_type m_inv_2hTri; 00096 00097 }; // End class WViscosity 00098 00099 } // namespace sph 00100 } // namespace OpenTissue 00101 00102 // OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_VISCOSITY_H 00103 #endif
1.7.1