yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10695
[Branch ~yade-pkg/yade/git-trunk] Rev 3895: Implement some more kernel functions.
------------------------------------------------------------
revno: 3895
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-09 16:03:15 +0200
message:
Implement some more kernel functions.
modified:
pkg/common/SPHEngine.cpp
pkg/common/SPHEngine.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.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 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp 2014-04-09 09:11:49 +0000
+++ pkg/common/SPHEngine.cpp 2014-04-09 14:03:15 +0000
@@ -20,6 +20,23 @@
b->rho0 = rho0;
}
Real rho = 0;
+
+ // Pointer to kernel function
+ Real (*kernelFunctionCurrent)(const double & r, const double & h);
+ if (KernFunctionDensity==Poly6) {
+ kernelFunctionCurrent = smoothkernelPoly6;
+ } else if (KernFunctionDensity==Spiky) {
+ kernelFunctionCurrent = smoothkernelSpiky;
+ } else if (KernFunctionDensity==Visco) {
+ kernelFunctionCurrent = smoothkernelVisco;
+ } else if (KernFunctionDensity==Lucy) {
+ kernelFunctionCurrent = smoothkernelLucy;
+ } else if (KernFunctionDensity==Monaghan) {
+ kernelFunctionCurrent = smoothkernelMonaghan;
+ } else {
+ throw runtime_error("Kernel types can only have the following types: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5.");
+ }
+
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
const shared_ptr<Body> b2 = Body::byId((*it).first,scene);
Sphere* s=dynamic_cast<Sphere*>(b->shape.get());
@@ -36,8 +53,8 @@
const Real SmoothDist = -geom.penetrationDepth + phys.h;
- // [Mueller2003], (3)
- rho += b2->state->mass*smoothkernelPoly6(SmoothDist, phys.h);
+ // [Mueller2003], (3)
+ rho += b2->state->mass*kernelFunctionCurrent(SmoothDist, phys.h);
}
// [Mueller2003], (3), we need to consider the density of the current body (?)
rho += b->state->mass*smoothkernelPoly6(0.0, s->radius);
@@ -45,52 +62,137 @@
b->rho = rho;
}
-Real smoothkernelPoly6(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 315/(64*M_PI*pow(h,9))*pow((h*h - rrj*rrj), 3); // [Mueller2003], (20)
- } else {
- return 0;
- }
-}
-
-Real smoothkernelPoly6Grad(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 315/(64*M_PI*pow(h,9))*(-6*rrj)*pow((h*h - rrj*rrj), 2);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelPoly6Lapl(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 315/(64*M_PI*pow(h,9))*(-6)*(h*h - rrj*rrj)*(h*h - 5*rrj*rrj);
- } else {
- return 0;
- }
-}
-
-Real smoothkernelSpiky(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 15/(M_PI*pow(h,6))*pow((h - rrj), 3); // [Mueller2003], (21)
- } else {
- return 0;
- }
-}
-
-Real smoothkernelVisco(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 15/(2*M_PI*pow(h,3))*(-rrj*rrj*rrj/(2*h*h*h) + rrj*rrj/(h*h) + h/(2*rrj) - 1); // [Mueller2003], (21)
- } else {
- return 0;
- }
-}
-
-Real smoothkernelViscoLapl(const double & rrj, const double & h) {
- if (rrj<=h) {
- return 45/(M_PI*pow(h,6))*(h - rrj); // [Mueller2003], (22+)
- } else {
- return 0;
- }
+Real smoothkernelPoly6(const double & r, const double & h) {
+ if (r<=h) {
+ return 315/(64*M_PI*pow(h,9))*pow((h*h - r*r), 3);
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelPoly6Grad(const double & r, const double & h) {
+ if (r<=h) {;
+ return -315/(64*M_PI*pow(h,9))*(-6*r*pow((h*h-r*r), 2));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelPoly6Lapl(const double & r, const double & h) {
+ if (r<=h) {
+ return 315/(64*M_PI*pow(h,9))*(-6*(h*h-r*r)*(h*h - 5*r*r));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelSpiky(const double & r, const double & h) {
+ if (r<=h) {
+ return 15/(M_PI*pow(h,6))*(pow((h-r), 3)); // [Mueller2003], (21)
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelSpikyGrad(const double & r, const double & h) {
+ if (r<=h) {
+ return -15/(M_PI*pow(h,6))*(-3*pow((h-r),2));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelSpikyLapl(const double & r, const double & h) {
+ if (r<=h) {
+ return 15/(M_PI*pow(h,6))*(6*(h-r));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelVisco(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ return 15/(2*M_PI*pow(h,3))*(-r*r*r/(2*h*h*h) + r*r/(h*h) + h/(2*r) -1); // [Mueller2003], (21)
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelViscoGrad(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ return -15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelViscoLapl(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ //return 45/(M_PI*pow(h,6))*(h - rrj); // [Mueller2003], (22+)
+ return 15/(2*M_PI*pow(h,3))*(-3*r*r/(2*h*h*h) + 2*r/(h*h) - h/(2*r*r));
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelLucy(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ return 5/(9*M_PI*pow(h,2))*(1+3*r/h)*pow((1-r/h),3);
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelLucyGrad(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ return -5/(9*M_PI*pow(h,2))*(-12*r/(h*h))*pow((1-r/h),2);
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelLucyLapl(const double & r, const double & h) {
+ if (r<=h and r!=0 and h!=0) {
+ return 5/(9*M_PI*pow(h,2))*(-12/(h*h))*(1-r/h)*(1-3*r/h);
+ } else {
+ return 0;
+ }
+}
+
+Real smoothkernelMonaghan(const double & r, const double & h) {
+ Real ret = 0.0;
+ if (h!=0) {
+ if (r/h<0.5) {
+ ret = 40/(7*M_PI*h*h)*(1 - 6*pow((r/h),2) + 6*pow((r/h),3));
+ } else {
+ ret = 80/(7*M_PI*h*h)*pow((1 - (r/h)), 3);
+ }
+ }
+ return ret;
+}
+
+Real smoothkernelMonaghanGrad(const double & r, const double & h) {
+ Real ret = 0.0;
+ if (h!=0) {
+ if (r/h<0.5) {
+ ret = 40/(7*M_PI*h*h)*( - 6*r/(h*h))*(2 - 3 * r/(h*h*h));
+ } else {
+ ret = 80/(7*M_PI*h*h)*( -3/h)*pow((1 -r/h), 2);
+ }
+ }
+ return ret;
+}
+
+Real smoothkernelMonaghanLapl(const double & r, const double & h) {
+ Real ret = 0.0;
+ if (h!=0) {
+ if (r/h<0.5) {
+ ret = 40/(7*M_PI*h*h)*( - 12/(h*h))*(1 - 3 * r/(h*h*h*h*h));
+ } else {
+ ret = 80/(7*M_PI*h*h)*( 6/(h*h))*(1 -r/h);
+ }
+ }
+ return ret;
}
void computeForceSPH(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force) {
=== modified file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp 2014-04-09 09:11:49 +0000
+++ pkg/common/SPHEngine.hpp 2014-04-09 14:03:15 +0000
@@ -4,6 +4,7 @@
#include<yade/core/PartialEngine.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
+enum KernFunctions {Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5};
class SPHEngine: public PartialEngine{
public:
void calculateSPHRho(const shared_ptr<Body>& b);
@@ -12,17 +13,26 @@
((int, mask,-1,, "Bitmask for SPH-particles."))
((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)
+ ((KernFunctions,KernFunctionDensity, Poly6,, "Kernel function for density calculation (by default - Poly6). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
);
};
REGISTER_SERIALIZABLE(SPHEngine);
-Real smoothkernelPoly6(const double & rrj, const double & h); // [Mueller2003] (20)
-Real smoothkernelPoly6Grad(const double & rrj, const double & h);
-Real smoothkernelPoly6Lapl(const double & rrj, const double & h);
-Real smoothkernelSpiky(const double & rrj, const double & h); // [Mueller2003] (21)
-Real smoothkernelVisco(const double & rrj, const double & h); // [Mueller2003] (22)
-Real smoothkernelViscoLapl(const double & rrj, const double & h); // [Mueller2003] (22+)
+Real smoothkernelPoly6(const double & r, const double & h); // [Mueller2003] (20)
+Real smoothkernelPoly6Grad(const double & r, const double & h);
+Real smoothkernelPoly6Lapl(const double & r, const double & h);
+Real smoothkernelSpiky(const double & r, const double & h); // [Mueller2003] (21)
+Real smoothkernelSpikyGrad(const double & r, const double & h);
+Real smoothkernelSpikyLapl(const double & r, const double & h);
+Real smoothkernelVisco(const double & r, const double & h); // [Mueller2003] (22)
+Real smoothkernelViscoGrad(const double & r, const double & h); // [Mueller2003] (22)
+Real smoothkernelViscoLapl(const double & r, const double & h);
+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 smoothkernelMonaghan(const double & r, const double & h);
+Real smoothkernelMonaghanGrad(const double & r, const double & h);
+Real smoothkernelMonaghanLapl(const double & r, const double & h);
void computeForceSPH(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force);
-
#endif
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2014-04-09 09:11:49 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-04-09 14:03:15 +0000
@@ -240,6 +240,35 @@
phys->SPHmode=true;
phys->mu=(mat1->mu+mat2->mu)/2.0;
}
+
+ if (mat1->KernFunctionPressure==Poly6 and mat2->KernFunctionPressure==Poly6) {
+ phys->kernelFunctionCurrentPressure = smoothkernelPoly6Grad;
+ } else if (mat1->KernFunctionPressure==Spiky and mat2->KernFunctionPressure==Spiky) {
+ phys->kernelFunctionCurrentPressure = smoothkernelSpikyGrad;
+ } else if (mat1->KernFunctionPressure==Visco and mat2->KernFunctionPressure==Visco) {
+ phys->kernelFunctionCurrentPressure = smoothkernelViscoGrad;
+ } else if (mat1->KernFunctionPressure==Lucy and mat2->KernFunctionPressure==Lucy) {
+ phys->kernelFunctionCurrentPressure = smoothkernelLucyGrad;
+ } else if (mat1->KernFunctionPressure==Monaghan and mat2->KernFunctionPressure==Monaghan) {
+ phys->kernelFunctionCurrentPressure = smoothkernelMonaghanGrad;
+ } else {
+ throw runtime_error("Kernel types can only have the following types: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5. And for all contacting materials they should be equal.");
+ }
+
+ if (mat1->KernFunctionVisco==Poly6 and mat2->KernFunctionVisco==Poly6) {
+ phys->kernelFunctionCurrentVisco = smoothkernelPoly6Lapl;
+ } else if (mat1->KernFunctionVisco==Spiky and mat2->KernFunctionVisco==Spiky) {
+ phys->kernelFunctionCurrentVisco = smoothkernelSpikyLapl;
+ } else if (mat1->KernFunctionVisco==Visco and mat2->KernFunctionVisco==Visco) {
+ phys->kernelFunctionCurrentVisco = smoothkernelViscoLapl;
+ } else if (mat1->KernFunctionVisco==Lucy and mat2->KernFunctionVisco==Lucy) {
+ phys->kernelFunctionCurrentVisco = smoothkernelLucyLapl;
+ } else if (mat1->KernFunctionVisco==Monaghan and mat2->KernFunctionVisco==Monaghan) {
+ phys->kernelFunctionCurrentVisco = smoothkernelMonaghanLapl;
+ } else {
+ throw runtime_error("Kernel types can only have the following types: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5. And for all contacting materials they should be equal.");
+ }
+
#endif
}
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-04-09 09:11:49 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-04-09 14:03:15 +0000
@@ -35,6 +35,8 @@
#ifdef YADE_SPH
((bool,SPHmode,false,,"True, if SPH-mode is enabled."))
((Real,mu,-1,, "Viscosity. See Mueller [Mueller2003]_ .")) // [Mueller2003], (14)
+ ((KernFunctions,KernFunctionPressure,Spiky,, "Kernel function for pressure calculation (by default - Spiky). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
+ ((KernFunctions,KernFunctionVisco, Visco,, "Kernel function for viscosity calculation (by default - Visco). The following kernel functions are available: Poly6=1, Spiky=2, Visco=3, Lucy=4, Monaghan=5."))
#endif
((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_.")),
createIndex();
@@ -60,6 +62,8 @@
((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_")),
createIndex();
)
+ Real (*kernelFunctionCurrentPressure)(const double & r, const double & h);
+ Real (*kernelFunctionCurrentVisco)(const double & r, const double & h);
REGISTER_CLASS_INDEX(ViscElPhys,FrictPhys);
};
REGISTER_SERIALIZABLE(ViscElPhys);
