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 1D 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 1 x p.numComp x 1 */
27
/* B = 1 x numEqu x p.numComp */
28
/* C = p.numEqu x 1 x p.numComp */
29
/* D = p.numEqu x p.numComp */
30
/* X = p.numEqu x 1 */
33
/**************************************************************/
41
/**************************************************************/
43
void Dudley_Assemble_PDE_System2_1D(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;
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
double *S = p.row_jac->BasisFunctions->S;
67
dim_t len_EM_S = p.row_numShapes * p.row_numShapes * p.numEqu * p.numComp;
68
dim_t len_EM_F = p.row_numShapes * p.numEqu;
70
#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,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.row_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)
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
Vol = &(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal, 1)]);
96
DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapes, DIM, p.numQuadTotal, 1)]);
97
for (q = 0; q < len_EM_S; ++q)
99
for (q = 0; q < len_EM_F; ++q)
104
/**************************************************************/
106
/**************************************************************/
113
A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
114
for (s = 0; s < p.row_numShapes; s++)
116
for (r = 0; r < p.row_numShapes; r++)
118
for (k = 0; k < p.numEqu; k++)
120
for (m = 0; m < p.numComp; m++)
123
for (q = 0; q < p.numQuadTotal; q++)
125
rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
126
A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] *
127
DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
129
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
137
for (s = 0; s < p.row_numShapes; s++)
139
for (r = 0; r < p.row_numShapes; r++)
142
for (q = 0; q < p.numQuadTotal; q++)
144
vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
145
DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
146
for (k = 0; k < p.numEqu; k++)
148
for (m = 0; m < p.numComp; m++)
150
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
151
rtmp * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)];
158
/**************************************************************/
160
/**************************************************************/
167
B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
168
for (s = 0; s < p.row_numShapes; s++)
170
for (r = 0; r < p.row_numShapes; r++)
172
for (k = 0; k < p.numEqu; k++)
174
for (m = 0; m < p.numComp; m++)
177
for (q = 0; q < p.numQuadTotal; q++)
179
rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
180
B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] *
181
S[INDEX2(r, q, p.row_numShapes)];
183
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
191
for (s = 0; s < p.row_numShapes; s++)
193
for (r = 0; r < p.row_numShapes; r++)
196
for (q = 0; q < p.numQuadTotal; q++)
198
vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
199
S[INDEX2(r, q, p.row_numShapes)];
200
for (k = 0; k < p.numEqu; k++)
202
for (m = 0; m < p.numComp; m++)
204
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
205
rtmp * B_p[INDEX3(k, 0, m, p.numEqu, DIM)];
212
/**************************************************************/
214
/**************************************************************/
221
C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
222
for (s = 0; s < p.row_numShapes; s++)
224
for (r = 0; r < p.row_numShapes; r++)
226
for (k = 0; k < p.numEqu; k++)
228
for (m = 0; m < p.numComp; m++)
231
for (q = 0; q < p.numQuadTotal; q++)
234
vol * S[INDEX2(s, q, p.row_numShapes)] *
235
C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
236
DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
238
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
246
for (s = 0; s < p.row_numShapes; s++)
248
for (r = 0; r < p.row_numShapes; r++)
251
for (q = 0; q < p.numQuadTotal; q++)
253
vol * S[INDEX2(s, q, p.row_numShapes)] *
254
DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
255
for (k = 0; k < p.numEqu; k++)
257
for (m = 0; m < p.numComp; m++)
259
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
260
rtmp * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)];
267
/************************************************************* */
269
/**************************************************************/
276
D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
277
for (s = 0; s < p.row_numShapes; s++)
279
for (r = 0; r < p.row_numShapes; r++)
281
for (k = 0; k < p.numEqu; k++)
283
for (m = 0; m < p.numComp; m++)
286
for (q = 0; q < p.numQuadTotal; q++)
289
vol * S[INDEX2(s, q, p.row_numShapes)] *
290
D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
291
S[INDEX2(r, q, p.row_numShapes)];
293
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
302
for (s = 0; s < p.row_numShapes; s++)
304
for (r = 0; r < p.row_numShapes; r++)
307
for (q = 0; q < p.numQuadTotal; q++)
309
vol * S[INDEX2(s, q, p.row_numShapes)] *
310
S[INDEX2(r, q, p.row_numShapes)];
311
for (k = 0; k < p.numEqu; k++)
313
for (m = 0; m < p.numComp; m++)
315
EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
316
rtmp * D_p[INDEX2(k, m, p.numEqu)];
323
/**************************************************************/
325
/**************************************************************/
332
X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
333
for (s = 0; s < p.row_numShapes; s++)
335
for (k = 0; k < p.numEqu; k++)
338
for (q = 0; q < p.numQuadTotal; q++)
340
vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
341
X_q[INDEX3(k, 0, q, p.numEqu, DIM)];
342
EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
348
for (s = 0; s < p.row_numShapes; s++)
351
for (q = 0; q < p.numQuadTotal; q++)
352
rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];
353
for (k = 0; k < p.numEqu; k++)
354
EM_F[INDEX2(k, s, p.numEqu)] += rtmp * X_p[INDEX2(k, 0, p.numEqu)];
358
/**************************************************************/
360
/**************************************************************/
367
Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
368
for (s = 0; s < p.row_numShapes; s++)
370
for (k = 0; k < p.numEqu; k++)
373
for (q = 0; q < p.numQuadTotal; q++)
375
vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
376
EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
382
for (s = 0; s < p.row_numShapes; s++)
385
for (q = 0; q < p.numQuadTotal; q++)
386
rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
387
for (k = 0; k < p.numEqu; k++)
388
EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
392
/***********************************************************************************************/
393
/* add the element matrices onto the matrix and right hand side */
394
/***********************************************************************************************/
396
for (q = 0; q < p.row_numShapes; q++)
397
row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
400
Dudley_Util_AddScatter(p.row_numShapes, row_index, p.numEqu, EM_F, F_p,
401
p.row_DOF_UpperBound);
403
Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapes, row_index, p.numEqu,
404
p.row_numShapes, row_index, p.numComp, EM_S);
406
} /* end color check */
407
} /* end element loop */
408
} /* end color loop */
410
THREAD_MEMFREE(EM_S);
411
THREAD_MEMFREE(EM_F);
412
THREAD_MEMFREE(row_index);
414
} /* end of pointer check */
415
} /* end parallel region */