00001 #ifndef OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_SPIKY_H 00002 #define OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_SPIKY_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 00024 template< typename Types, 00025 struct OpenTissue::utility::RuntimeType<typename Types::real_type>* Radius, 00026 bool CheckRange > 00027 class WSpiky : public FixedSmoothingKernel<Types, CheckRange> 00028 { 00029 public: 00030 typedef FixedSmoothingKernel<Types, CheckRange> base_type; 00031 typedef typename Types::real_type real_type; 00032 typedef typename Types::vector vector; 00033 public: 00037 WSpiky() : base_type(*Radius) 00038 { 00039 m_k = 15./(math::detail::pi<real_type>()*pow(base_type::m_radius, 6)); 00040 m_l = -45./(math::detail::pi<real_type>()*pow(base_type::m_radius, 6)); 00041 m_m = -90./(math::detail::pi<real_type>()*pow(base_type::m_radius, 6)); 00042 } 00043 00044 public: 00045 00049 real_type evaluate(const vector& r) const 00050 { 00051 if (!checkRange(r)) 00052 return 0.; 00053 register real_type res = base_type::m_radius- length(r); 00054 res *= res*res*m_k; 00055 return res; 00056 } 00057 00061 vector gradient(const vector& r) const 00062 { 00063 if (!checkRange(r)) 00064 return vector(0); 00065 register real_type tmp = length(r); 00066 if (tmp <= 0.) { 00067 vector rnd; 00068 random(rnd,-0.0001, 0.0001); 00069 return vector(m_l*rnd); 00070 } 00071 const register real_type tmp2 = base_type::m_radius-tmp; 00072 tmp = (tmp2*tmp2)/tmp; 00073 return vector((m_l*tmp)*r); 00074 } 00075 00079 real_type laplacian(const vector& r) const 00080 { 00081 if (!checkRange(r)) 00082 return 0.; 00083 register real_type tmp = length(r); 00084 if (tmp <= 0.) 00085 return -1e38; // -oo 00086 tmp = ((base_type::m_radius-tmp)*(base_type::m_radius-2*tmp))/tmp; 00087 return real_type(m_m*tmp); 00088 } 00089 00090 protected: 00091 real_type m_k; 00092 real_type m_l; 00093 real_type m_m; 00094 00095 }; // End class WSpiky 00096 00097 } // namespace sph 00098 } // namespace OpenTissue 00099 00100 // OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_SPIKY_H 00101 #endif
1.7.1