/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_System2_1D.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_1D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 12950 byte(s)
like sand though the hourglass
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 /************************************************************************************/
17
18 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
19 /* the shape functions for test and solution must be identical */
20
21 /* -(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 */
22
23 /* u has p.numComp components in a 1D domain. The shape functions for test and solution must be identical */
24 /* and row_NS == row_NN */
25
26 /* Shape of the coefficients: */
27
28 /* A = p.numEqu x 1 x p.numComp x 1 */
29 /* B = 1 x numEqu x p.numComp */
30 /* C = p.numEqu x 1 x p.numComp */
31 /* D = p.numEqu x p.numComp */
32 /* X = p.numEqu x 1 */
33 /* Y = p.numEqu */
34
35 /************************************************************************************/
36
37 #include "Assemble.h"
38 #include "Util.h"
39 #ifdef _OPENMP
40 #include <omp.h>
41 #endif
42
43 /************************************************************************************/
44
45 void Dudley_Assemble_PDE_System2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46 Paso_SystemMatrix * Mat, escriptDataC * F,
47 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48 escriptDataC * X, escriptDataC * Y)
49 {
50
51 #define DIM 1
52 index_t color;
53 dim_t e;
54 __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;
55 double *EM_S, *EM_F, *DSDX;
56 index_t *row_index;
57 register dim_t q, s, r, k, m;
58 register double rtmp;
59 bool_t add_EM_F, add_EM_S;
60
61 bool_t extendedA = isExpanded(A);
62 bool_t extendedB = isExpanded(B);
63 bool_t extendedC = isExpanded(C);
64 bool_t extendedD = isExpanded(D);
65 bool_t extendedX = isExpanded(X);
66 bool_t extendedY = isExpanded(Y);
67 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68 double *S = p.row_jac->BasisFunctions->S;
69 dim_t len_EM_S = p.row_numShapes * p.row_numShapes * p.numEqu * p.numComp;
70 dim_t len_EM_F = p.row_numShapes * p.numEqu;
71
72 #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)
73 {
74 EM_S = new double[len_EM_S];
75 EM_F = new double[len_EM_F];
76 row_index = new index_t[p.row_numShapes];
77
78 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79 {
80
81 for (color = elements->minColor; color <= elements->maxColor; color++)
82 {
83 /* open loop over all elements: */
84 #pragma omp for private(e) schedule(static)
85 for (e = 0; e < elements->numElements; e++)
86 {
87 if (elements->Color[e] == color)
88 {
89
90 A_p = getSampleDataRO(A, e);
91 B_p = getSampleDataRO(B, e);
92 C_p = getSampleDataRO(C, e);
93 D_p = getSampleDataRO(D, e);
94 X_p = getSampleDataRO(X, e);
95 Y_p = getSampleDataRO(Y, e);
96 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
97 Vol = &(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal, 1)]);
98 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapes, DIM, p.numQuadTotal, 1)]);
99 for (q = 0; q < len_EM_S; ++q)
100 EM_S[q] = 0;
101 for (q = 0; q < len_EM_F; ++q)
102 EM_F[q] = 0;
103 add_EM_F = FALSE;
104 add_EM_S = FALSE;
105
106 /************************************************************************************/
107 /* process A: */
108 /************************************************************************************/
109
110 if (NULL != A_p)
111 {
112 add_EM_S = TRUE;
113 if (extendedA)
114 {
115 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
116 for (s = 0; s < p.row_numShapes; s++)
117 {
118 for (r = 0; r < p.row_numShapes; r++)
119 {
120 for (k = 0; k < p.numEqu; k++)
121 {
122 for (m = 0; m < p.numComp; m++)
123 {
124 rtmp = 0.;
125 for (q = 0; q < p.numQuadTotal; q++)
126 {
127 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
128 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] *
129 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
130 }
131 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
132 }
133 }
134 }
135 }
136 }
137 else
138 {
139 for (s = 0; s < p.row_numShapes; s++)
140 {
141 for (r = 0; r < p.row_numShapes; r++)
142 {
143 rtmp = 0;
144 for (q = 0; q < p.numQuadTotal; q++)
145 rtmp +=
146 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
147 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
148 for (k = 0; k < p.numEqu; k++)
149 {
150 for (m = 0; m < p.numComp; m++)
151 {
152 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
153 rtmp * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)];
154 }
155 }
156 }
157 }
158 }
159 }
160 /************************************************************************************/
161 /* process B: */
162 /************************************************************************************/
163
164 if (NULL != B_p)
165 {
166 add_EM_S = TRUE;
167 if (extendedB)
168 {
169 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
170 for (s = 0; s < p.row_numShapes; s++)
171 {
172 for (r = 0; r < p.row_numShapes; r++)
173 {
174 for (k = 0; k < p.numEqu; k++)
175 {
176 for (m = 0; m < p.numComp; m++)
177 {
178 rtmp = 0.;
179 for (q = 0; q < p.numQuadTotal; q++)
180 {
181 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
182 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] *
183 S[INDEX2(r, q, p.row_numShapes)];
184 }
185 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
186 }
187 }
188 }
189 }
190 }
191 else
192 {
193 for (s = 0; s < p.row_numShapes; s++)
194 {
195 for (r = 0; r < p.row_numShapes; r++)
196 {
197 rtmp = 0;
198 for (q = 0; q < p.numQuadTotal; q++)
199 rtmp +=
200 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
201 S[INDEX2(r, q, p.row_numShapes)];
202 for (k = 0; k < p.numEqu; k++)
203 {
204 for (m = 0; m < p.numComp; m++)
205 {
206 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
207 rtmp * B_p[INDEX3(k, 0, m, p.numEqu, DIM)];
208 }
209 }
210 }
211 }
212 }
213 }
214 /************************************************************************************/
215 /* process C: */
216 /************************************************************************************/
217
218 if (NULL != C_p)
219 {
220 add_EM_S = TRUE;
221 if (extendedC)
222 {
223 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
224 for (s = 0; s < p.row_numShapes; s++)
225 {
226 for (r = 0; r < p.row_numShapes; r++)
227 {
228 for (k = 0; k < p.numEqu; k++)
229 {
230 for (m = 0; m < p.numComp; m++)
231 {
232 rtmp = 0;
233 for (q = 0; q < p.numQuadTotal; q++)
234 {
235 rtmp +=
236 vol * S[INDEX2(s, q, p.row_numShapes)] *
237 C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
238 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
239 }
240 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
241 }
242 }
243 }
244 }
245 }
246 else
247 {
248 for (s = 0; s < p.row_numShapes; s++)
249 {
250 for (r = 0; r < p.row_numShapes; r++)
251 {
252 rtmp = 0;
253 for (q = 0; q < p.numQuadTotal; q++)
254 rtmp +=
255 vol * S[INDEX2(s, q, p.row_numShapes)] *
256 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
257 for (k = 0; k < p.numEqu; k++)
258 {
259 for (m = 0; m < p.numComp; m++)
260 {
261 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
262 rtmp * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)];
263 }
264 }
265 }
266 }
267 }
268 }
269 /*********************************************************************************** */
270 /* process D */
271 /************************************************************************************/
272
273 if (NULL != D_p)
274 {
275 add_EM_S = TRUE;
276 if (extendedD)
277 {
278 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
279 for (s = 0; s < p.row_numShapes; s++)
280 {
281 for (r = 0; r < p.row_numShapes; r++)
282 {
283 for (k = 0; k < p.numEqu; k++)
284 {
285 for (m = 0; m < p.numComp; m++)
286 {
287 rtmp = 0;
288 for (q = 0; q < p.numQuadTotal; q++)
289 {
290 rtmp +=
291 vol * S[INDEX2(s, q, p.row_numShapes)] *
292 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
293 S[INDEX2(r, q, p.row_numShapes)];
294 }
295 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
296
297 }
298 }
299 }
300 }
301 }
302 else
303 {
304 for (s = 0; s < p.row_numShapes; s++)
305 {
306 for (r = 0; r < p.row_numShapes; r++)
307 {
308 rtmp = 0;
309 for (q = 0; q < p.numQuadTotal; q++)
310 rtmp +=
311 vol * S[INDEX2(s, q, p.row_numShapes)] *
312 S[INDEX2(r, q, p.row_numShapes)];
313 for (k = 0; k < p.numEqu; k++)
314 {
315 for (m = 0; m < p.numComp; m++)
316 {
317 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
318 rtmp * D_p[INDEX2(k, m, p.numEqu)];
319 }
320 }
321 }
322 }
323 }
324 }
325 /************************************************************************************/
326 /* process X: */
327 /************************************************************************************/
328
329 if (NULL != X_p)
330 {
331 add_EM_F = TRUE;
332 if (extendedX)
333 {
334 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
335 for (s = 0; s < p.row_numShapes; s++)
336 {
337 for (k = 0; k < p.numEqu; k++)
338 {
339 rtmp = 0;
340 for (q = 0; q < p.numQuadTotal; q++)
341 rtmp +=
342 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
343 X_q[INDEX3(k, 0, q, p.numEqu, DIM)];
344 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
345 }
346 }
347 }
348 else
349 {
350 for (s = 0; s < p.row_numShapes; s++)
351 {
352 rtmp = 0;
353 for (q = 0; q < p.numQuadTotal; q++)
354 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];
355 for (k = 0; k < p.numEqu; k++)
356 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * X_p[INDEX2(k, 0, p.numEqu)];
357 }
358 }
359 }
360 /************************************************************************************/
361 /* process Y: */
362 /************************************************************************************/
363
364 if (NULL != Y_p)
365 {
366 add_EM_F = TRUE;
367 if (extendedY)
368 {
369 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
370 for (s = 0; s < p.row_numShapes; s++)
371 {
372 for (k = 0; k < p.numEqu; k++)
373 {
374 rtmp = 0;
375 for (q = 0; q < p.numQuadTotal; q++)
376 rtmp +=
377 vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
378 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
379 }
380 }
381 }
382 else
383 {
384 for (s = 0; s < p.row_numShapes; s++)
385 {
386 rtmp = 0;
387 for (q = 0; q < p.numQuadTotal; q++)
388 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
389 for (k = 0; k < p.numEqu; k++)
390 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
391 }
392 }
393 }
394 /*********************************************************************************************************************/
395 /* add the element matrices onto the matrix and right hand side */
396 /*********************************************************************************************************************/
397
398 for (q = 0; q < p.row_numShapes; q++)
399 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
400
401 if (add_EM_F)
402 Dudley_Util_AddScatter(p.row_numShapes, row_index, p.numEqu, EM_F, F_p,
403 p.row_DOF_UpperBound);
404 if (add_EM_S)
405 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapes, row_index, p.numEqu,
406 p.row_numShapes, row_index, p.numComp, EM_S);
407
408 } /* end color check */
409 } /* end element loop */
410 } /* end color loop */
411
412 delete[] EM_S;
413 delete[] EM_F;
414 delete[] row_index;
415
416 } /* end of pointer check */
417 } /* end parallel region */
418 }
419
420 /*
421 * $Log$
422 */

  ViewVC Help
Powered by ViewVC 1.1.26