yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11565
[Branch ~yade-pkg/yade/git-trunk] Rev 3504: Add one more kernel function in SPH.
------------------------------------------------------------
revno: 3504
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
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);
