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_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
21
/* 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: */
33
/**************************************************************/
41
/**************************************************************/
43
void Dudley_Assemble_PDE_Single2_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;
56
register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
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;
68
dim_t len_EM_F = p.numShapes;
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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,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);
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)
99
for (q = 0; q < len_EM_F; ++q)
104
/**************************************************************/
106
/**************************************************************/
112
A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)]);
113
for (s = 0; s < p.numShapes; s++)
115
for (r = 0; r < p.numShapes; r++)
118
for (q = 0; q < p.numQuad; q++)
121
vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
122
A_q[INDEX3(0, 0, q, DIM, DIM)] *
123
DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
124
DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
125
A_q[INDEX3(0, 1, q, DIM, DIM)] *
126
DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
127
DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
128
A_q[INDEX3(0, 2, q, DIM, DIM)] *
129
DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
130
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
131
A_q[INDEX3(1, 0, q, DIM, DIM)] *
132
DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
133
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
134
A_q[INDEX3(1, 1, q, DIM, DIM)] *
135
DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
136
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
137
A_q[INDEX3(1, 2, q, DIM, DIM)] *
138
DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
139
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
140
A_q[INDEX3(2, 0, q, DIM, DIM)] *
141
DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
142
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
143
A_q[INDEX3(2, 1, q, DIM, DIM)] *
144
DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
145
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
146
A_q[INDEX3(2, 2, q, DIM, DIM)] *
147
DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
149
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
155
for (s = 0; s < p.numShapes; s++)
157
for (r = 0; r < p.numShapes; r++)
168
for (q = 0; q < p.numQuad; q++)
171
rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
172
rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
173
rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
174
rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
176
rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
177
rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
178
rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
179
rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
181
rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
182
rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
183
rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
184
rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
186
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
187
rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
188
rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
189
rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
190
rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
191
rtmp22 * A_p[INDEX2(2, 2, DIM)];
196
/**************************************************************/
198
/**************************************************************/
204
B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
205
for (s = 0; s < p.numShapes; s++)
207
for (r = 0; r < p.numShapes; r++)
210
for (q = 0; q < p.numQuad; q++)
212
rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
213
(DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
214
B_q[INDEX2(0, q, DIM)] +
215
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
216
B_q[INDEX2(1, q, DIM)] +
217
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * B_q[INDEX2(2, q, DIM)]);
219
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
225
for (s = 0; s < p.numShapes; s++)
227
for (r = 0; r < p.numShapes; r++)
232
for (q = 0; q < p.numQuad; q++)
234
rtmp = vol * S[INDEX2(r, q, p.numShapes)];
235
rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
236
rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
237
rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
239
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
240
rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
245
/**************************************************************/
247
/**************************************************************/
253
C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
254
for (s = 0; s < p.numShapes; s++)
256
for (r = 0; r < p.numShapes; r++)
259
for (q = 0; q < p.numQuad; q++)
261
rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
262
(C_q[INDEX2(0, q, DIM)] *
263
DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
264
C_q[INDEX2(1, q, DIM)] *
265
DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
266
C_q[INDEX2(2, q, DIM)] * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
268
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
274
for (s = 0; s < p.numShapes; s++)
276
for (r = 0; r < p.numShapes; r++)
281
for (q = 0; q < p.numQuad; q++)
283
rtmp = vol * S[INDEX2(s, q, p.numShapes)];
284
rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
285
rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
286
rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
288
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
289
rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
294
/************************************************************* */
296
/**************************************************************/
302
D_q = &(D_p[INDEX2(0, 0, p.numQuad)]);
303
for (s = 0; s < p.numShapes; s++)
305
for (r = 0; r < p.numShapes; r++)
308
for (q = 0; q < p.numQuad; q++)
310
vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
311
S[INDEX2(r, q, p.numShapes)];
312
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
318
for (s = 0; s < p.numShapes; s++)
320
for (r = 0; r < p.numShapes; r++)
323
for (q = 0; q < p.numQuad; q++)
324
rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
325
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp * D_p[0];
330
/**************************************************************/
332
/**************************************************************/
338
X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
339
for (s = 0; s < p.numShapes; s++)
342
for (q = 0; q < p.numQuad; q++)
345
vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
346
X_q[INDEX2(0, q, DIM)] +
347
DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
348
X_q[INDEX2(1, q, DIM)] +
349
DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * X_q[INDEX2(2, q, DIM)]);
351
EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
356
for (s = 0; s < p.numShapes; s++)
361
for (q = 0; q < p.numQuad; q++)
363
rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
364
rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
365
rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
367
EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
371
/**************************************************************/
373
/**************************************************************/
379
Y_q = &(Y_p[INDEX2(0, 0, p.numQuad)]);
380
for (s = 0; s < p.numShapes; s++)
383
for (q = 0; q < p.numQuad; q++)
384
rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
385
EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
390
for (s = 0; s < p.numShapes; s++)
393
for (q = 0; q < p.numQuad; q++)
394
rtmp += vol * S[INDEX2(s, q, p.numShapes)];
395
EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
399
/***********************************************************************************************/
400
/* add the element matrices onto the matrix and right hand side */
401
/***********************************************************************************************/
403
for (q = 0; q < p.numShapes; q++)
404
row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
407
Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
409
Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
410
p.numShapes, row_index, p.numComp, EM_S);
412
} /* end color check */
413
} /* end element loop */
414
} /* end color loop */
416
THREAD_MEMFREE(EM_S);
417
THREAD_MEMFREE(EM_F);
418
THREAD_MEMFREE(row_index);
420
} /* end of pointer check */
421
} /* end parallel region */