1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2016 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Hammad Mazhar, Radu Serban
13 // =============================================================================
14 //
15 // Utilities for fluid SPH kernels
16 //
17 // =============================================================================
18
19 #ifndef CH_FLUID_KERNELS_H
20 #define CH_FLUID_KERNELS_H
21
22 #include "chrono/multicore_math/ChMulticoreMath.h"
23
24 namespace chrono {
25
26 #define F_PI 3.141592653589793238462643383279
27 #define INVPI (1 / F_PI)
28
29 #define KERNEL poly6
30 #define GRAD_KERNEL unormalized_grad_spiky
31 #define GRAD2_KERNEL grad2_viscosity
32
33 #define H2 h* h
34 #define H3 h* h* h
35 #define H6 H3* H3
36 #define H9 H3* H3* H3
37
38 ///
39 #define CPOLY6 315.0 / (64.0 * F_PI * H9)
40 #define KPOLY6 CPOLY6 * Pow((H2 - dist * dist), 3)
41
42 #define CGPOLY6 -945.0 / (32.0 * F_PI * H9)
43 #define KGPOLY6 CGPOLY6 * Pow((H2 - dist * dist), 2)
44
45 #define CLPOLY6 945.0 / (32.0 * F_PI * H9)
46 #define KLPOLY6 CLPOLY6 * (H2 - dist * dist) * (7 * dist * dist - 3 * H2)
47
48 ///
49 #define CGSPIKY -45.0 / (F_PI * H6)
50 #define KGSPIKY CGSPIKY * Pow(h - dist, 2)
51
52 #define CLVISC 45.0 / (F_PI * H6)
53 #define KLVISC CLVISC * (h - dist)
54
N(const real & dist,const real & h)55 inline real N(const real& dist, const real& h) {
56 real x = Abs(dist) / h;
57 if (Abs(x) < real(1.0)) {
58 return real(0.5) * Cube(Abs(x)) - Sqr(x) + 2.0 / 3.0;
59 } else if (Abs(x) < real(2.0)) {
60 return -1.0 / 6.0 * Cube(Abs(x)) + Sqr(x) - real(2.0) * Abs(x) + 4.0 / 3.0;
61 }
62 return real(0.0);
63 }
64
65 // Cubic spline kernel
66 // d is positive. h is the sph particle radius (i.e. h in the document) d is the distance of 2 particles
cubic_spline(const real & dist,const real & h)67 inline real cubic_spline(const real& dist, const real& h) {
68 real q = Abs(dist) / h;
69 if (q < 1) {
70 return (0.25f / (F_PI * h * h * h) * (Pow(2 - q, 3) - 4 * Pow(1 - q, 3)));
71 }
72 if (q < 2) {
73 return (0.25f / (F_PI * h * h * h) * Pow(2 - q, 3));
74 }
75 return 0;
76 }
77 // d is positive. r is the sph particles
grad_cubic_spline(const real3 & dist,const real d,const real & h)78 inline real3 grad_cubic_spline(const real3& dist, const real d, const real& h) {
79 real q = d / h;
80
81 if (q < 1) {
82 return (3 * q - 4) * .75 * INVPI * Pow(h, -5) * dist;
83 }
84 if (q < 2) {
85 return (-q + 4.0 - 4.0 / q) * .75 * INVPI * Pow(h, -5) * dist;
86 }
87 return real3(0);
88 }
poly6(const real & dist,const real & h)89 inline real poly6(const real& dist, const real& h) {
90 return /* (dist <= h)* */ KPOLY6;
91 }
92
grad_poly6(const real3 & xij,const real d,const real & h)93 inline real3 grad_poly6(const real3& xij, const real d, const real& h) {
94 return (d <= h) * -945.0 / (32.0 * F_PI * Pow(h, 9)) * Pow((h * h - d * d), 2) * xij;
95 }
96
spiky(const real & dist,const real & h)97 inline real spiky(const real& dist, const real& h) {
98 return (dist <= h) * 15.0 / (F_PI * Pow(h, 6)) * Pow(h - dist, 3);
99 }
grad_spiky(const real3 & xij,const real dist,const real & h)100 inline real3 grad_spiky(const real3& xij, const real dist, const real& h) {
101 return (dist <= h) * KGSPIKY * xij;
102 }
103
unormalized_spiky(const real & dist,const real & h)104 inline real unormalized_spiky(const real& dist, const real& h) {
105 const real k = 15.0 / (F_PI * h * h * h);
106 return k * Sqr(1.0 - dist / h);
107 }
108
unormalized_grad_spiky(const real3 & xij,const real d,const real & h)109 inline real3 unormalized_grad_spiky(const real3& xij, const real d, const real& h) {
110 const real k = 15.0 / (F_PI * h * h * h);
111 return -k * (1.0 - d / h) / h * xij / d;
112 }
113
viscosity(const real3 & xij,const real d,const real & h)114 inline real3 viscosity(const real3& xij, const real d, const real& h) {
115 return (d <= h) * 15.0 / (2 * F_PI * Pow(h, 3)) *
116 (-(d * d * d) / (2 * h * h * h) + (d * d) / (h * h) + (h) / (2 * d) - 1) * xij;
117 }
grad2_viscosity(const real3 & xij,const real d,const real & h)118 inline real3 grad2_viscosity(const real3& xij, const real d, const real& h) {
119 return real3((d <= h) * 45.0 / (F_PI * Pow(h, 6)) * (h - d));
120 }
121
122 // kernel from constraint fluid approximation paper/code
kernel(const real & dist,const real & h)123 inline real kernel(const real& dist, const real& h) {
124 if (dist > h) {
125 return 0;
126 }
127
128 return Pow(1 - Pow(dist / h, 2), 3);
129 }
130
131 // laplacian operator for poly6
grad2_poly6(const real & dist,const real & h)132 inline real grad2_poly6(const real& dist, const real& h) {
133 if (dist > h) {
134 return 0;
135 }
136 return 945.0 / (32.0 * F_PI * Pow(h, 9)) * (h * h - dist * dist) * (7 * dist * dist - 3 * h * h);
137 }
138
139 } // end namespace chrono
140
141 #endif
142