00001 #ifndef OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_GAUSSIAN_H 00002 #define OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_GAUSSIAN_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 { 00033 template< typename Types, 00034 struct OpenTissue::utility::RuntimeType<typename Types::real_type>* Radius, 00035 bool CheckRange > 00036 class WFixedGaussian : public FixedSmoothingKernel<Types, CheckRange> 00037 { 00038 public: 00039 typedef FixedSmoothingKernel<Types, CheckRange> base_type; 00040 typedef typename Types::real_type real_type; 00041 typedef typename Types::vector vector; 00042 public: 00046 WFixedGaussian() : base_type(*Radius) 00047 { 00048 m_inv_2hSqr = 1./(2.*this->m_radiusSqr); 00049 m_3hSqr = 3.*this->m_radiusSqr; 00050 00051 m_k = 1./pow(2*math::detail::pi<real_type>()*this->m_radiusSqr, 1.5); 00052 m_l = sqrt(2.)/(4.*math::detail::pi<real_type>()*this->m_radiusSqr*this->m_radiusSqr*sqrt(math::detail::pi<real_type>()*this->m_radiusSqr)); 00053 m_m = sqrt(2.)/(4.*math::detail::pi<real_type>()*this->m_radiusSqr*this->m_radiusSqr*this->m_radiusSqr*sqrt(math::detail::pi<real_type>()*this->m_radiusSqr)); 00054 } 00055 00059 ~WFixedGaussian() 00060 { 00061 } 00062 00063 public: 00067 real_type evaluate(const vector& r) const 00068 { 00069 real_type res = 0; 00070 if (!checkRange(r)) 00071 return res; 00072 res = (r*r)*m_inv_2hSqr; 00073 res = m_k*::exp(-res); 00074 return res; 00075 } 00076 00080 vector gradient(const vector& r) const 00081 { 00082 vector res; 00083 res *= 0; 00084 if (!checkRange(r)) 00085 return res; 00086 real_type tmp = (r*r)*m_inv_2hSqr; 00087 tmp = m_l*::exp(-tmp); 00088 res = tmp*r; 00089 return res; 00090 } 00091 00095 real_type laplacian(const vector& r) const 00096 { 00097 real_type res = 0; 00098 if (!checkRange(r)) 00099 return res; 00100 const real_type tmp = r*r; 00101 res = tmp*m_inv_2hSqr; 00102 res = m_m*::exp(-res); 00103 res *= tmp-m_3hSqr; 00104 return res; 00105 } 00106 protected: 00107 real_type m_k; 00108 real_type m_l; 00109 real_type m_m; 00110 real_type m_inv_2hSqr; 00111 real_type m_3hSqr; 00112 00113 }; // End class WFixedGaussian 00114 00115 } // namespace sph 00116 } // namespace OpenTissue 00117 00118 // OPENTISSUE_DYNAMICS_SPH_KERNELS_SPH_GAUSSIAN_H 00119 #endif