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
|
#include "fem.h"
void StartUp3d(Mesh *mesh){
#if p_N==1
#include "data3dN01.h"
#elif p_N==2
#include "data3dN02.h"
#elif p_N==3
#include "data3dN03.h"
#elif p_N==4
#include "data3dN04.h"
#elif p_N==5
#include "data3dN05.h"
#elif p_N==6
#include "data3dN06.h"
#elif p_N==7
#include "data3dN07.h"
#elif p_N==8
#include "data3dN08.h"
#elif p_N==9
#include "data3dN09.h"
#endif
int n, m, k;
/* load r, s, t */
mesh->r = BuildVector(p_Np);
mesh->s = BuildVector(p_Np);
mesh->t = BuildVector(p_Np);
for(n=0;n<p_Np;++n){
mesh->r[n] = p_r[n];
mesh->s[n] = p_s[n];
mesh->t[n] = p_t[n];
}
/* load Dr, Ds, Dt */
mesh->Dr = BuildMatrix(p_Np, p_Np);
mesh->Ds = BuildMatrix(p_Np, p_Np);
mesh->Dt = BuildMatrix(p_Np, p_Np);
for(n=0;n<p_Np;++n){
for(m=0;m<p_Np;++m){
mesh->Dr[n][m] = p_Dr[n][m];
mesh->Ds[n][m] = p_Ds[n][m];
mesh->Dt[n][m] = p_Dt[n][m];
}
}
/* load LIFT */
mesh->LIFT = BuildMatrix(p_Np, p_Nfp*p_Nfaces);
for(n=0;n<p_Np;++n){
for(m=0;m<p_Nfp*p_Nfaces;++m){
mesh->LIFT[n][m] = p_LIFT[n][m];
}
}
mesh->Fmask = BuildIntMatrix(p_Nfaces, p_Nfp);
for(n=0;n<p_Nfaces;++n){
for(m=0;m<p_Nfp;++m){
mesh->Fmask[n][m] = p_Fmask[n][m];
}
}
/* low storage RK coefficients */
mesh->rk4a = BuildVector(5);
mesh->rk4a[0] = 0.0;
mesh->rk4a[1] = -567301805773.0 / 1357537059087.0;
mesh->rk4a[2] = -2404267990393.0 / 2016746695238.0;
mesh->rk4a[3] = -3550918686646.0 / 2091501179385.0;
mesh->rk4a[4] = -1275806237668.0 / 842570457699.0;
mesh->rk4b = BuildVector(5);
mesh->rk4b[0] = 1432997174477.0 / 9575080441755.0;
mesh->rk4b[1] = 5161836677717.0 / 13612068292357.0;
mesh->rk4b[2] = 1720146321549.0 / 2090206949498.0;
mesh->rk4b[3] = 3134564353537.0 / 4481467310338.0;
mesh->rk4b[4] = 2277821191437.0 / 14882151754819.0;
mesh->rk4c = BuildVector(6);
mesh->rk4c[0] = 0.0;
mesh->rk4c[1] = 1432997174477.0 / 9575080441755.0;
mesh->rk4c[2] = 2526269341429.0 / 6820363962896.0;
mesh->rk4c[3] = 2006345519317.0 / 3224310063776.0;
mesh->rk4c[4] = 2802321613138.0 / 2924317926251.0;
mesh->rk4c[5] = 1.0;
/* build coordinates */
mesh->x = BuildMatrix(mesh->K, p_Np);
mesh->y = BuildMatrix(mesh->K, p_Np);
mesh->z = BuildMatrix(mesh->K, p_Np);
for(k=0;k<mesh->K;++k){
for(n=0;n<p_Np;++n){
double r = mesh->r[n];
double s = mesh->s[n];
double t = mesh->t[n];
mesh->x[k][n] = 0.5*( -mesh->GX[k][0]*(r+s+t+1) +
mesh->GX[k][1]*(1.+r) +
mesh->GX[k][2]*(1.+s) +
mesh->GX[k][3]*(1.+t)
);
mesh->y[k][n] = 0.5*( -mesh->GY[k][0]*(r+s+t+1) +
mesh->GY[k][1]*(1.+r) +
mesh->GY[k][2]*(1.+s) +
mesh->GY[k][3]*(1.+t)
);
mesh->z[k][n] = 0.5*( -mesh->GZ[k][0]*(r+s+t+1) +
mesh->GZ[k][1]*(1.+r) +
mesh->GZ[k][2]*(1.+s) +
mesh->GZ[k][3]*(1.+t)
);
}
}
/* build node-node connectivity maps */
void BuildMaps3d(Mesh *mesh);
BuildMaps3d(mesh);
}
|
