-
Notifications
You must be signed in to change notification settings - Fork 16
/
KERNELCUBIC.cpp
67 lines (53 loc) · 1.2 KB
/
KERNELCUBIC.cpp
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
/*
* KERNELCUBIC.cpp
*
* Created on: Jul 8, 2011
* Author: tony
*/
#include "KERNELCUBIC.h"
#include "math.h"
namespace std {
KERNEL_CUBIC::KERNEL_CUBIC(STORAGE *store):KERNEL_BASE(store) {
a1 = 0.25/atan(1);
a2 = a1/(store->h*store->h*store->h);
aa = a1/(store->h*store->h*store->h*store->h);
a24 = 0.25*a2;
c1 = -3*aa;
d1 = 9*aa/4;
c2 = -3*aa/4;
deltap=1/1.5;
Wdeltap=a2*(1.-1.5*deltap*deltap+0.75*deltap*deltap*deltap);
od_Wdeltap=1./Wdeltap;
store->adh = a2;
}
KERNEL_CUBIC::~KERNEL_CUBIC() {
// TODO Auto-generated destructor stub
}
void KERNEL_CUBIC::kernel(double drx,double dry,double drz,int i,int j,int j1,int j2,double drr2)
{
store->index_tensile= 1;
//c Cubic 3D Kernel
double rad = sqrt(drr2);
double qq = rad/store->h;
if(drr2 > store->h*store->h)
{
double tmq = 2.-qq;
Wab = a24*tmq*tmq*tmq;
double fac = c2*tmq*tmq /rad;
store->frx = fac * drx;
if (store->dim == 3)
store->fry = fac * dry;
store->frz = fac * drz;
}
else
{
Wab = a2*( 1. - 1.5*qq*qq + 0.75*qq*qq*qq );
double fac = (c1*qq + d1*qq*qq) /rad;
store->frx = fac * drx;
if (store->dim == 3)
store->fry = fac * dry;
store->frz = fac * drz;
}
return;
}
} /* namespace std */