/[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 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 16757 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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(Dudley_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 const double *S = p.shapeFns;
67 dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
68 dim_t len_EM_F = p.numShapes * p.numEqu;
69
70 #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,k,m,rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,add_EM_F, add_EM_S)
71 {
72 EM_S = THREAD_MEMALLOC(len_EM_S, double);
73 EM_F = THREAD_MEMALLOC(len_EM_F, double);
74 row_index = THREAD_MEMALLOC(p.numShapes, index_t);
75
76 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
77 {
78
79 for (color = elements->minColor; color <= elements->maxColor; color++)
80 {
81 /* open loop over all elements: */
82 #pragma omp for private(e) schedule(static)
83 for (e = 0; e < elements->numElements; e++)
84 {
85 if (elements->Color[e] == color)
86 {
87
88 A_p = getSampleDataRO(A, e);
89 B_p = getSampleDataRO(B, e);
90 C_p = getSampleDataRO(C, e);
91 D_p = getSampleDataRO(D, e);
92 X_p = getSampleDataRO(X, e);
93 Y_p = getSampleDataRO(Y, e);
94 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
95
96 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
97 for (q = 0; q < len_EM_S; ++q)
98 EM_S[q] = 0;
99 for (q = 0; q < len_EM_F; ++q)
100 EM_F[q] = 0;
101 add_EM_F = FALSE;
102 add_EM_S = FALSE;
103
104 /**************************************************************/
105 /* process A: */
106 /**************************************************************/
107 if (NULL != A_p)
108 {
109 add_EM_S = TRUE;
110 if (extendedA)
111 {
112 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)]);
113 for (s = 0; s < p.numShapes; s++)
114 {
115 for (r = 0; r < p.numShapes; r++)
116 {
117 for (k = 0; k < p.numEqu; k++)
118 {
119 for (m = 0; m < p.numComp; m++)
120 {
121 rtmp = 0;
122 for (q = 0; q < p.numQuad; q++)
123 {
124 rtmp +=
125 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
126 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
127 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
128 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
130 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
131 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
132 A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
133 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
134 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
136 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
137 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
138 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
139 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
140 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
141 A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
142 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
143 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
144 A_q[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
145 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
146 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
147 A_q[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
148 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
149 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
150 A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
151 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
152
153 }
154 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
155 }
156 }
157 }
158 }
159 }
160 else
161 {
162 for (s = 0; s < p.numShapes; s++)
163 {
164 for (r = 0; r < p.numShapes; r++)
165 {
166 rtmp00 = 0;
167 rtmp01 = 0;
168 rtmp02 = 0;
169 rtmp10 = 0;
170 rtmp11 = 0;
171 rtmp12 = 0;
172 rtmp20 = 0;
173 rtmp21 = 0;
174 rtmp22 = 0;
175 for (q = 0; q < p.numQuad; q++)
176 {
177 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
178 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
179 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
180 rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
181
182 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
183 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
184 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
185 rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
186
187 rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
188 rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
189 rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
190 rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
191 }
192 for (k = 0; k < p.numEqu; k++)
193 {
194 for (m = 0; m < p.numComp; m++)
195 {
196 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
197 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +
198 rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +
199 rtmp02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +
200 rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)] +
201 rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)] +
202 rtmp12 * A_p[INDEX4(k, 1, m, 2, p.numEqu, DIM, p.numComp)] +
203 rtmp20 * A_p[INDEX4(k, 2, m, 0, p.numEqu, DIM, p.numComp)] +
204 rtmp21 * A_p[INDEX4(k, 2, m, 1, p.numEqu, DIM, p.numComp)] +
205 rtmp22 * A_p[INDEX4(k, 2, m, 2, p.numEqu, DIM, p.numComp)];
206 }
207 }
208 }
209 }
210 }
211 }
212 /**************************************************************/
213 /* process B: */
214 /**************************************************************/
215 if (NULL != B_p)
216 {
217 add_EM_S = TRUE;
218 if (extendedB)
219 {
220 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]);
221 for (s = 0; s < p.numShapes; s++)
222 {
223 for (r = 0; r < p.numShapes; r++)
224 {
225 for (k = 0; k < p.numEqu; k++)
226 {
227 for (m = 0; m < p.numComp; m++)
228 {
229 rtmp = 0;
230 for (q = 0; q < p.numQuad; q++)
231 {
232 rtmp +=
233 vol * S[INDEX2(r, q, p.numShapes)] *
234 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
235 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
236 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
237 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +
238 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
239 B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);
240 }
241 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
242 }
243 }
244 }
245 }
246 }
247 else
248 {
249 for (s = 0; s < p.numShapes; s++)
250 {
251 for (r = 0; r < p.numShapes; r++)
252 {
253 rtmp0 = 0;
254 rtmp1 = 0;
255 rtmp2 = 0;
256 for (q = 0; q < p.numQuad; q++)
257 {
258 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
259 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
260 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
261 rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
262 }
263 for (k = 0; k < p.numEqu; k++)
264 {
265 for (m = 0; m < p.numComp; m++)
266 {
267 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
268 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
269 rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +
270 rtmp2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];
271 }
272 }
273 }
274 }
275 }
276 }
277 /**************************************************************/
278 /* process C: */
279 /**************************************************************/
280 if (NULL != C_p)
281 {
282 add_EM_S = TRUE;
283 if (extendedC)
284 {
285 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]);
286 for (s = 0; s < p.numShapes; s++)
287 {
288 for (r = 0; r < p.numShapes; r++)
289 {
290 for (k = 0; k < p.numEqu; k++)
291 {
292 for (m = 0; m < p.numComp; m++)
293 {
294 rtmp = 0;
295 for (q = 0; q < p.numQuad; q++)
296 {
297 rtmp +=
298 vol * S[INDEX2(s, q, p.numShapes)] *
299 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
300 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
301 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
302 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
303 C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *
304 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
305 }
306 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
307 }
308 }
309 }
310 }
311 }
312 else
313 {
314 for (s = 0; s < p.numShapes; s++)
315 {
316 for (r = 0; r < p.numShapes; r++)
317 {
318 rtmp0 = 0;
319 rtmp1 = 0;
320 rtmp2 = 0;
321 for (q = 0; q < p.numQuad; q++)
322 {
323 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
324 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
325 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
326 rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
327 }
328 for (k = 0; k < p.numEqu; k++)
329 {
330 for (m = 0; m < p.numComp; m++)
331 {
332 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
333 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
334 rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +
335 rtmp2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];
336 }
337 }
338 }
339 }
340 }
341 }
342 /************************************************************* */
343 /* process D */
344 /**************************************************************/
345 if (NULL != D_p)
346 {
347 add_EM_S = TRUE;
348 if (extendedD)
349 {
350 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
351 for (s = 0; s < p.numShapes; s++)
352 {
353 for (r = 0; r < p.numShapes; r++)
354 {
355 for (k = 0; k < p.numEqu; k++)
356 {
357 for (m = 0; m < p.numComp; m++)
358 {
359 rtmp = 0;
360 for (q = 0; q < p.numQuad; q++)
361 {
362 rtmp +=
363 vol * S[INDEX2(s, q, p.numShapes)] *
364 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
365 S[INDEX2(r, q, p.numShapes)];
366 }
367 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
368 }
369 }
370 }
371 }
372 }
373 else
374 {
375 for (s = 0; s < p.numShapes; s++)
376 {
377 for (r = 0; r < p.numShapes; r++)
378 {
379 rtmp = 0;
380 for (q = 0; q < p.numQuad; q++)
381 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
382 for (k = 0; k < p.numEqu; k++)
383 {
384 for (m = 0; m < p.numComp; m++)
385 {
386 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
387 rtmp * D_p[INDEX2(k, m, p.numEqu)];
388 }
389 }
390 }
391 }
392 }
393 }
394 /**************************************************************/
395 /* process X: */
396 /**************************************************************/
397 if (NULL != X_p)
398 {
399 add_EM_F = TRUE;
400 if (extendedX)
401 {
402 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
403 for (s = 0; s < p.numShapes; s++)
404 {
405 for (k = 0; k < p.numEqu; k++)
406 {
407 rtmp = 0;
408 for (q = 0; q < p.numQuad; q++)
409 {
410 rtmp +=
411 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
412 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
413 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
414 X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
415 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
416 X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
417 }
418 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
419 }
420 }
421 }
422 else
423 {
424 for (s = 0; s < p.numShapes; s++)
425 {
426 rtmp0 = 0;
427 rtmp1 = 0;
428 rtmp2 = 0;
429 for (q = 0; q < p.numQuad; q++)
430 {
431 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
432 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
433 rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
434 }
435 for (k = 0; k < p.numEqu; k++)
436 {
437 EM_F[INDEX2(k, s, p.numEqu)] += rtmp0 * X_p[INDEX2(k, 0, p.numEqu)]
438 + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)] + rtmp2 * X_p[INDEX2(k, 2, p.numEqu)];
439 }
440 }
441 }
442 }
443 /**************************************************************/
444 /* process Y: */
445 /**************************************************************/
446 if (NULL != Y_p)
447 {
448 add_EM_F = TRUE;
449 if (extendedY)
450 {
451 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
452 for (s = 0; s < p.numShapes; s++)
453 {
454 for (k = 0; k < p.numEqu; k++)
455 {
456 rtmp = 0.;
457 for (q = 0; q < p.numQuad; q++)
458 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
459 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
460 }
461 }
462 }
463 else
464 {
465 for (s = 0; s < p.numShapes; s++)
466 {
467 rtmp = 0;
468 for (q = 0; q < p.numQuad; q++)
469 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
470 for (k = 0; k < p.numEqu; k++)
471 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
472 }
473 }
474 }
475
476 /***********************************************************************************************/
477 /* add the element matrices onto the matrix and right hand side */
478 /***********************************************************************************************/
479 for (q = 0; q < p.numShapes; q++)
480 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
481
482 if (add_EM_F)
483 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
484 if (add_EM_S)
485 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
486 p.numShapes, row_index, p.numComp, EM_S);
487 } /* end color check */
488 } /* end element loop */
489 } /* end color loop */
490
491 THREAD_MEMFREE(EM_S);
492 THREAD_MEMFREE(EM_F);
493 THREAD_MEMFREE(row_index);
494
495 } /* end of pointer check */
496 } /* end parallel region */
497 }

  ViewVC Help
Powered by ViewVC 1.1.26