/[escript]/trunk/dudley/src/Assemble_PDE_Single2_3D.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_PDE_Single2_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 14532 byte(s)
Round 1 of copyright fixes
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 /************************************************************************************/
17
18 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
19 /* the shape functions for test and solution must be identical */
20
21 /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
22
23 /* in a 3D domain. The shape functions for test and solution must be identical */
24 /* and row_NS == row_NN */
25
26 /* Shape of the coefficients: */
27
28 /* A = 3 x 3 */
29 /* B = 3 */
30 /* C = 3 */
31 /* D = scalar */
32 /* X = 3 */
33 /* Y = scalar */
34
35 /************************************************************************************/
36
37 #include "Assemble.h"
38 #include "Util.h"
39 #ifdef _OPENMP
40 #include <omp.h>
41 #endif
42
43 /************************************************************************************/
44
45 void Dudley_Assemble_PDE_Single2_3D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46 Paso_SystemMatrix * Mat, escriptDataC * F,
47 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48 escriptDataC * X, escriptDataC * Y)
49 {
50
51 #define DIM 3
52 index_t color;
53 dim_t e;
54 __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
55 double *EM_S, *EM_F, *DSDX;
56 index_t *row_index;
57 register dim_t q, s, r;
58 register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
59 bool_t add_EM_F, add_EM_S;
60
61 bool_t extendedA = isExpanded(A);
62 bool_t extendedB = isExpanded(B);
63 bool_t extendedC = isExpanded(C);
64 bool_t extendedD = isExpanded(D);
65 bool_t extendedX = isExpanded(X);
66 bool_t extendedY = isExpanded(Y);
67 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68 const double *S = p.shapeFns;
69 dim_t len_EM_S = p.numShapes * p.numShapes;
70 dim_t len_EM_F = p.numShapes;
71
72 #pragma omp parallel private(color,EM_S, EM_F, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q,row_index,q, s,r,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
73 {
74 EM_S = THREAD_MEMALLOC(len_EM_S, double);
75 EM_F = THREAD_MEMALLOC(len_EM_F, double);
76 row_index = THREAD_MEMALLOC(p.numShapes, index_t);
77
78 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79 {
80
81 for (color = elements->minColor; color <= elements->maxColor; color++)
82 {
83 /* open loop over all elements: */
84 #pragma omp for private(e) schedule(static)
85 for (e = 0; e < elements->numElements; e++)
86 {
87 if (elements->Color[e] == color)
88 {
89 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
90
91 A_p = getSampleDataRO(A, e);
92 B_p = getSampleDataRO(B, e);
93 C_p = getSampleDataRO(C, e);
94 D_p = getSampleDataRO(D, e);
95 X_p = getSampleDataRO(X, e);
96 Y_p = getSampleDataRO(Y, e);
97
98 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
99 for (q = 0; q < len_EM_S; ++q)
100 EM_S[q] = 0;
101 for (q = 0; q < len_EM_F; ++q)
102 EM_F[q] = 0;
103 add_EM_F = FALSE;
104 add_EM_S = FALSE;
105
106 /************************************************************************************/
107 /* process A: */
108 /************************************************************************************/
109 if (NULL != A_p)
110 {
111 add_EM_S = TRUE;
112 if (extendedA)
113 {
114 A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)]);
115 for (s = 0; s < p.numShapes; s++)
116 {
117 for (r = 0; r < p.numShapes; r++)
118 {
119 rtmp = 0;
120 for (q = 0; q < p.numQuad; q++)
121 {
122 rtmp +=
123 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
124 A_q[INDEX3(0, 0, q, DIM, DIM)] *
125 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
126 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
127 A_q[INDEX3(0, 1, q, DIM, DIM)] *
128 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
129 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
130 A_q[INDEX3(0, 2, q, DIM, DIM)] *
131 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
132 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
133 A_q[INDEX3(1, 0, q, DIM, DIM)] *
134 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
135 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
136 A_q[INDEX3(1, 1, q, DIM, DIM)] *
137 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
138 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
139 A_q[INDEX3(1, 2, q, DIM, DIM)] *
140 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
141 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
142 A_q[INDEX3(2, 0, q, DIM, DIM)] *
143 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
144 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
145 A_q[INDEX3(2, 1, q, DIM, DIM)] *
146 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
147 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
148 A_q[INDEX3(2, 2, q, DIM, DIM)] *
149 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
150 }
151 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
152 }
153 }
154 }
155 else
156 {
157 for (s = 0; s < p.numShapes; s++)
158 {
159 for (r = 0; r < p.numShapes; r++)
160 {
161 rtmp00 = 0;
162 rtmp01 = 0;
163 rtmp02 = 0;
164 rtmp10 = 0;
165 rtmp11 = 0;
166 rtmp12 = 0;
167 rtmp20 = 0;
168 rtmp21 = 0;
169 rtmp22 = 0;
170 for (q = 0; q < p.numQuad; q++)
171 {
172
173 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
174 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
175 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
176 rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
177
178 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
179 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
180 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
181 rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
182
183 rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
184 rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
185 rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
186 rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
187 }
188 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
189 rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
190 rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
191 rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
192 rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
193 rtmp22 * A_p[INDEX2(2, 2, DIM)];
194 }
195 }
196 }
197 }
198 /************************************************************************************/
199 /* process B: */
200 /************************************************************************************/
201 if (NULL != B_p)
202 {
203 add_EM_S = TRUE;
204 if (extendedB)
205 {
206 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
207 for (s = 0; s < p.numShapes; s++)
208 {
209 for (r = 0; r < p.numShapes; r++)
210 {
211 rtmp = 0;
212 for (q = 0; q < p.numQuad; q++)
213 {
214 rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
215 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
216 B_q[INDEX2(0, q, DIM)] +
217 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
218 B_q[INDEX2(1, q, DIM)] +
219 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * B_q[INDEX2(2, q, DIM)]);
220 }
221 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
222 }
223 }
224 }
225 else
226 {
227 for (s = 0; s < p.numShapes; s++)
228 {
229 for (r = 0; r < p.numShapes; r++)
230 {
231 rtmp0 = 0;
232 rtmp1 = 0;
233 rtmp2 = 0;
234 for (q = 0; q < p.numQuad; q++)
235 {
236 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
237 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
238 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
239 rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
240 }
241 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
242 rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
243 }
244 }
245 }
246 }
247 /************************************************************************************/
248 /* process C: */
249 /************************************************************************************/
250 if (NULL != C_p)
251 {
252 add_EM_S = TRUE;
253 if (extendedC)
254 {
255 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
256 for (s = 0; s < p.numShapes; s++)
257 {
258 for (r = 0; r < p.numShapes; r++)
259 {
260 rtmp = 0;
261 for (q = 0; q < p.numQuad; q++)
262 {
263 rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
264 (C_q[INDEX2(0, q, DIM)] *
265 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
266 C_q[INDEX2(1, q, DIM)] *
267 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
268 C_q[INDEX2(2, q, DIM)] * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
269 }
270 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
271 }
272 }
273 }
274 else
275 {
276 for (s = 0; s < p.numShapes; s++)
277 {
278 for (r = 0; r < p.numShapes; r++)
279 {
280 rtmp0 = 0;
281 rtmp1 = 0;
282 rtmp2 = 0;
283 for (q = 0; q < p.numQuad; q++)
284 {
285 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
286 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
287 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
288 rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
289 }
290 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
291 rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
292 }
293 }
294 }
295 }
296 /*********************************************************************************** */
297 /* process D */
298 /************************************************************************************/
299 if (NULL != D_p)
300 {
301 add_EM_S = TRUE;
302 if (extendedD)
303 {
304 D_q = &(D_p[INDEX2(0, 0, p.numQuad)]);
305 for (s = 0; s < p.numShapes; s++)
306 {
307 for (r = 0; r < p.numShapes; r++)
308 {
309 rtmp = 0;
310 for (q = 0; q < p.numQuad; q++)
311 rtmp +=
312 vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
313 S[INDEX2(r, q, p.numShapes)];
314 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
315 }
316 }
317 }
318 else
319 {
320 for (s = 0; s < p.numShapes; s++)
321 {
322 for (r = 0; r < p.numShapes; r++)
323 {
324 rtmp = 0;
325 for (q = 0; q < p.numQuad; q++)
326 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
327 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp * D_p[0];
328 }
329 }
330 }
331 }
332 /************************************************************************************/
333 /* process X: */
334 /************************************************************************************/
335 if (NULL != X_p)
336 {
337 add_EM_F = TRUE;
338 if (extendedX)
339 {
340 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
341 for (s = 0; s < p.numShapes; s++)
342 {
343 rtmp = 0;
344 for (q = 0; q < p.numQuad; q++)
345 {
346 rtmp +=
347 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
348 X_q[INDEX2(0, q, DIM)] +
349 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
350 X_q[INDEX2(1, q, DIM)] +
351 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * X_q[INDEX2(2, q, DIM)]);
352 }
353 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
354 }
355 }
356 else
357 {
358 for (s = 0; s < p.numShapes; s++)
359 {
360 rtmp0 = 0;
361 rtmp1 = 0;
362 rtmp2 = 0;
363 for (q = 0; q < p.numQuad; q++)
364 {
365 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
366 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
367 rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
368 }
369 EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
370 }
371 }
372 }
373 /************************************************************************************/
374 /* process Y: */
375 /************************************************************************************/
376 if (NULL != Y_p)
377 {
378 add_EM_F = TRUE;
379 if (extendedY)
380 {
381 Y_q = &(Y_p[INDEX2(0, 0, p.numQuad)]);
382 for (s = 0; s < p.numShapes; s++)
383 {
384 rtmp = 0;
385 for (q = 0; q < p.numQuad; q++)
386 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
387 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
388 }
389 }
390 else
391 {
392 for (s = 0; s < p.numShapes; s++)
393 {
394 rtmp = 0;
395 for (q = 0; q < p.numQuad; q++)
396 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
397 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
398 }
399 }
400 }
401 /*********************************************************************************************************************/
402 /* add the element matrices onto the matrix and right hand side */
403 /*********************************************************************************************************************/
404
405 for (q = 0; q < p.numShapes; q++)
406 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
407
408 if (add_EM_F)
409 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
410 if (add_EM_S)
411 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
412 p.numShapes, row_index, p.numComp, EM_S);
413
414 } /* end color check */
415 } /* end element loop */
416 } /* end color loop */
417
418 THREAD_MEMFREE(EM_S);
419 THREAD_MEMFREE(EM_F);
420 THREAD_MEMFREE(row_index);
421
422 } /* end of pointer check */
423 } /* end parallel region */
424 }

  ViewVC Help
Powered by ViewVC 1.1.26