2
/*******************************************************
4
* Copyright (c) 2003-2010 by University of Queensland
5
* Earth Systems Science Computational Center (ESSCC)
6
* http://www.uq.edu.au/esscc
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
12
*******************************************************/
14
/**************************************************************/
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 */
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 */
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 */
24
/* Shape of the coefficients: */
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 */
33
/**************************************************************/
41
/**************************************************************/
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)
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;
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;
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;
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)
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);
76
if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79
for (color = elements->minColor; color <= elements->maxColor; color++)
81
/* open loop over all elements: */
82
#pragma omp for private(e) schedule(static)
83
for (e = 0; e < elements->numElements; e++)
85
if (elements->Color[e] == color)
87
double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
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);
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)
100
for (q = 0; q < len_EM_F; ++q)
105
/**************************************************************/
107
/**************************************************************/
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++)
116
for (r = 0; r < p.numShapes; r++)
118
for (k = 0; k < p.numEqu; k++)
120
for (m = 0; m < p.numComp; m++)
123
for (q = 0; q < p.numQuad; q++)
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)]);
155
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
163
for (s = 0; s < p.numShapes; s++)
165
for (r = 0; r < p.numShapes; r++)
176
for (q = 0; q < p.numQuad; q++)
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)];
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)];
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)];
193
for (k = 0; k < p.numEqu; k++)
195
for (m = 0; m < p.numComp; m++)
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)];
213
/**************************************************************/
215
/**************************************************************/
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++)
224
for (r = 0; r < p.numShapes; r++)
226
for (k = 0; k < p.numEqu; k++)
228
for (m = 0; m < p.numComp; m++)
231
for (q = 0; q < p.numQuad; q++)
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)]);
242
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
250
for (s = 0; s < p.numShapes; s++)
252
for (r = 0; r < p.numShapes; r++)
257
for (q = 0; q < p.numQuad; q++)
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)];
264
for (k = 0; k < p.numEqu; k++)
266
for (m = 0; m < p.numComp; m++)
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)];
278
/**************************************************************/
280
/**************************************************************/
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++)
289
for (r = 0; r < p.numShapes; r++)
291
for (k = 0; k < p.numEqu; k++)
293
for (m = 0; m < p.numComp; m++)
296
for (q = 0; q < p.numQuad; q++)
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)]);
307
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
315
for (s = 0; s < p.numShapes; s++)
317
for (r = 0; r < p.numShapes; r++)
322
for (q = 0; q < p.numQuad; q++)
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)];
329
for (k = 0; k < p.numEqu; k++)
331
for (m = 0; m < p.numComp; m++)
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)];
343
/************************************************************* */
345
/**************************************************************/
351
D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
352
for (s = 0; s < p.numShapes; s++)
354
for (r = 0; r < p.numShapes; r++)
356
for (k = 0; k < p.numEqu; k++)
358
for (m = 0; m < p.numComp; m++)
361
for (q = 0; q < p.numQuad; q++)
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)];
368
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
376
for (s = 0; s < p.numShapes; s++)
378
for (r = 0; r < p.numShapes; r++)
381
for (q = 0; q < p.numQuad; q++)
382
rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
383
for (k = 0; k < p.numEqu; k++)
385
for (m = 0; m < p.numComp; m++)
387
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
388
rtmp * D_p[INDEX2(k, m, p.numEqu)];
395
/**************************************************************/
397
/**************************************************************/
403
X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
404
for (s = 0; s < p.numShapes; s++)
406
for (k = 0; k < p.numEqu; k++)
409
for (q = 0; q < p.numQuad; q++)
412
vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
413
X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
414
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
415
X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
416
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
417
X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
419
EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
425
for (s = 0; s < p.numShapes; s++)
430
for (q = 0; q < p.numQuad; q++)
432
rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
433
rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
434
rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
436
for (k = 0; k < p.numEqu; k++)
438
EM_F[INDEX2(k, s, p.numEqu)] += rtmp0 * X_p[INDEX2(k, 0, p.numEqu)]
439
+ rtmp1 * X_p[INDEX2(k, 1, p.numEqu)] + rtmp2 * X_p[INDEX2(k, 2, p.numEqu)];
444
/**************************************************************/
446
/**************************************************************/
452
Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
453
for (s = 0; s < p.numShapes; s++)
455
for (k = 0; k < p.numEqu; k++)
458
for (q = 0; q < p.numQuad; q++)
459
rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
460
EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
466
for (s = 0; s < p.numShapes; s++)
469
for (q = 0; q < p.numQuad; q++)
470
rtmp += vol * S[INDEX2(s, q, p.numShapes)];
471
for (k = 0; k < p.numEqu; k++)
472
EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
477
/***********************************************************************************************/
478
/* add the element matrices onto the matrix and right hand side */
479
/***********************************************************************************************/
480
for (q = 0; q < p.numShapes; q++)
481
row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
484
Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
486
Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
487
p.numShapes, row_index, p.numComp, EM_S);
488
} /* end color check */
489
} /* end element loop */
490
} /* end color loop */
492
THREAD_MEMFREE(EM_S);
493
THREAD_MEMFREE(EM_F);
494
THREAD_MEMFREE(row_index);
496
} /* end of pointer check */
497
} /* end parallel region */