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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
/* -*- c++ -*- ----------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
#ifndef SMD_KERNEL_FUNCTIONS_H_
#define SMD_KERNEL_FUNCTIONS_H_
namespace SMD_Kernels {
static inline double Kernel_Wendland_Quintic_NotNormalized(const double r, const double h) {
if (r < h) {
double q = 2.0 * r / h;
return pow(1.0 - 0.5 * q, 4) * (2.0 * q + 1.0);
} else {
return 0.0;
}
}
static inline double Kernel_Cubic_Spline(const double r, const double h) {
double q = 2.0 * r / h;
if (q > 2.0) {
return 0.0;
} else if ((q <= 2.0) && (q > 1.0)) {
return pow(2.0 - q, 3.0) / 6.0;
} else if ((q >= 0.0) && (q <= 1.0)) {
return 2. / 3. - q * q + 0.5 * q * q * q;
} else {
return 0.0;
}
}
static inline double Kernel_Barbara(const double r, const double h) {
double arg = (1.570796327 * (r + h)) / h;
double hsq = h * h;
//wf = (1.680351548 * (cos(arg) + 1.)) / hsq;
return -2.639490040 * sin(arg) / (hsq * h);
}
static inline void spiky_kernel_and_derivative(const double h, const double r, const int dimension, double &wf, double &wfd) {
/*
* Spiky kernel
*/
if (r > h) {
printf("r=%f > h=%f in Spiky kernel\n", r, h);
wf = wfd = 0.0;
return;
}
double hr = h - r; // [m]
if (dimension == 2) {
double n = 0.3141592654e0 * h * h * h * h * h; // [m^5]
wfd = -3.0e0 * hr * hr / n; // [m*m/m^5] = [1/m^3] ==> correct for dW/dr in 2D
wf = -0.333333333333e0 * hr * wfd; // [m/m^3] ==> [1/m^2] correct for W in 2D
} else {
wfd = -14.0323944878e0 * hr * hr / (h * h * h * h * h * h); // [1/m^4] ==> correct for dW/dr in 3D
wf = -0.333333333333e0 * hr * wfd; // [m/m^4] ==> [1/m^3] correct for W in 3D
}
// alternative formulation
// double hr = h - r;
//
// /*
// * Spiky kernel
// */
//
// if (domain->dimension == 2) {
// double h5 = h * h * h * h * h;
// wf = 3.183098861e0 * hr * hr * hr / h5;
// wfd = -9.549296583 * hr * hr / h5;
//
// } else {
// double h6 = h * h * h * h * h * h;
// wf = 4.774648292 * hr * hr * hr / h6;
// wfd = -14.32394487 * hr * hr / h6;
// }
// }
}
static inline void barbara_kernel_and_derivative(const double h, const double r, const int dimension, double &wf, double &wfd) {
/*
* Barbara kernel
*/
double arg = (1.570796327 * (r + h)) / h;
double hsq = h * h;
if (r > h) {
printf("r = %f > h = %f in barbara kernel function\n", r, h);
exit(1);
//wf = wfd = 0.0;
//return;
}
if (dimension == 2) {
wf = (1.680351548 * (cos(arg) + 1.)) / hsq;
wfd = -2.639490040 * sin(arg) / (hsq * h);
} else {
wf = 2.051578323 * (cos(arg) + 1.) / (hsq * h);
wfd = -3.222611694 * sin(arg) / (hsq * hsq);
}
}
/*
* compute a normalized smoothing kernel only
*/
static inline void Poly6Kernel(const double hsq, const double h, const double rsq, const int dimension, double &wf) {
double tmp = hsq - rsq;
if (dimension == 2) {
wf = tmp * tmp * tmp / (0.7853981635e0 * hsq * hsq * hsq * hsq);
} else {
wf = tmp * tmp * tmp / (0.6382918409e0 * hsq * hsq * hsq * hsq * h);
}
}
/*
* M4 Prime Kernel
*/
static inline void M4PrimeKernel(const double s, double &wf) {
if (s < 1.0) {
//wf = 1.0 - 2.5 * s * s + (3./2.) * s * s * s;
wf = 1.0 - s * s *(2.5 -1.5 *s);
} else if (s < 2.0) {
//wf = 0.5 * (1.0 - s) * ((2.0 - s) * (2.0 - s));
wf = 2.0 + (-4.0 + (2.5 - 0.5 * s)*s)*s;
} else {
wf = 0.0;
}
}
}
#endif /* SMD_KERNEL_FUNCTIONS_H_ */