yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3518
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Mon 2014-11-03 14:32:49 +0100
message:
  Implement second variant of BSpline kernel function.
  
  Thanks to Sebastian Borrmann.
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-10-24 08:50:30 +0000
+++ pkg/common/SPHEngine.cpp	2014-11-03 13:32:49 +0000
@@ -54,9 +54,9 @@
 
 Real smoothkernelLucy(const double & r, const double & h) {
   if (r<=h && h>0) {
-    // Lucy Kernel function, [Lucy1977] (27) 
+    // Lucy Kernel function, [Lucy1977] (27)
     const Real r_h = r / h;
-    return 105./(16.*M_PI*h*h*h) * (1 + 3 * r_h) * pow((1.0 - r_h), 3);
+    return 105./(16.*M_PI*h*h*h) * (1. + 3. * r_h) * std::pow((1. - r_h), 3);
   } else {
     return 0;
   }
@@ -65,7 +65,7 @@
 Real smoothkernelLucyGrad(const double & r, const double & h) {
   if (r<=h && h>0) {
     // 1st derivative of Lucy Kernel function, [Lucy1977] (27)
-    return 105./(16.*M_PI*h*h*h) * (-12 * r) * pow((h - r), 2) / ( h * h * h * h );
+    return 105./(16.*M_PI*h*h*h) * (-12. * r) * std::pow((h - r), 2) / ( h * h * h * h );
   } else {
     return 0;
   }
@@ -74,52 +74,99 @@
 Real smoothkernelLucyLapl(const double & r, const double & h) {
   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) {
+    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 smoothkernelBSpline1(const double & r, const double & h) {
+  // BSpline Kernel function, [Monaghan1985] (21)
+  if (r<=2.0*h && h>0) {
+    const Real coefA = 3. / (2. * M_PI * h * h * h);
+    const Real r_h = r / h;
+    if (r<=h) {
+      return coefA * (2. /3. - r_h * r_h + 1. / 2. * r_h * r_h * r_h);
+    } else {
+      return coefA / 6. * std::pow((2. - r_h), 3);
+    }
+  } else {
+    return 0;
+  }
+}
+
+Real smoothkernelBSpline1Grad(const double & r, const double & h) {
+  // 1st derivative of BSpline Kernel function, [Monaghan1985] (21)
+  if (r<=2.*h && h>0) {
+    const Real coefA = 3. / (2. * M_PI * h * h * h);
+    const Real r_h = r / h;
+    if (r<=h) {
+      return coefA * (-2. * r_h ) * (1. - 3. / 2. * r_h);
+    } else {
+      return coefA * (-1. / 2.) * std::pow((2. - r_h), 2);
+    }
+  } else {
+    return 0;
+  }
+}
+
+Real smoothkernelBSpline1Lapl(const double & r, const double & h) {
+  // 2nd derivative of BSpline Kernel function, [Monaghan1985] (21)
+  if (r<=2.0*h && h>0) {
+    const Real coefA = 3. / (2. * M_PI * h * h * h);
+    const Real r_h = r / h;
+    if (r<=h) {
+      return coefA * (-2. + 3. * r_h);
+    } else {
+      return coefA * (2. - r_h);
+    }
+  } else {
+    return 0;
+  }
+}
+
+//=========================================================================
+
+Real smoothkernelBSpline2(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);
+    const Real coefA = 3. / (4. * M_PI * h * h * h);
+    const Real r_h = r / h;
+    if (r<=h) {
+      return coefA * (10. /3. - 7. * r_h * r_h + 4 * r_h * r_h * r_h);
+    } else {
+      return coefA * std::pow((2. - r_h), 2)*((5. - 4. * r_h) / 3.);
+    }
+  } else {
+    return 0;
+  }
+}
+
+Real smoothkernelBSpline2Grad(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. / (4. * 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 smoothkernelBSpline2Lapl(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. / (4. * 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;
@@ -145,13 +192,23 @@
     } else {
       KERNELFUNCDESCR
     }
-  } else if (a==BSpline) {
-    if (typeF==Norm) {
-      return smoothkernelBSpline;
-    } else if (typeF==Grad) {
-      return smoothkernelBSplineGrad;
-    } else if (typeF==Lapl) {
-      return smoothkernelBSplineLapl;
+  } else if (a==BSpline1) {
+    if (typeF==Norm) {
+      return smoothkernelBSpline1;
+    } else if (typeF==Grad) {
+      return smoothkernelBSpline1Grad;
+    } else if (typeF==Lapl) {
+      return smoothkernelBSpline1Lapl;
+    } else {
+      KERNELFUNCDESCR
+    }
+  } else if (a==BSpline2) {
+    if (typeF==Norm) {
+      return smoothkernelBSpline2;
+    } else if (typeF==Grad) {
+      return smoothkernelBSpline2Grad;
+    } else if (typeF==Lapl) {
+      return smoothkernelBSpline2Lapl;
     } else {
       KERNELFUNCDESCR
     }

=== modified file 'pkg/common/SPHEngine.hpp'
--- pkg/common/SPHEngine.hpp	2014-10-23 13:46:20 +0000
+++ pkg/common/SPHEngine.hpp	2014-11-03 13:32:49 +0000
@@ -6,8 +6,8 @@
 
 typedef Real (* KernelFunction)(const double & r, const double & h);
 
-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 KernFunctions {Lucy=1,BSpline1=2,BSpline2=3};
+#define KERNELFUNCDESCR throw runtime_error("Type of kernel function undefined! The following kernel functions are available: Lucy=1 ([Lucy1977]_ (27)), BSpline1=2 ([Monaghan1985]_ (21)), BSpline2=3 ([Monaghan1985]_ (22)).");
 
 enum typeKernFunctions {Norm, Grad, Lapl};
 class SPHEngine: public PartialEngine{
@@ -19,16 +19,19 @@
     ((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, BSpline=2."))
+    ((int,KernFunctionDensity, Lucy,, "Kernel function for density calculation (by default - Lucy). The following kernel functions are available: Lucy=1 ([Lucy1977]_ (27)), BSpline1=2 ([Monaghan1985]_ (21)), BSpline2=3 ([Monaghan1985]_ (22))."))
   );
 };
 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);
+Real smoothkernelBSpline1(const double & r, const double & h);         
+Real smoothkernelBSpline1Grad(const double & r, const double & h);     
+Real smoothkernelBSpline1Lapl(const double & r, const double & h);
+Real smoothkernelBSpline2(const double & r, const double & h);         
+Real smoothkernelBSpline2Grad(const double & r, const double & h);     
+Real smoothkernelBSpline2Lapl(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);