/[escript]/branches/domexper/dudley/src/Assemble_PDE_System2_3D.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_PDE_System2_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3205 - (show annotations)
Fri Sep 24 00:30:43 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 16826 byte(s)
Removing -total from field names

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

  ViewVC Help
Powered by ViewVC 1.1.26