Merge pull request #308 from aperijake/add_cardinal_b_spline_weightin…

…g_function

Add the cardinal cubic b-spline as a weighting function option

cmake/Compadre_Version.txt

@@ -1 +1 @@
-1.5.8
+1.5.9

pycompadre/examples/pycompadre_kernels.py

@@ -61,16 +61,17 @@ def approximate(wt, h, p, n):
 t_power = plt.text(0.02, 0.05, '$$\\left(1-\\left(\\frac{|r|}{h}\\right)^N\\right)^P$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
 t_gaussian = plt.text(0.02, 0.05, '$$\\frac{1}{\\frac{h}{P} \\sqrt(2 \\pi)} e^{\\frac{-1}{2}r^2/{\\left(\\frac{h}{P}\\right)^2}}$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
 t_cubic = plt.text(0.02, 0.05, '$$((1-\\frac{r}{h})+\\frac{r}{h} (1-\\frac{r}{h}) (1-2 \\frac{r}{h}))$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
+t_cubic_b = plt.text(0.02, 0.05, '$$(1 + 6x^2(-1 + x)), \\frac{r}{h} \\le 0.5\\\\;\\; 2(1 + x(-3 + 3x - x^2)), 0.5< \\frac{r}{h} < 1.0$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
 t_cosine = plt.text(0.02, 0.05, '$$\\cos(\\frac{\\pi r}{2h})$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
 t_sigmoid = plt.text(0.02, 0.05, '$$\\frac{1}{e^{Pr} + e^{-Pr} + N}$$', fontsize=fontsize, transform=plt.gcf().transFigure,usetex=True)
 
-t_labels = {'Power':t_power, 'Gaussian':t_gaussian, 'Cubic Spl.':t_cubic, 'Cosine':t_cosine, 'Sigmoid':t_sigmoid}
+t_labels = {'Power':t_power, 'Gaussian':t_gaussian, 'Cubic Spl.':t_cubic, 'Cubic B-Spl.':t_cubic_b, 'Cosine':t_cosine, 'Sigmoid':t_sigmoid}
 for item in t_labels.keys():
     if item!='Power':
         t_labels[item].update({'visible':False})
 
 #radios
-rad_weighting_type = RadioButtons(ax_weighting_type, ('Power', 'Cubic Spl.', 'Cosine', 'Gaussian', 'Sigmoid'), active=0)
+rad_weighting_type = RadioButtons(ax_weighting_type, ('Power', 'Cubic Spl.', 'Cubic B-Spl.', 'Cosine', 'Gaussian', 'Sigmoid'), active=0)
 
 def update(val):
     global wt
@@ -83,7 +84,7 @@ def update(val):
 sl_n.on_changed(update)
 
 def weighting_type_update(label):
-    weighting_type_dict = {'Power': pycompadre.WeightingFunctionType.Power, 'Cubic Spl.': pycompadre.WeightingFunctionType.CubicSpline, 'Cosine': pycompadre.WeightingFunctionType.Cosine, 'Gaussian': pycompadre.WeightingFunctionType.Gaussian, 'Sigmoid': pycompadre.WeightingFunctionType.Sigmoid}
+    weighting_type_dict = {'Power': pycompadre.WeightingFunctionType.Power, 'Cubic Spl.': pycompadre.WeightingFunctionType.CubicSpline, 'Cosine': pycompadre.WeightingFunctionType.Cosine, 'Gaussian': pycompadre.WeightingFunctionType.Gaussian, 'Sigmoid': pycompadre.WeightingFunctionType.Sigmoid, 'Cubic B-Spl.': pycompadre.WeightingFunctionType.CardinalCubicBSpline}
     global wt
     wt = weighting_type_dict[label]
     for item in t_labels.keys():
@@ -95,6 +96,9 @@ def weighting_type_update(label):
     elif (label=='Cubic Spl.'):
         ax_p.set_visible(False)
         ax_n.set_visible(False)
+    elif (label=='Cubic B-Spl.'):
+        ax_p.set_visible(False)
+        ax_n.set_visible(False)
     elif (label=='Cosine'):
         ax_p.set_visible(False)
         ax_n.set_visible(False)

pycompadre/pycompadre.cpp

@@ -784,6 +784,7 @@ Implementation details at:
     .value("Power", WeightingFunctionType::Power)
     .value("Gaussian", WeightingFunctionType::Gaussian)
     .value("CubicSpline", WeightingFunctionType::CubicSpline)
+    .value("CardinalCubicBSpline", WeightingFunctionType::CardinalCubicBSpline)
     .value("Cosine", WeightingFunctionType::Cosine)
     .value("Sigmoid", WeightingFunctionType::Sigmoid)
     .export_values();

src/Compadre_GMLS.hpp

@@ -556,6 +556,14 @@ friend class Evaluator;
             // invariant to p and n
             double x = std::abs(r/h);
             return ((1-x)+x*(1-x)*(1-2*x)) * double(x<=1.0);
+        } else if (weighting_type == WeightingFunctionType::CardinalCubicBSpline) {
+            // compactly supported on [0,h]
+            // invariant to p and n
+            // Calculate the value using a cardinal cubic b-spline kernel (often just called cubic b spline)
+            double x = std::abs(r/h);
+            if (x < 0.5) return 1.0 + 6.0 * x * x * (-1.0 + x);
+            if (x < 1.0) return 2.0 * (1.0 + x * (-3.0 + 3.0 * x - 1.0 * x * x));
+            return 0.0;
         } else if (weighting_type == WeightingFunctionType::Cosine) {
             // compactly supported on [0,h]
             double pi = 3.14159265358979323846;
@@ -1161,6 +1169,8 @@ friend class Evaluator;
             _weighting_type = WeightingFunctionType::Gaussian;
         } else if (wt_to_lower == "cubicspline") {
             _weighting_type = WeightingFunctionType::CubicSpline;
+        } else if (wt_to_lower == "cardinalcubicbspline") {
+            _weighting_type = WeightingFunctionType::CardinalCubicBSpline;
         } else if (wt_to_lower == "cosine") {
             _weighting_type = WeightingFunctionType::Cosine;
         } else if (wt_to_lower == "sigmoid") {
@@ -1188,6 +1198,8 @@ friend class Evaluator;
             _curvature_weighting_type = WeightingFunctionType::Gaussian;
         } else if (wt_to_lower == "cubicspline") {
             _curvature_weighting_type = WeightingFunctionType::CubicSpline;
+        } else if (wt_to_lower == "cardinalcubicbspline") {
+            _curvature_weighting_type = WeightingFunctionType::CardinalCubicBSpline;
         } else if (wt_to_lower == "cosine") {
             _curvature_weighting_type = WeightingFunctionType::Cosine;
         } else if (wt_to_lower == "sigmoid") {

src/Compadre_Operators.hpp

@@ -235,6 +235,7 @@ namespace Compadre {
         Power,
         Gaussian,
         CubicSpline,
+        CardinalCubicBSpline,
         Cosine,
         Sigmoid
     };