/[escript]/trunk/ripley/src/Assemble_Interpolation_3D.c
ViewVC logotype

Contents of /trunk/ripley/src/Assemble_Interpolation_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3542 - (show annotations)
Fri Jul 22 06:57:41 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/plain
File size: 26594 byte(s)
first attempts toward a multi grid solver
1 Ripley_Assemble_Interpolation_3D(Ripley_Grid *grid, Escript in, Escript out)
2 {
3 /* GENERATOR SNIP TOP */
4 if (out_data_type==RIPLEY_ELEMENTS) {
5 const double tmp0_0 = 0.69856850115866751967;
6 const double tmp0_4 = 0.39364916731037084426;
7 const double tmp0_1 = 0.088729833462074168852;
8 const double tmp0_8 = 0.028175416344814577871;
9 const double tmp0_5 = 0.050000000000000000000;
10 const double tmp0_2 = 0.011270166537925831148;
11 const double tmp0_9 = 0.12500000000000000000;
12 const double tmp0_6 = 0.0063508326896291557410;
13 const double tmp0_3 = 0.0014314988413324803339;
14 const double tmp0_7 = 0.22182458365518542213;
15 #pragma omp parallel for private(i,k2,k1,k0)
16 for (k2 =0; k2 < N2; ++k2) {
17 for (k1 =0; k1 < N1; ++k1) {
18 for (k0 =0; k0 < N0; ++k0) {
19 for (i =0; i < NCOMP; ++i) {
20 register const double f_000 = in[INDEX2(i,INDEX3(k0,k1,k2, M0,M1),NCOMP)];
21 register const double f_001 = in[INDEX2(i,INDEX3(k0,k1,k2+1, M0,M1),NCOMP)];
22 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,k1,k2+1, M0,M1),NCOMP)];
23 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,k1,k2, M0,M1),NCOMP)];
24 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,k1+1,k2, M0,M1),NCOMP)];
25 register const double f_011 = in[INDEX2(i,INDEX3(k0,k1+1,k2+1, M0,M1),NCOMP)];
26 register const double f_010 = in[INDEX2(i,INDEX3(k0,k1+1,k2, M0,M1),NCOMP)];
27 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,k1+1,k2+1, M0,M1),NCOMP)];
28 out[INDEX3(i,0,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_000*tmp0_0 + f_111*tmp0_3 + tmp0_1*(f_001 + f_010 + f_100) + tmp0_2*(f_011 + f_101 + f_110);
29 out[INDEX3(i,1,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_000 + f_100) + tmp0_5*(f_001 + f_010 + f_101 + f_110) + tmp0_6*(f_011 + f_111);
30 out[INDEX3(i,2,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_011*tmp0_3 + f_100*tmp0_0 + tmp0_1*(f_000 + f_101 + f_110) + tmp0_2*(f_001 + f_010 + f_111);
31 out[INDEX3(i,3,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_000 + f_010) + tmp0_5*(f_001 + f_011 + f_100 + f_110) + tmp0_6*(f_101 + f_111);
32 out[INDEX3(i,4,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_000 + f_010 + f_100 + f_110) + tmp0_8*(f_001 + f_011 + f_101 + f_111);
33 out[INDEX3(i,5,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_100 + f_110) + tmp0_5*(f_000 + f_010 + f_101 + f_111) + tmp0_6*(f_001 + f_011);
34 out[INDEX3(i,6,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_010*tmp0_0 + f_101*tmp0_3 + tmp0_1*(f_000 + f_011 + f_110) + tmp0_2*(f_001 + f_100 + f_111);
35 out[INDEX3(i,7,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_010 + f_110) + tmp0_5*(f_000 + f_011 + f_100 + f_111) + tmp0_6*(f_001 + f_101);
36 out[INDEX3(i,8,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_001*tmp0_3 + f_110*tmp0_0 + tmp0_1*(f_010 + f_100 + f_111) + tmp0_2*(f_000 + f_011 + f_101);
37 out[INDEX3(i,9,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_000 + f_001) + tmp0_5*(f_010 + f_011 + f_100 + f_101) + tmp0_6*(f_110 + f_111);
38 out[INDEX3(i,10,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_000 + f_001 + f_100 + f_101) + tmp0_8*(f_010 + f_011 + f_110 + f_111);
39 out[INDEX3(i,11,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_100 + f_101) + tmp0_5*(f_000 + f_001 + f_110 + f_111) + tmp0_6*(f_010 + f_011);
40 out[INDEX3(i,12,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_000 + f_001 + f_010 + f_011) + tmp0_8*(f_100 + f_101 + f_110 + f_111);
41 out[INDEX3(i,13,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_9*(f_000 + f_001 + f_010 + f_011 + f_100 + f_101 + f_110 + f_111);
42 out[INDEX3(i,14,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_100 + f_101 + f_110 + f_111) + tmp0_8*(f_000 + f_001 + f_010 + f_011);
43 out[INDEX3(i,15,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_010 + f_011) + tmp0_5*(f_000 + f_001 + f_110 + f_111) + tmp0_6*(f_100 + f_101);
44 out[INDEX3(i,16,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_010 + f_011 + f_110 + f_111) + tmp0_8*(f_000 + f_001 + f_100 + f_101);
45 out[INDEX3(i,17,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_110 + f_111) + tmp0_5*(f_010 + f_011 + f_100 + f_101) + tmp0_6*(f_000 + f_001);
46 out[INDEX3(i,18,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_001*tmp0_0 + f_110*tmp0_3 + tmp0_1*(f_000 + f_011 + f_101) + tmp0_2*(f_010 + f_100 + f_111);
47 out[INDEX3(i,19,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_001 + f_101) + tmp0_5*(f_000 + f_011 + f_100 + f_111) + tmp0_6*(f_010 + f_110);
48 out[INDEX3(i,20,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_010*tmp0_3 + f_101*tmp0_0 + tmp0_1*(f_001 + f_100 + f_111) + tmp0_2*(f_000 + f_011 + f_110);
49 out[INDEX3(i,21,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_001 + f_011) + tmp0_5*(f_000 + f_010 + f_101 + f_111) + tmp0_6*(f_100 + f_110);
50 out[INDEX3(i,22,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_7*(f_001 + f_011 + f_101 + f_111) + tmp0_8*(f_000 + f_010 + f_100 + f_110);
51 out[INDEX3(i,23,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_101 + f_111) + tmp0_5*(f_001 + f_011 + f_100 + f_110) + tmp0_6*(f_000 + f_010);
52 out[INDEX3(i,24,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_011*tmp0_0 + f_100*tmp0_3 + tmp0_1*(f_001 + f_010 + f_111) + tmp0_2*(f_000 + f_101 + f_110);
53 out[INDEX3(i,25,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = tmp0_4*(f_011 + f_111) + tmp0_5*(f_001 + f_010 + f_101 + f_110) + tmp0_6*(f_000 + f_100);
54 out[INDEX3(i,26,INDEX3(k0,k1,k2,N0,N1),NCOMP,27)] = f_000*tmp0_3 + f_111*tmp0_0 + tmp0_1*(f_011 + f_101 + f_110) + tmp0_2*(f_001 + f_010 + f_100);
55 } /* close component loop i */
56 } /* close k0 loop */
57 } /* close k1 loop */
58 } /* close k2 loop */
59 } else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {
60 const double tmp0_0 = 0.12500000000000000000;
61 #pragma omp parallel for private(i,k2,k1,k0)
62 for (k2 =0; k2 < N2; ++k2) {
63 for (k1 =0; k1 < N1; ++k1) {
64 for (k0 =0; k0 < N0; ++k0) {
65 for (i =0; i < NCOMP; ++i) {
66 register const double f_000 = in[INDEX2(i,INDEX3(k0,k1,k2, M0,M1),NCOMP)];
67 register const double f_001 = in[INDEX2(i,INDEX3(k0,k1,k2+1, M0,M1),NCOMP)];
68 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,k1,k2+1, M0,M1),NCOMP)];
69 register const double f_011 = in[INDEX2(i,INDEX3(k0,k1+1,k2+1, M0,M1),NCOMP)];
70 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,k1,k2, M0,M1),NCOMP)];
71 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,k1+1,k2, M0,M1),NCOMP)];
72 register const double f_010 = in[INDEX2(i,INDEX3(k0,k1+1,k2, M0,M1),NCOMP)];
73 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,k1+1,k2+1, M0,M1),NCOMP)];
74 out[INDEX3(i,0,INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_000 + f_001 + f_010 + f_011 + f_100 + f_101 + f_110 + f_111);
75 } /* close component loop i */
76 } /* close k0 loop */
77 } /* close k1 loop */
78 } /* close k2 loop */
79 } else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {
80 if (face_offset(0)>-1) {
81 const double tmp0_0 = 0.78729833462074168852;
82 const double tmp0_4 = 0.056350832689629155741;
83 const double tmp0_1 = 0.10000000000000000000;
84 const double tmp0_5 = 0.25000000000000000000;
85 const double tmp0_2 = 0.012701665379258311482;
86 const double tmp0_3 = 0.44364916731037084426;
87 #pragma omp parallel for private(i,k2,k1)
88 for (k2 =0; k2 < N2; ++k2) {
89 for (k1 =0; k1 < N1; ++k1) {
90 for (i =0; i < NCOMP; ++i) {
91 register const double f_000 = in[INDEX2(i,INDEX3(0,k1,k2, M0,M1),NCOMP)];
92 register const double f_001 = in[INDEX2(i,INDEX3(0,k1,k2+1, M0,M1),NCOMP)];
93 register const double f_011 = in[INDEX2(i,INDEX3(0,k1+1,k2+1, M0,M1),NCOMP)];
94 register const double f_010 = in[INDEX2(i,INDEX3(0,k1+1,k2, M0,M1),NCOMP)];
95 out[INDEX3(i,0,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_0 + f_011*tmp0_2 + tmp0_1*(f_001 + f_010);
96 out[INDEX3(i,1,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_010) + tmp0_4*(f_001 + f_011);
97 out[INDEX3(i,2,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_2 + f_010*tmp0_0 + tmp0_1*(f_000 + f_011);
98 out[INDEX3(i,3,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_001) + tmp0_4*(f_010 + f_011);
99 out[INDEX3(i,4,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_000 + f_001 + f_010 + f_011);
100 out[INDEX3(i,5,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_010 + f_011) + tmp0_4*(f_000 + f_001);
101 out[INDEX3(i,6,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_0 + f_010*tmp0_2 + tmp0_1*(f_000 + f_011);
102 out[INDEX3(i,7,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_001 + f_011) + tmp0_4*(f_000 + f_010);
103 out[INDEX3(i,8,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_2 + f_011*tmp0_0 + tmp0_1*(f_001 + f_010);
104 } /* close component loop i */
105 } /* close k1 loop */
106 } /* close k2 loop */
107 } /* end of face 0 */
108 if (face_offset(1)>-1) {
109 const double tmp0_0 = 0.10000000000000000000;
110 const double tmp0_4 = 0.44364916731037084426;
111 const double tmp0_1 = 0.78729833462074168852;
112 const double tmp0_5 = 0.25000000000000000000;
113 const double tmp0_2 = 0.012701665379258311482;
114 const double tmp0_3 = 0.056350832689629155741;
115 #pragma omp parallel for private(i,k2,k1)
116 for (k2 =0; k2 < N2; ++k2) {
117 for (k1 =0; k1 < N1; ++k1) {
118 for (i =0; i < NCOMP; ++i) {
119 register const double f_101 = in[INDEX2(i,INDEX3(M0-1,k1,k2+1, M0,M1),NCOMP)];
120 register const double f_100 = in[INDEX2(i,INDEX3(M0-1,k1,k2, M0,M1),NCOMP)];
121 register const double f_110 = in[INDEX2(i,INDEX3(M0-1,k1+1,k2, M0,M1),NCOMP)];
122 register const double f_111 = in[INDEX2(i,INDEX3(M0-1,k1+1,k2+1, M0,M1),NCOMP)];
123 out[INDEX3(i,0,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_100*tmp0_1 + f_111*tmp0_2 + tmp0_0*(f_101 + f_110);
124 out[INDEX3(i,1,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_101 + f_111) + tmp0_4*(f_100 + f_110);
125 out[INDEX3(i,2,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_101*tmp0_2 + f_110*tmp0_1 + tmp0_0*(f_100 + f_111);
126 out[INDEX3(i,3,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_110 + f_111) + tmp0_4*(f_100 + f_101);
127 out[INDEX3(i,4,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_100 + f_101 + f_110 + f_111);
128 out[INDEX3(i,5,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_100 + f_101) + tmp0_4*(f_110 + f_111);
129 out[INDEX3(i,6,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_101*tmp0_1 + f_110*tmp0_2 + tmp0_0*(f_100 + f_111);
130 out[INDEX3(i,7,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_100 + f_110) + tmp0_4*(f_101 + f_111);
131 out[INDEX3(i,8,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_100*tmp0_2 + f_111*tmp0_1 + tmp0_0*(f_101 + f_110);
132 } /* close component loop i */
133 } /* close k1 loop */
134 } /* close k2 loop */
135 } /* end of face 1 */
136 if (face_offset(2)>-1) {
137 const double tmp0_0 = 0.78729833462074168852;
138 const double tmp0_4 = 0.056350832689629155741;
139 const double tmp0_1 = 0.10000000000000000000;
140 const double tmp0_5 = 0.25000000000000000000;
141 const double tmp0_2 = 0.012701665379258311482;
142 const double tmp0_3 = 0.44364916731037084426;
143 #pragma omp parallel for private(i,k2,k0)
144 for (k2 =0; k2 < N2; ++k2) {
145 for (k0 =0; k0 < N0; ++k0) {
146 for (i =0; i < NCOMP; ++i) {
147 register const double f_000 = in[INDEX2(i,INDEX3(k0,0,k2, M0,M1),NCOMP)];
148 register const double f_001 = in[INDEX2(i,INDEX3(k0,0,k2+1, M0,M1),NCOMP)];
149 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,0,k2+1, M0,M1),NCOMP)];
150 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,0,k2, M0,M1),NCOMP)];
151 out[INDEX3(i,0,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_0 + f_101*tmp0_2 + tmp0_1*(f_001 + f_100);
152 out[INDEX3(i,1,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_100) + tmp0_4*(f_001 + f_101);
153 out[INDEX3(i,2,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_2 + f_100*tmp0_0 + tmp0_1*(f_000 + f_101);
154 out[INDEX3(i,3,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_001) + tmp0_4*(f_100 + f_101);
155 out[INDEX3(i,4,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_000 + f_001 + f_100 + f_101);
156 out[INDEX3(i,5,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_100 + f_101) + tmp0_4*(f_000 + f_001);
157 out[INDEX3(i,6,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_0 + f_100*tmp0_2 + tmp0_1*(f_000 + f_101);
158 out[INDEX3(i,7,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_001 + f_101) + tmp0_4*(f_000 + f_100);
159 out[INDEX3(i,8,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_2 + f_101*tmp0_0 + tmp0_1*(f_001 + f_100);
160 } /* close component loop i */
161 } /* close k0 loop */
162 } /* close k2 loop */
163 } /* end of face 2 */
164 if (face_offset(3)>-1) {
165 const double tmp0_0 = 0.10000000000000000000;
166 const double tmp0_4 = 0.44364916731037084426;
167 const double tmp0_1 = 0.78729833462074168852;
168 const double tmp0_5 = 0.25000000000000000000;
169 const double tmp0_2 = 0.012701665379258311482;
170 const double tmp0_3 = 0.056350832689629155741;
171 #pragma omp parallel for private(i,k2,k0)
172 for (k2 =0; k2 < N2; ++k2) {
173 for (k0 =0; k0 < N0; ++k0) {
174 for (i =0; i < NCOMP; ++i) {
175 register const double f_011 = in[INDEX2(i,INDEX3(k0,M1-1,k2+1, M0,M1),NCOMP)];
176 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,M1-1,k2, M0,M1),NCOMP)];
177 register const double f_010 = in[INDEX2(i,INDEX3(k0,M1-1,k2, M0,M1),NCOMP)];
178 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,M1-1,k2+1, M0,M1),NCOMP)];
179 out[INDEX3(i,0,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_010*tmp0_1 + f_111*tmp0_2 + tmp0_0*(f_011 + f_110);
180 out[INDEX3(i,1,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_011 + f_111) + tmp0_4*(f_010 + f_110);
181 out[INDEX3(i,2,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_011*tmp0_2 + f_110*tmp0_1 + tmp0_0*(f_010 + f_111);
182 out[INDEX3(i,3,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_110 + f_111) + tmp0_4*(f_010 + f_011);
183 out[INDEX3(i,4,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_010 + f_011 + f_110 + f_111);
184 out[INDEX3(i,5,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_010 + f_011) + tmp0_4*(f_110 + f_111);
185 out[INDEX3(i,6,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_011*tmp0_1 + f_110*tmp0_2 + tmp0_0*(f_010 + f_111);
186 out[INDEX3(i,7,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_010 + f_110) + tmp0_4*(f_011 + f_111);
187 out[INDEX3(i,8,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_010*tmp0_2 + f_111*tmp0_1 + tmp0_0*(f_011 + f_110);
188 } /* close component loop i */
189 } /* close k0 loop */
190 } /* close k2 loop */
191 } /* end of face 3 */
192 if (face_offset(4)>-1) {
193 const double tmp0_0 = 0.78729833462074168852;
194 const double tmp0_4 = 0.056350832689629155741;
195 const double tmp0_1 = 0.10000000000000000000;
196 const double tmp0_5 = 0.25000000000000000000;
197 const double tmp0_2 = 0.012701665379258311482;
198 const double tmp0_3 = 0.44364916731037084426;
199 #pragma omp parallel for private(i,k1,k0)
200 for (k1 =0; k1 < N1; ++k1) {
201 for (k0 =0; k0 < N0; ++k0) {
202 for (i =0; i < NCOMP; ++i) {
203 register const double f_000 = in[INDEX2(i,INDEX3(k0,k1,0, M0,M1),NCOMP)];
204 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,k1,0, M0,M1),NCOMP)];
205 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,k1+1,0, M0,M1),NCOMP)];
206 register const double f_010 = in[INDEX2(i,INDEX3(k0,k1+1,0, M0,M1),NCOMP)];
207 out[INDEX3(i,0,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_0 + f_110*tmp0_2 + tmp0_1*(f_010 + f_100);
208 out[INDEX3(i,1,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_100) + tmp0_4*(f_010 + f_110);
209 out[INDEX3(i,2,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_010*tmp0_2 + f_100*tmp0_0 + tmp0_1*(f_000 + f_110);
210 out[INDEX3(i,3,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_000 + f_010) + tmp0_4*(f_100 + f_110);
211 out[INDEX3(i,4,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_000 + f_010 + f_100 + f_110);
212 out[INDEX3(i,5,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_100 + f_110) + tmp0_4*(f_000 + f_010);
213 out[INDEX3(i,6,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_010*tmp0_0 + f_100*tmp0_2 + tmp0_1*(f_000 + f_110);
214 out[INDEX3(i,7,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_010 + f_110) + tmp0_4*(f_000 + f_100);
215 out[INDEX3(i,8,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_000*tmp0_2 + f_110*tmp0_0 + tmp0_1*(f_010 + f_100);
216 } /* close component loop i */
217 } /* close k0 loop */
218 } /* close k1 loop */
219 } /* end of face 4 */
220 if (face_offset(5)>-1) {
221 const double tmp0_0 = 0.78729833462074168852;
222 const double tmp0_4 = 0.056350832689629155741;
223 const double tmp0_1 = 0.10000000000000000000;
224 const double tmp0_5 = 0.25000000000000000000;
225 const double tmp0_2 = 0.012701665379258311482;
226 const double tmp0_3 = 0.44364916731037084426;
227 #pragma omp parallel for private(i,k1,k0)
228 for (k1 =0; k1 < N1; ++k1) {
229 for (k0 =0; k0 < N0; ++k0) {
230 for (i =0; i < NCOMP; ++i) {
231 register const double f_001 = in[INDEX2(i,INDEX3(k0,k1,M2-1, M0,M1),NCOMP)];
232 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,k1,M2-1, M0,M1),NCOMP)];
233 register const double f_011 = in[INDEX2(i,INDEX3(k0,k1+1,M2-1, M0,M1),NCOMP)];
234 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,k1+1,M2-1, M0,M1),NCOMP)];
235 out[INDEX3(i,0,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_0 + f_111*tmp0_2 + tmp0_1*(f_011 + f_101);
236 out[INDEX3(i,1,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_001 + f_101) + tmp0_4*(f_011 + f_111);
237 out[INDEX3(i,2,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_011*tmp0_2 + f_101*tmp0_0 + tmp0_1*(f_001 + f_111);
238 out[INDEX3(i,3,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_001 + f_011) + tmp0_4*(f_101 + f_111);
239 out[INDEX3(i,4,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_5*(f_001 + f_011 + f_101 + f_111);
240 out[INDEX3(i,5,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_101 + f_111) + tmp0_4*(f_001 + f_011);
241 out[INDEX3(i,6,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_011*tmp0_0 + f_101*tmp0_2 + tmp0_1*(f_001 + f_111);
242 out[INDEX3(i,7,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = tmp0_3*(f_011 + f_111) + tmp0_4*(f_001 + f_101);
243 out[INDEX3(i,8,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,9)] = f_001*tmp0_2 + f_111*tmp0_0 + tmp0_1*(f_011 + f_101);
244 } /* close component loop i */
245 } /* close k0 loop */
246 } /* close k1 loop */
247 } /* end of face 5 */
248 } else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {
249 if (face_offset(0)>-1) {
250 const double tmp0_0 = 0.25000000000000000000;
251 #pragma omp parallel for private(i,k2,k1)
252 for (k2 =0; k2 < N2; ++k2) {
253 for (k1 =0; k1 < N1; ++k1) {
254 for (i =0; i < NCOMP; ++i) {
255 register const double f_011 = in[INDEX2(i,INDEX3(0,k1+1,k2+1, M0,M1),NCOMP)];
256 register const double f_010 = in[INDEX2(i,INDEX3(0,k1+1,k2, M0,M1),NCOMP)];
257 register const double f_001 = in[INDEX2(i,INDEX3(0,k1,k2+1, M0,M1),NCOMP)];
258 register const double f_000 = in[INDEX2(i,INDEX3(0,k1,k2, M0,M1),NCOMP)];
259 out[INDEX3(i,0,face_offset(0)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_000 + f_001 + f_010 + f_011);
260 } /* close component loop i */
261 } /* close k1 loop */
262 } /* close k2 loop */
263 } /* end of face 0 */
264 if (face_offset(1)>-1) {
265 const double tmp0_0 = 0.25000000000000000000;
266 #pragma omp parallel for private(i,k2,k1)
267 for (k2 =0; k2 < N2; ++k2) {
268 for (k1 =0; k1 < N1; ++k1) {
269 for (i =0; i < NCOMP; ++i) {
270 register const double f_110 = in[INDEX2(i,INDEX3(M0-1,k1+1,k2, M0,M1),NCOMP)];
271 register const double f_100 = in[INDEX2(i,INDEX3(M0-1,k1,k2, M0,M1),NCOMP)];
272 register const double f_101 = in[INDEX2(i,INDEX3(M0-1,k1,k2+1, M0,M1),NCOMP)];
273 register const double f_111 = in[INDEX2(i,INDEX3(M0-1,k1+1,k2+1, M0,M1),NCOMP)];
274 out[INDEX3(i,0,face_offset(1)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_100 + f_101 + f_110 + f_111);
275 } /* close component loop i */
276 } /* close k1 loop */
277 } /* close k2 loop */
278 } /* end of face 1 */
279 if (face_offset(2)>-1) {
280 const double tmp0_0 = 0.25000000000000000000;
281 #pragma omp parallel for private(i,k2,k0)
282 for (k2 =0; k2 < N2; ++k2) {
283 for (k0 =0; k0 < N0; ++k0) {
284 for (i =0; i < NCOMP; ++i) {
285 register const double f_000 = in[INDEX2(i,INDEX3(k0,0,k2, M0,M1),NCOMP)];
286 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,0,k2, M0,M1),NCOMP)];
287 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,0,k2+1, M0,M1),NCOMP)];
288 register const double f_001 = in[INDEX2(i,INDEX3(k0,0,k2+1, M0,M1),NCOMP)];
289 out[INDEX3(i,0,face_offset(2)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_000 + f_001 + f_100 + f_101);
290 } /* close component loop i */
291 } /* close k0 loop */
292 } /* close k2 loop */
293 } /* end of face 2 */
294 if (face_offset(3)>-1) {
295 const double tmp0_0 = 0.25000000000000000000;
296 #pragma omp parallel for private(i,k2,k0)
297 for (k2 =0; k2 < N2; ++k2) {
298 for (k0 =0; k0 < N0; ++k0) {
299 for (i =0; i < NCOMP; ++i) {
300 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,M1-1,k2, M0,M1),NCOMP)];
301 register const double f_011 = in[INDEX2(i,INDEX3(k0,M1-1,k2+1, M0,M1),NCOMP)];
302 register const double f_010 = in[INDEX2(i,INDEX3(k0,M1-1,k2, M0,M1),NCOMP)];
303 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,M1-1,k2+1, M0,M1),NCOMP)];
304 out[INDEX3(i,0,face_offset(3)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_010 + f_011 + f_110 + f_111);
305 } /* close component loop i */
306 } /* close k0 loop */
307 } /* close k2 loop */
308 } /* end of face 3 */
309 if (face_offset(4)>-1) {
310 const double tmp0_0 = 0.25000000000000000000;
311 #pragma omp parallel for private(i,k1,k0)
312 for (k1 =0; k1 < N1; ++k1) {
313 for (k0 =0; k0 < N0; ++k0) {
314 for (i =0; i < NCOMP; ++i) {
315 register const double f_110 = in[INDEX2(i,INDEX3(k0+1,k1+1,0, M0,M1),NCOMP)];
316 register const double f_010 = in[INDEX2(i,INDEX3(k0,k1+1,0, M0,M1),NCOMP)];
317 register const double f_100 = in[INDEX2(i,INDEX3(k0+1,k1,0, M0,M1),NCOMP)];
318 register const double f_000 = in[INDEX2(i,INDEX3(k0,k1,0, M0,M1),NCOMP)];
319 out[INDEX3(i,0,face_offset(4)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_000 + f_010 + f_100 + f_110);
320 } /* close component loop i */
321 } /* close k0 loop */
322 } /* close k1 loop */
323 } /* end of face 4 */
324 if (face_offset(5)>-1) {
325 const double tmp0_0 = 0.25000000000000000000;
326 #pragma omp parallel for private(i,k1,k0)
327 for (k1 =0; k1 < N1; ++k1) {
328 for (k0 =0; k0 < N0; ++k0) {
329 for (i =0; i < NCOMP; ++i) {
330 register const double f_011 = in[INDEX2(i,INDEX3(k0,k1+1,M2-1, M0,M1),NCOMP)];
331 register const double f_001 = in[INDEX2(i,INDEX3(k0,k1,M2-1, M0,M1),NCOMP)];
332 register const double f_101 = in[INDEX2(i,INDEX3(k0+1,k1,M2-1, M0,M1),NCOMP)];
333 register const double f_111 = in[INDEX2(i,INDEX3(k0+1,k1+1,M2-1, M0,M1),NCOMP)];
334 out[INDEX3(i,0,face_offset(5)+INDEX3(k0,k1,k2,N0,N1),NCOMP,1)] = tmp0_0*(f_001 + f_011 + f_101 + f_111);
335 } /* close component loop i */
336 } /* close k0 loop */
337 } /* close k1 loop */
338 } /* end of face 5 */
339 } /* end of out_data_type branching
340 /* GENERATOR SNIP BOTTOM */
341 }

  ViewVC Help
Powered by ViewVC 1.1.26