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 1D 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_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;
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_numShapesTotal * p.row_numShapesTotal;
68
dim_t len_EM_F = p.row_numShapesTotal;
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,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_numShapesTotal, 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
C_p = getSampleDataRO(C, e);
90
B_p = getSampleDataRO(B, e);
91
D_p = getSampleDataRO(D, e);
92
X_p = getSampleDataRO(X, e);
93
Y_p = getSampleDataRO(Y, e);
95
double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
96
DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
97
for (q = 0; q < len_EM_S; ++q)
99
for (q = 0; q < len_EM_F; ++q)
103
/**************************************************************/
105
/**************************************************************/
111
A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
112
for (s = 0; s < p.row_numShapes; s++)
114
for (r = 0; r < p.row_numShapes; r++)
117
for (q = 0; q < p.numQuadTotal; q++)
120
vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
121
A_q[INDEX3(0, 0, q, DIM, DIM)] *
122
DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
124
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
130
for (s = 0; s < p.row_numShapes; s++)
132
for (r = 0; r < p.row_numShapes; r++)
135
for (q = 0; q < p.numQuadTotal; q++)
137
vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
138
DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
139
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
140
rtmp * A_p[INDEX2(0, 0, DIM)];
145
/**************************************************************/
147
/**************************************************************/
153
B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
154
for (s = 0; s < p.row_numShapes; s++)
156
for (r = 0; r < p.row_numShapes; r++)
159
for (q = 0; q < p.numQuadTotal; q++)
162
vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
163
B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
165
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
171
for (s = 0; s < p.row_numShapes; s++)
173
for (r = 0; r < p.row_numShapes; r++)
176
for (q = 0; q < p.numQuadTotal; q++)
178
vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
179
S[INDEX2(r, q, p.row_numShapes)];
180
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
186
/**************************************************************/
188
/**************************************************************/
194
C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
195
for (s = 0; s < p.row_numShapes; s++)
197
for (r = 0; r < p.row_numShapes; r++)
200
for (q = 0; q < p.numQuadTotal; q++)
203
vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
204
DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
206
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
212
for (s = 0; s < p.row_numShapes; s++)
214
for (r = 0; r < p.row_numShapes; r++)
217
for (q = 0; q < p.numQuadTotal; q++)
219
vol * S[INDEX2(s, q, p.row_numShapes)] *
220
DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
221
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
227
/************************************************************* */
229
/**************************************************************/
235
D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
236
for (s = 0; s < p.row_numShapes; s++)
238
for (r = 0; r < p.row_numShapes; r++)
241
for (q = 0; q < p.numQuadTotal; q++)
244
vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
245
S[INDEX2(r, q, p.row_numShapes)];
247
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
253
for (s = 0; s < p.row_numShapes; s++)
255
for (r = 0; r < p.row_numShapes; r++)
258
for (q = 0; q < p.numQuadTotal; q++)
260
vol * S[INDEX2(s, q, p.row_numShapes)] *
261
S[INDEX2(r, q, p.row_numShapes)];
262
EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
268
/**************************************************************/
270
/**************************************************************/
276
X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
277
for (s = 0; s < p.row_numShapes; s++)
280
for (q = 0; q < p.numQuadTotal; q++)
282
vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
283
X_q[INDEX2(0, q, DIM)];
284
EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
289
for (s = 0; s < p.row_numShapes; s++)
292
for (q = 0; q < p.numQuadTotal; q++)
293
rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
294
EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
298
/**************************************************************/
300
/**************************************************************/
306
Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
307
for (s = 0; s < p.row_numShapes; s++)
310
for (q = 0; q < p.numQuadTotal; q++)
311
rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
312
EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
317
for (s = 0; s < p.row_numShapes; s++)
320
for (q = 0; q < p.numQuadTotal; q++)
321
rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
322
EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
326
/***********************************************************************************************/
327
/* add the element matrices onto the matrix and right hand side */
328
/***********************************************************************************************/
329
for (q = 0; q < p.row_numShapesTotal; q++)
330
row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
333
Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
334
p.row_DOF_UpperBound);
336
Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
337
p.row_numShapesTotal, row_index, p.numComp, EM_S);
339
} /* end color check */
340
} /* end element loop */
341
} /* end color loop */
343
THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
344
THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
345
THREAD_MEMFREE(row_index);
347
} /* end of pointer check */
348
} /* end parallel region */