yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3894
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-09 11:11:49 +0200
message:
  Add gradients and laplacian of the kernel functions.
modified:
  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 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp	2014-04-09 09:11:49 +0000
+++ pkg/common/SPHEngine.cpp	2014-04-09 09:11:49 +0000
@@ -35,9 +35,11 @@
       if (Mass == 0) Mass = b->state->mass;
       
       const Real SmoothDist = -geom.penetrationDepth + phys.h;
-      
+     
+      // [Mueller2003], (3) 
       rho += b2->state->mass*smoothkernelPoly6(SmoothDist, phys.h);
     }
+    // [Mueller2003], (3), we need to consider the density of the current body (?)
     rho += b->state->mass*smoothkernelPoly6(0.0, s->radius);
   }
   b->rho = rho;
@@ -51,6 +53,22 @@
   }
 }
 
+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)
@@ -125,7 +143,7 @@
   
   Real Rho = bodies[id2]->rho;
   if (Rho==0.0 and bodies[id1]->rho!=0.0) {
-    Rho = bodies[id1]->rho!=0.0;
+    Rho = bodies[id1]->rho;
   }
   
   const Vector3r xixj = (-geom.penetrationDepth + phys.h)*geom.normal;
@@ -133,18 +151,19 @@
   if (xixj.norm() < phys.h) {
     Real fpressure = 0.0;
     if (Rho!=0.0) {
-      fpressure = Mass * 
+      // [Mueller2003], (10)
+      fpressure = -Mass * 
                   (bodies[id1]->press + bodies[id2]->press)/(2.0*Rho) *
-                  smoothkernelSpiky(xixj.norm(), phys.h);
+                  smoothkernelPoly6Grad(xixj.norm(), phys.h);
     }
     
     Vector3r fvisc = Vector3r::Zero();
     if (Rho!=0.0) {
       fvisc = phys.mu * Mass * 
                   normalVelocity*geom.normal/Rho *
-                  smoothkernelViscoLapl(xixj.norm(), phys.h);
+                  smoothkernelPoly6Lapl(xixj.norm(), phys.h);
     }
-    force = fpressure*geom.normal - fvisc;
+    force = fpressure*geom.normal + fvisc;
     return;
   } else {
     return;

=== modified file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp	2014-04-09 09:11:49 +0000
+++ pkg/common/SPHEngine.hpp	2014-04-09 09:11:49 +0000
@@ -16,6 +16,8 @@
 };
 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+)