sparkl2d_core/dynamics/solver/
kernel.rs

1use crate::math::{Point, Real, Vector, DIM};
2use na::vector;
3#[cfg(not(feature = "std"))]
4use na::ComplexField;
5
6pub struct QuadraticKernel;
7
8impl QuadraticKernel {
9    #[inline(always)]
10    pub fn inv_d(cell_width: Real) -> Real {
11        4.0 / (cell_width * cell_width)
12    }
13
14    #[inline(always)]
15    pub fn eval_all(x: Real) -> [Real; 3] {
16        [
17            0.5 * (1.5 - x).powi(2),
18            0.75 - (x - 1.0).powi(2),
19            0.5 * (x - 0.5).powi(2),
20        ]
21    }
22
23    // #[inline(always)]
24    // pub fn eval_all_simd(x: SimdReal) -> [SimdReal; 3] {
25    //     let a = SimdReal::splat(1.5) - x;
26    //     let b = x - SimdReal::splat(1.0);
27    //     let c = x - SimdReal::splat(0.5);
28    //
29    //     [
30    //         SimdReal::splat(0.5) * a * a,
31    //         SimdReal::splat(0.75) * b * b,
32    //         SimdReal::splat(0.5) * c * c,
33    //     ]
34    // }
35
36    #[inline(always)]
37    pub fn eval(x: Real) -> Real {
38        let x_abs = x.abs();
39
40        if x_abs < 0.5 {
41            3.0 / 4.0 - x_abs.powi(2)
42        } else if x_abs < 3.0 / 2.0 {
43            0.5 * (3.0 / 2.0 - x_abs).powi(2)
44        } else {
45            0.0
46        }
47    }
48
49    #[inline(always)]
50    pub fn eval_derivative(x: Real) -> Real {
51        let x_abs = x.abs();
52
53        if x_abs < 0.5 {
54            -2.0 * x.signum() * x_abs
55        } else if x_abs < 3.0 / 2.0 {
56            -x.signum() * (3.0 / 2.0 - x_abs)
57        } else {
58            0.0
59        }
60    }
61
62    #[inline(always)]
63    pub fn precompute_weights(
64        ref_elt_pos_minus_particle_pos: Vector<Real>,
65        h: Real,
66    ) -> [[Real; 3]; DIM] {
67        [
68            Self::eval_all(-ref_elt_pos_minus_particle_pos.x / h),
69            Self::eval_all(-ref_elt_pos_minus_particle_pos.y / h),
70            #[cfg(feature = "dim3")]
71            Self::eval_all(-ref_elt_pos_minus_particle_pos.z / h),
72        ]
73    }
74
75    // #[inline(always)]
76    // pub fn precompute_weights_simd(
77    //     ref_elt_pos_minus_particle_pos: Vector<SimdReal>,
78    //     h: SimdReal,
79    // ) -> [[SimdReal; 3]; DIM] {
80    //     [
81    //         Self::eval_all_simd(-ref_elt_pos_minus_particle_pos.x / h),
82    //         Self::eval_all_simd(-ref_elt_pos_minus_particle_pos.y / h),
83    //         #[cfg(feature = "dim3")]
84    //         Self::eval_all_simd(-ref_elt_pos_minus_particle_pos.z / h),
85    //     ]
86    // }
87
88    #[inline(always)]
89    pub fn stencil_with_dir(elt_pos_minus_particle_pos: Vector<Real>, h: Real) -> Real {
90        let dpt = -elt_pos_minus_particle_pos / h;
91        #[cfg(feature = "dim2")]
92        return Self::eval(dpt.x) * Self::eval(dpt.y);
93        #[cfg(feature = "dim3")]
94        return Self::eval(dpt.x) * Self::eval(dpt.y) * Self::eval(dpt.z);
95    }
96
97    #[inline(always)]
98    pub fn stencil(elt_pos: Point<Real>, particle_pos: Point<Real>, h: Real) -> Real {
99        Self::stencil_with_dir(elt_pos - particle_pos, h)
100    }
101
102    #[inline(always)]
103    pub fn stencil_gradient_with_dir(
104        elt_pos_minus_particle_pos: Vector<Real>,
105        h: Real,
106    ) -> Vector<Real> {
107        let inv_h = 1.0 / h;
108        let dpt = -elt_pos_minus_particle_pos * inv_h;
109        let val_x = Self::eval(dpt.x);
110        let val_y = Self::eval(dpt.y);
111
112        #[cfg(feature = "dim3")]
113        let val_z = Self::eval(dpt.z);
114
115        #[cfg(feature = "dim2")]
116        return vector![
117            inv_h * Self::eval_derivative(dpt.x) * val_y,
118            inv_h * val_x * Self::eval_derivative(dpt.y)
119        ];
120        #[cfg(feature = "dim3")]
121        return vector![
122            inv_h * Self::eval_derivative(dpt.x) * val_y * val_z,
123            inv_h * val_x * Self::eval_derivative(dpt.y) * val_z,
124            inv_h * val_x * val_y * Self::eval_derivative(dpt.z)
125        ];
126    }
127
128    #[inline(always)]
129    pub fn stencil_gradient(
130        elt_pos: Point<Real>,
131        particle_pos: Point<Real>,
132        h: Real,
133    ) -> Vector<Real> {
134        Self::stencil_gradient_with_dir(elt_pos - particle_pos, h)
135    }
136}
sparkl2d_core/dynamics/solver/kernel.rs

sparkl2d_core/dynamics/solver/
kernel.rs
