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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3641 - (show annotations)
Thu Oct 27 02:16:12 2011 UTC (7 years, 11 months ago) by gross
File MIME type: text/plain
File size: 16760 byte(s)
more work on the assemblage
1 Ripley_Assemble_Integration_2D(Ripley_Grid *grid, Escript in, double *out)
2 {
3 /* GENERATOR SNIP TOP */
4 if (out_data_type==RIPLEY_ELEMENTS) {
5 const double tmp0_16 = -0.88729833462074168852/h1;
6 const double tmp0_18 = 0.11270166537925831148/h1;
7 const double tmp0_13 = -0.88729833462074168852/h1;
8 const double tmp0_0 = 0.88729833462074168852/h0;
9 const double tmp0_17 = 0.88729833462074168852/h1;
10 const double tmp0_4 = 0.5/h0;
11 const double tmp0_19 = -0.11270166537925831148/h1;
12 const double tmp0_10 = -0.11270166537925831148/h1;
13 const double tmp0_1 = 0.11270166537925831148/h0;
14 const double tmp0_8 = -0.88729833462074168852/h0;
15 const double tmp0_14 = -0.5/h1;
16 const double tmp0_5 = -0.5/h0;
17 const double tmp0_11 = 0.11270166537925831148/h1;
18 const double tmp0_2 = -0.11270166537925831148/h0;
19 const double tmp0_9 = -0.11270166537925831148/h0;
20 const double tmp0_15 = 0.5/h1;
21 const double tmp0_6 = 0.11270166537925831148/h0;
22 const double tmp0_3 = -0.88729833462074168852/h0;
23 const double tmp0_12 = 0.88729833462074168852/h1;
24 const double tmp0_7 = 0.88729833462074168852/h0;
25 #pragma omp parallel for private(i,k1,k0)
26 for (k1 =0; k1 < N1; ++k1) {
27 for (k0 =0; k0 < N0; ++k0) {
28 for (i =0; i < NCOMP; ++i) {
29 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,k1, M0),NCOMP)];
30 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,k1+1, M0),NCOMP)];
31 const register double f_01 = in[INDEX2(i,INDEX2(k0,k1+1, M0),NCOMP)];
32 const register double f_00 = in[INDEX2(i,INDEX2(k0,k1, M0),NCOMP)];
33 out[INDEX4(i,0,0,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_3 + f_01*tmp0_2 + f_10*tmp0_0 + f_11*tmp0_1;
34 out[INDEX4(i,1,0,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_13 + f_01*tmp0_12 + f_10*tmp0_10 + f_11*tmp0_11;
35 out[INDEX4(i,0,1,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_3 + f_01*tmp0_2 + f_10*tmp0_0 + f_11*tmp0_1;
36 out[INDEX4(i,1,1,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_14*(f_00 + f_10) + tmp0_15*(f_01 + f_11);
37 out[INDEX4(i,0,2,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_3 + f_01*tmp0_2 + f_10*tmp0_0 + f_11*tmp0_1;
38 out[INDEX4(i,1,2,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_19 + f_01*tmp0_18 + f_10*tmp0_16 + f_11*tmp0_17;
39 out[INDEX4(i,0,3,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_4*(f_10 + f_11) + tmp0_5*(f_00 + f_01);
40 out[INDEX4(i,1,3,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_13 + f_01*tmp0_12 + f_10*tmp0_10 + f_11*tmp0_11;
41 out[INDEX4(i,0,4,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_4*(f_10 + f_11) + tmp0_5*(f_00 + f_01);
42 out[INDEX4(i,1,4,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_14*(f_00 + f_10) + tmp0_15*(f_01 + f_11);
43 out[INDEX4(i,0,5,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_4*(f_10 + f_11) + tmp0_5*(f_00 + f_01);
44 out[INDEX4(i,1,5,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_19 + f_01*tmp0_18 + f_10*tmp0_16 + f_11*tmp0_17;
45 out[INDEX4(i,0,6,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_9 + f_01*tmp0_8 + f_10*tmp0_6 + f_11*tmp0_7;
46 out[INDEX4(i,1,6,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_13 + f_01*tmp0_12 + f_10*tmp0_10 + f_11*tmp0_11;
47 out[INDEX4(i,0,7,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_9 + f_01*tmp0_8 + f_10*tmp0_6 + f_11*tmp0_7;
48 out[INDEX4(i,1,7,INDEX2(k0,k1,N0),NCOMP,2,9)] = tmp0_14*(f_00 + f_10) + tmp0_15*(f_01 + f_11);
49 out[INDEX4(i,0,8,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_9 + f_01*tmp0_8 + f_10*tmp0_6 + f_11*tmp0_7;
50 out[INDEX4(i,1,8,INDEX2(k0,k1,N0),NCOMP,2,9)] = f_00*tmp0_19 + f_01*tmp0_18 + f_10*tmp0_16 + f_11*tmp0_17;
51 } /* close component loop i */
52 } /* close k0 loop */
53 } /* close k1 loop */
54 } else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {
55 const double tmp0_3 = 0.5/h1;
56 const double tmp0_2 = -0.5/h1;
57 const double tmp0_1 = -0.5/h0;
58 const double tmp0_0 = 0.5/h0;
59 #pragma omp parallel for private(i,k1,k0)
60 for (k1 =0; k1 < N1; ++k1) {
61 for (k0 =0; k0 < N0; ++k0) {
62 for (i =0; i < NCOMP; ++i) {
63 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,k1, M0),NCOMP)];
64 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,k1+1, M0),NCOMP)];
65 const register double f_01 = in[INDEX2(i,INDEX2(k0,k1+1, M0),NCOMP)];
66 const register double f_00 = in[INDEX2(i,INDEX2(k0,k1, M0),NCOMP)];
67 out[INDEX4(i,0,0,INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_0*(f_10 + f_11) + tmp0_1*(f_00 + f_01);
68 out[INDEX4(i,1,0,INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_2*(f_00 + f_10) + tmp0_3*(f_01 + f_11);
69 } /* close component loop i */
70 } /* close k0 loop */
71 } /* close k1 loop */
72 } else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {
73 if (face_offset(0)>-1) {
74 const double tmp0_0 = 0.88729833462074168852/h0;
75 const double tmp0_4 = 0.5/h0;
76 const double tmp0_10 = 1.0/h1;
77 const double tmp0_1 = 0.11270166537925831148/h0;
78 const double tmp0_8 = -0.88729833462074168852/h0;
79 const double tmp0_5 = -0.5/h0;
80 const double tmp0_11 = -1/h1;
81 const double tmp0_2 = -0.11270166537925831148/h0;
82 const double tmp0_9 = -0.11270166537925831148/h0;
83 const double tmp0_6 = 0.11270166537925831148/h0;
84 const double tmp0_3 = -0.88729833462074168852/h0;
85 const double tmp0_7 = 0.88729833462074168852/h0;
86 #pragma omp parallel for private(i,k1)
87 for (k1 =0; k1 < N1; ++k1) {
88 for (i =0; i < NCOMP; ++i) {
89 const register double f_10 = in[INDEX2(i,INDEX2(1,k1, M0),NCOMP)];
90 const register double f_11 = in[INDEX2(i,INDEX2(1,k1+1, M0),NCOMP)];
91 const register double f_01 = in[INDEX2(i,INDEX2(0,k1+1, M0),NCOMP)];
92 const register double f_00 = in[INDEX2(i,INDEX2(0,k1, M0),NCOMP)];
93 out[INDEX4(i,0,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_3 + f_01*tmp0_2 + f_10*tmp0_0 + f_11*tmp0_1;
94 out[INDEX4(i,1,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_11 + f_01*tmp0_10;
95 out[INDEX4(i,0,1,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = tmp0_4*(f_10 + f_11) + tmp0_5*(f_00 + f_01);
96 out[INDEX4(i,1,1,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_11 + f_01*tmp0_10;
97 out[INDEX4(i,0,2,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_9 + f_01*tmp0_8 + f_10*tmp0_6 + f_11*tmp0_7;
98 out[INDEX4(i,1,2,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_11 + f_01*tmp0_10;
99 } /* close component loop i */
100 } /* close k1 loop */
101 } /* end of face 0 */
102 if (face_offset(1)>-1) {
103 const double tmp0_0 = 0.88729833462074168852/h0;
104 const double tmp0_4 = 0.5/h0;
105 const double tmp0_10 = -1/h1;
106 const double tmp0_1 = 0.11270166537925831148/h0;
107 const double tmp0_8 = -0.88729833462074168852/h0;
108 const double tmp0_5 = -0.5/h0;
109 const double tmp0_11 = 1.0/h1;
110 const double tmp0_2 = -0.11270166537925831148/h0;
111 const double tmp0_9 = -0.11270166537925831148/h0;
112 const double tmp0_6 = 0.11270166537925831148/h0;
113 const double tmp0_3 = -0.88729833462074168852/h0;
114 const double tmp0_7 = 0.88729833462074168852/h0;
115 #pragma omp parallel for private(i,k1)
116 for (k1 =0; k1 < N1; ++k1) {
117 for (i =0; i < NCOMP; ++i) {
118 const register double f_10 = in[INDEX2(i,INDEX2(M0-1,k1, M0),NCOMP)];
119 const register double f_11 = in[INDEX2(i,INDEX2(M0-1,k1+1, M0),NCOMP)];
120 const register double f_01 = in[INDEX2(i,INDEX2(M0-2,k1+1, M0),NCOMP)];
121 const register double f_00 = in[INDEX2(i,INDEX2(M0-2,k1, M0),NCOMP)];
122 out[INDEX4(i,0,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_3 + f_01*tmp0_2 + f_10*tmp0_0 + f_11*tmp0_1;
123 out[INDEX4(i,1,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_10*tmp0_10 + f_11*tmp0_11;
124 out[INDEX4(i,0,1,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = tmp0_4*(f_10 + f_11) + tmp0_5*(f_00 + f_01);
125 out[INDEX4(i,1,1,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_10*tmp0_10 + f_11*tmp0_11;
126 out[INDEX4(i,0,2,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_9 + f_01*tmp0_8 + f_10*tmp0_6 + f_11*tmp0_7;
127 out[INDEX4(i,1,2,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_10*tmp0_10 + f_11*tmp0_11;
128 } /* close component loop i */
129 } /* close k1 loop */
130 } /* end of face 1 */
131 if (face_offset(2)>-1) {
132 const double tmp0_0 = 1.0/h0;
133 const double tmp0_4 = 0.11270166537925831148/h1;
134 const double tmp0_10 = 0.88729833462074168852/h1;
135 const double tmp0_1 = -1/h0;
136 const double tmp0_8 = -0.88729833462074168852/h1;
137 const double tmp0_5 = 0.88729833462074168852/h1;
138 const double tmp0_11 = 0.11270166537925831148/h1;
139 const double tmp0_2 = -0.11270166537925831148/h1;
140 const double tmp0_9 = -0.11270166537925831148/h1;
141 const double tmp0_6 = -0.5/h1;
142 const double tmp0_3 = -0.88729833462074168852/h1;
143 const double tmp0_7 = 0.5/h1;
144 #pragma omp parallel for private(i,k0)
145 for (k0 =0; k0 < N0; ++k0) {
146 for (i =0; i < NCOMP; ++i) {
147 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,0, M0),NCOMP)];
148 const register double f_00 = in[INDEX2(i,INDEX2(k0,0, M0),NCOMP)];
149 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,1, M0),NCOMP)];
150 const register double f_01 = in[INDEX2(i,INDEX2(k0,1, M0),NCOMP)];
151 out[INDEX4(i,0,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_1 + f_10*tmp0_0;
152 out[INDEX4(i,1,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_3 + f_01*tmp0_5 + f_10*tmp0_2 + f_11*tmp0_4;
153 out[INDEX4(i,0,1,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_1 + f_10*tmp0_0;
154 out[INDEX4(i,1,1,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = tmp0_6*(f_00 + f_10) + tmp0_7*(f_01 + f_11);
155 out[INDEX4(i,0,2,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_1 + f_10*tmp0_0;
156 out[INDEX4(i,1,2,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_9 + f_01*tmp0_11 + f_10*tmp0_8 + f_11*tmp0_10;
157 } /* close component loop i */
158 } /* close k0 loop */
159 } /* end of face 2 */
160 if (face_offset(3)>-1) {
161 const double tmp0_0 = 1.0/h0;
162 const double tmp0_4 = -0.11270166537925831148/h1;
163 const double tmp0_10 = -0.88729833462074168852/h1;
164 const double tmp0_1 = -1/h0;
165 const double tmp0_8 = 0.88729833462074168852/h1;
166 const double tmp0_5 = -0.88729833462074168852/h1;
167 const double tmp0_11 = -0.11270166537925831148/h1;
168 const double tmp0_2 = 0.11270166537925831148/h1;
169 const double tmp0_9 = 0.11270166537925831148/h1;
170 const double tmp0_6 = 0.5/h1;
171 const double tmp0_3 = 0.88729833462074168852/h1;
172 const double tmp0_7 = -0.5/h1;
173 #pragma omp parallel for private(i,k0)
174 for (k0 =0; k0 < N0; ++k0) {
175 for (i =0; i < NCOMP; ++i) {
176 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,M1-1, M0),NCOMP)];
177 const register double f_01 = in[INDEX2(i,INDEX2(k0,M1-1, M0),NCOMP)];
178 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,M1-2, M0),NCOMP)];
179 const register double f_00 = in[INDEX2(i,INDEX2(k0,M1-2, M0),NCOMP)];
180 out[INDEX4(i,0,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_01*tmp0_1 + f_11*tmp0_0;
181 out[INDEX4(i,1,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_5 + f_01*tmp0_3 + f_10*tmp0_4 + f_11*tmp0_2;
182 out[INDEX4(i,0,1,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_01*tmp0_1 + f_11*tmp0_0;
183 out[INDEX4(i,1,1,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = tmp0_6*(f_01 + f_11) + tmp0_7*(f_00 + f_10);
184 out[INDEX4(i,0,2,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_01*tmp0_1 + f_11*tmp0_0;
185 out[INDEX4(i,1,2,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,3)] = f_00*tmp0_11 + f_01*tmp0_9 + f_10*tmp0_10 + f_11*tmp0_8;
186 } /* close component loop i */
187 } /* close k0 loop */
188 } /* end of face 3 */
189 } else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {
190 if (face_offset(0)>-1) {
191 const double tmp0_3 = -1/h1;
192 const double tmp0_2 = 1.0/h1;
193 const double tmp0_1 = -0.5/h0;
194 const double tmp0_0 = 0.5/h0;
195 #pragma omp parallel for private(i,k1)
196 for (k1 =0; k1 < N1; ++k1) {
197 for (i =0; i < NCOMP; ++i) {
198 const register double f_10 = in[INDEX2(i,INDEX2(1,k1, M0),NCOMP)];
199 const register double f_11 = in[INDEX2(i,INDEX2(1,k1+1, M0),NCOMP)];
200 const register double f_01 = in[INDEX2(i,INDEX2(0,k1+1, M0),NCOMP)];
201 const register double f_00 = in[INDEX2(i,INDEX2(0,k1, M0),NCOMP)];
202 out[INDEX4(i,0,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_0*(f_10 + f_11) + tmp0_1*(f_00 + f_01);
203 out[INDEX4(i,1,0,face_offset(0)+INDEX2(k0,k1,N0),NCOMP,2,1)] = f_00*tmp0_3 + f_01*tmp0_2;
204 } /* close component loop i */
205 } /* close k1 loop */
206 } /* end of face 0 */
207 if (face_offset(1)>-1) {
208 const double tmp0_3 = 1.0/h1;
209 const double tmp0_2 = -1/h1;
210 const double tmp0_1 = -0.5/h0;
211 const double tmp0_0 = 0.5/h0;
212 #pragma omp parallel for private(i,k1)
213 for (k1 =0; k1 < N1; ++k1) {
214 for (i =0; i < NCOMP; ++i) {
215 const register double f_10 = in[INDEX2(i,INDEX2(M0-1,k1, M0),NCOMP)];
216 const register double f_11 = in[INDEX2(i,INDEX2(M0-1,k1+1, M0),NCOMP)];
217 const register double f_01 = in[INDEX2(i,INDEX2(M0-2,k1+1, M0),NCOMP)];
218 const register double f_00 = in[INDEX2(i,INDEX2(M0-2,k1, M0),NCOMP)];
219 out[INDEX4(i,0,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_0*(f_10 + f_11) + tmp0_1*(f_00 + f_01);
220 out[INDEX4(i,1,0,face_offset(1)+INDEX2(k0,k1,N0),NCOMP,2,1)] = f_10*tmp0_2 + f_11*tmp0_3;
221 } /* close component loop i */
222 } /* close k1 loop */
223 } /* end of face 1 */
224 if (face_offset(2)>-1) {
225 const double tmp0_3 = 0.5/h1;
226 const double tmp0_2 = -0.5/h1;
227 const double tmp0_1 = -1/h0;
228 const double tmp0_0 = 1.0/h0;
229 #pragma omp parallel for private(i,k0)
230 for (k0 =0; k0 < N0; ++k0) {
231 for (i =0; i < NCOMP; ++i) {
232 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,0, M0),NCOMP)];
233 const register double f_00 = in[INDEX2(i,INDEX2(k0,0, M0),NCOMP)];
234 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,1, M0),NCOMP)];
235 const register double f_01 = in[INDEX2(i,INDEX2(k0,1, M0),NCOMP)];
236 out[INDEX4(i,0,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,1)] = f_00*tmp0_1 + f_10*tmp0_0;
237 out[INDEX4(i,1,0,face_offset(2)+INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_2*(f_00 + f_10) + tmp0_3*(f_01 + f_11);
238 } /* close component loop i */
239 } /* close k0 loop */
240 } /* end of face 2 */
241 if (face_offset(3)>-1) {
242 const double tmp0_3 = -0.5/h1;
243 const double tmp0_2 = 0.5/h1;
244 const double tmp0_1 = -1/h0;
245 const double tmp0_0 = 1.0/h0;
246 #pragma omp parallel for private(i,k0)
247 for (k0 =0; k0 < N0; ++k0) {
248 for (i =0; i < NCOMP; ++i) {
249 const register double f_11 = in[INDEX2(i,INDEX2(k0+1,M1-1, M0),NCOMP)];
250 const register double f_01 = in[INDEX2(i,INDEX2(k0,M1-1, M0),NCOMP)];
251 const register double f_10 = in[INDEX2(i,INDEX2(k0+1,M1-2, M0),NCOMP)];
252 const register double f_00 = in[INDEX2(i,INDEX2(k0,M1-2, M0),NCOMP)];
253 out[INDEX4(i,0,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,1)] = f_01*tmp0_1 + f_11*tmp0_0;
254 out[INDEX4(i,1,0,face_offset(3)+INDEX2(k0,k1,N0),NCOMP,2,1)] = tmp0_2*(f_01 + f_11) + tmp0_3*(f_00 + f_10);
255 } /* close component loop i */
256 } /* close k0 loop */
257 } /* end of face 3 */
258 } /* end of out_data_type branching
259 /* GENERATOR SNIP BOTTOM */
260 }

  ViewVC Help
Powered by ViewVC 1.1.26