------------------------------------------------------------ revno: 3504 committer: Anton Gladky <[email protected]> timestamp: Thu 2014-10-23 15:46:20 +0200 message: Add one more kernel function in SPH. modified: doc/references.bib pkg/common/SPHEngine.cpp pkg/common/SPHEngine.hpp
-- lp:yade https://code.launchpad.net/~yade-pkg/yade/git-trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'doc/references.bib' --- doc/references.bib 2014-10-20 13:42:59 +0000 +++ doc/references.bib 2014-10-23 13:46:20 +0000 @@ -839,3 +839,28 @@ doi = "http://dx.doi.org/10.1006/jcph.1997.5776";, url = "http://www.sciencedirect.com/science/article/pii/S0021999197957764";, } + +@article{Lucy1977, + author = {{Lucy}, L.~B.}, + title = "{A numerical approach to the testing of the fission hypothesis}", + journal = {\aj}, + year = 1977, + month = dec, + volume = 82, + pages = {1013-1024}, + doi = {10.1086/112164}, + url = {http://adsabs.harvard.edu/full/1977AJ.....82.1013L} +} + +@article{Monaghan1985, + author = {{Monaghan}, J.~J. and {Lattanzio}, J.~C.}, + title = "{A refined particle method for astrophysical problems}", + journal = {\aap}, + keywords = {Computational Astrophysics, Gravitational Collapse, Gravitational Fields, Many Body Problem, Molecular Clouds, Stellar Evolution, Angular Momentum, Binary Stars, Computational Grids, Interpolation, Kernel Functions, Particle Mass, Stellar Orbits}, + year = 1985, + month = aug, + volume = 149, + pages = {135-143}, + url = {http://adsabs.harvard.edu/abs/1985A%26A...149..135M} +} + === modified file 'pkg/common/SPHEngine.cpp' --- pkg/common/SPHEngine.cpp 2014-10-23 13:46:20 +0000 +++ pkg/common/SPHEngine.cpp 2014-10-23 13:46:20 +0000 @@ -53,8 +53,8 @@ } Real smoothkernelLucy(const double & r, const double & h) { - if (r<=h) { - // Lucy Kernel function, + if (r<=h && h>0) { + // Lucy Kernel function, [Lucy1977] (27) return 105./(16.*M_PI*h*h*h) * (1 + 3 * r / h) * pow((1.0 - r / h), 3); } else { return 0; @@ -62,8 +62,8 @@ } Real smoothkernelLucyGrad(const double & r, const double & h) { - if (r<=h) { - // 1st derivative of Lucy Kernel function, + if (r<=h && h>0) { + // 1st derivative of Lucy Kernel function, [Lucy1977] (27) return 105./(16.*M_PI*h*h*h) * (-12 * r) * (1/( h * h ) - 2 * r / (h * h * h ) + r * r / (h * h * h * h)); } else { return 0; @@ -71,14 +71,61 @@ } Real smoothkernelLucyLapl(const double & r, const double & h) { - if (r<=h) { - // 2nd derivative of Lucy Kernel function, + if (r<=h && h>0) { + // 2nd derivative of Lucy Kernel function, [Lucy1977] (27) return 105./(16.*M_PI*h*h*h) * (-12) / (h * h * h * h) * (h * h - 2 * r * h + 3 * r * r); } else { return 0; } } - +//========================================================================= + +Real smoothkernelBSpline(const double & r, const double & h) { + // BSpline Kernel function, [Monaghan1985] (22) + if (r<=2.0*h && h>0) { + const Real coefA = 3.0 / (4.0 * M_PI * h * h * h); + const Real r_h = r / h; + if (r<=h) { + return coefA * (10.0/3.0 - 7 * r_h * r_h + 4 * r_h * r_h * r_h); + } else { + return coefA * pow((2 - r_h), 2)*((5 - 4 * r_h) / 3.0); + } + } else { + return 0; + } +} + +Real smoothkernelBSplineGrad(const double & r, const double & h) { + // 1st derivative of BSpline Kernel function, [Monaghan1985] (22) + if (r<=2.0*h && h>0) { + const Real coefA = 3.0 / (4.0 * M_PI * h * h * h); + const Real r_h = r / h; + if (r<=h) { + return coefA * (-2) / (h * h) * ( 7 * r - 6 * r * r / h); + } else { + return coefA * 2 / h * (-6 + 7 * r / h - 2 * r * r / (h * h ) ); + } + } else { + return 0; + } +} + +Real smoothkernelBSplineLapl(const double & r, const double & h) { + // 2nd derivative of BSpline Kernel function, [Monaghan1985] (22) + if (r<=2.0*h && h>0) { + const Real coefA = 3.0 / (4.0 * M_PI * h * h * h); + const Real r_h = r / h; + if (r<=h) { + return coefA * (-2) / (h * h) * ( 7 - 12 * r_h); + } else { + return coefA * 2 / (h * h) * ( 7 - 4 * r_h); + } + } else { + return 0; + } +} + +//========================================================================= KernelFunction returnKernelFunction(const int a, const typeKernFunctions typeF) { return returnKernelFunction(a, a, typeF); } @@ -97,6 +144,16 @@ } else { KERNELFUNCDESCR } + } else if (a==BSpline) { + if (typeF==Norm) { + return smoothkernelBSpline; + } else if (typeF==Grad) { + return smoothkernelBSplineGrad; + } else if (typeF==Lapl) { + return smoothkernelBSplineLapl; + } else { + KERNELFUNCDESCR + } } else { KERNELFUNCDESCR } @@ -139,28 +196,32 @@ const Vector3r xixj = de2.pos - de1.pos; - if (xixj.norm() < phys.h) { - Real fpressure = 0.0; - if (Rho1!=0.0 and Rho2!=0.0) { - // from [Monaghan1992], (3.3), multiply by Mass2, because we need a force, not du/dt - fpressure = - Mass1 * Mass2 * ( - bodies[id1]->state->press/(Rho1*Rho1) + - bodies[id2]->state->press/(Rho2*Rho2) - ) - * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h); - } - - Vector3r fvisc = Vector3r::Zero(); - if (Rho1!=0.0 and Rho2!=0.0) { - // from [Morris1997], (22), multiply by Mass2, because we need a force, not du/dt - fvisc = phys.mu * Mass1 * Mass2 * - (-normalVelocity*geom.normal)/(Rho1*Rho2) * - 1 / (xixj.norm()) * - phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h); - //phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h); - } - force = fpressure*geom.normal + fvisc; - return true; + if ( phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h)) { + if (xixj.norm() < phys.h) { + Real fpressure = 0.0; + if (Rho1!=0.0 and Rho2!=0.0) { + // from [Monaghan1992], (3.3), multiply by Mass2, because we need a force, not du/dt + fpressure = - Mass1 * Mass2 * ( + bodies[id1]->state->press/(Rho1*Rho1) + + bodies[id2]->state->press/(Rho2*Rho2) + ) + * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h); + } + + Vector3r fvisc = Vector3r::Zero(); + if (Rho1!=0.0 and Rho2!=0.0) { + // from [Morris1997], (22), multiply by Mass2, because we need a force, not du/dt + fvisc = phys.mu * Mass1 * Mass2 * + (normalVelocity*geom.normal)/(Rho1*Rho2) * + 1 / (xixj.norm()) * + phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h); + //phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h); + } + force = fpressure*geom.normal + fvisc; + return true; + } else { + return false; + } } else { return false; } === modified file 'pkg/common/SPHEngine.hpp' --- pkg/common/SPHEngine.hpp 2014-10-17 17:14:29 +0000 +++ pkg/common/SPHEngine.hpp 2014-10-23 13:46:20 +0000 @@ -6,8 +6,8 @@ typedef Real (* KernelFunction)(const double & r, const double & h); -enum KernFunctions {Lucy=1}; -#define KERNELFUNCDESCR throw runtime_error("Type of kernel function undefined! The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."); +enum KernFunctions {Lucy=1,BSpline=2}; +#define KERNELFUNCDESCR throw runtime_error("Type of kernel function undefined! The following kernel functions are available: Lucy=1, BSpline=2."); enum typeKernFunctions {Norm, Grad, Lapl}; class SPHEngine: public PartialEngine{ @@ -19,13 +19,16 @@ ((Real,k,-1,, "Gas constant for SPH-interactions (only for SPH-model). See Mueller [Mueller2003]_ .")) // [Mueller2003], (11) ((Real,rho0,-1,, "Rest density. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1) ((Real,h,-1,, "Core radius. See Mueller [Mueller2003]_ .")) // [Mueller2003], (1) - ((int,KernFunctionDensity, Lucy,, "Kernel function for density calculation (by default - Lucy). The following kernel functions are available: Lucy=1.")) + ((int,KernFunctionDensity, Lucy,, "Kernel function for density calculation (by default - Lucy). The following kernel functions are available: Lucy=1, BSpline=2.")) ); }; REGISTER_SERIALIZABLE(SPHEngine); Real smoothkernelLucy(const double & r, const double & h); Real smoothkernelLucyGrad(const double & r, const double & h); Real smoothkernelLucyLapl(const double & r, const double & h); +Real smoothkernelBSpline(const double & r, const double & h); +Real smoothkernelBSplineGrad(const double & r, const double & h); +Real smoothkernelBSplineLapl(const double & r, const double & h); KernelFunction returnKernelFunction(const int a, const int b, const typeKernFunctions typeF); KernelFunction returnKernelFunction(const int a, const typeKernFunctions typeF);
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
