1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101

use crate::kernel::Kernel;
use crate::math::Real;
use na::RealField;

/// The cubic spline smoothing kernel.
///
/// See https://pysph.readthedocs.io/en/latest/reference/kernels.html
#[derive(Copy, Clone, Debug)]
pub struct CubicSplineKernel;

impl Kernel for CubicSplineKernel {
    fn scalar_apply(r: Real, h: Real) -> Real {
        assert!(r >= na::zero::<Real>());

        #[cfg(feature = "dim2")]
        let normalizer = na::convert::<_, Real>(40.0 / 7.0) / (Real::pi() * h * h);
        #[cfg(feature = "dim3")]
        let normalizer = na::convert::<_, Real>(8.0) / (Real::pi() * h * h * h);

        let _2: Real = na::convert::<_, Real>(2.0);
        let q = r / h;

        let rhs = if q <= na::convert::<_, Real>(0.5) {
            let q2 = q * q;
            na::one::<Real>() + (q2 * q - q2) * na::convert::<_, Real>(6.0)
        } else if q <= na::one::<Real>() {
            (na::one::<Real>() - q).powi(3) * _2
        } else {
            na::zero::<Real>()
        };

        normalizer * rhs

        /*
        let q = r / h;
        #[cfg(feature = "dim2")]
            let normalizer = na::convert::<_, Real>(10.0 / 7.0) / (Real::pi() * h * h);
        #[cfg(feature = "dim3")]
            let normalizer = na::one::<Real>() / (Real::pi() * h * h * h);

        let _2: Real = na::convert::<_, Real>(2.0);
        let _3: Real = na::convert::<_, Real>(3.0);
        let rhs = if q <= na::one::<Real>() {
            na::one::<Real>() - _3 / _2 * q * q * (na::one::<Real>() - q / _2)
        } else if q <= _2 {
            (_2 - q).powi(3) / na::convert::<_, Real>(4.0)
        } else {
            na::zero::<Real>()
        };

        normalizer * rhs
        */
    }

    fn scalar_apply_diff(r: Real, h: Real) -> Real {
        assert!(r >= na::zero::<Real>());

        #[cfg(feature = "dim2")]
        let normalizer = na::convert::<_, Real>(40.0 / 7.0) / (Real::pi() * h * h);
        #[cfg(feature = "dim3")]
        let normalizer = na::convert::<_, Real>(8.0) / (Real::pi() * h * h * h);

        let _1: Real = na::convert::<_, Real>(1.0);
        let _2: Real = na::convert::<_, Real>(2.0);
        let _3: Real = na::convert::<_, Real>(3.0);
        let _eps: Real = na::convert::<_, Real>(1.0e-5);
        let q = r / h;

        let rhs = if q > _1 || q <= _eps {
            na::zero::<Real>()
        } else if q <= na::convert::<_, Real>(0.5) {
            (q * _3 - _2) * q * na::convert::<_, Real>(6.0)
        } else {
            // 0.5 < q <= 1.0
            let one_q = _1 - q;
            -one_q * one_q * na::convert::<_, Real>(6.0)
        };

        normalizer * rhs / h

        /*
        let q = r / h;
        #[cfg(feature = "dim2")]
            let normalizer = na::convert::<_, Real>(10.0 / 7.0) / (Real::pi() * h * h);
        #[cfg(feature = "dim3")]
            let normalizer = na::one::<Real>() / (Real::pi() * h * h * h);

        let _2: Real = na::convert::<_, Real>(2.0);
        let _3: Real = na::convert::<_, Real>(3.0);
        let rhs = if q <= na::one::<Real>() {
            -_3 * q * (na::one::<Real>() - q * na::convert::<_, Real>(3.0 / 4.0))
        } else if q <= _2 {
            -(_2 - q).powi(2) * na::convert::<_, Real>(3.0 / 4.0)
        } else {
            na::zero::<Real>()
        };

        normalizer * rhs
        */
    }
}