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

Contents of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.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: 14770 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 2D 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 2 x p.numComp x 2 */
29 /* B = 2 x p.numEqu x p.numComp */
30 /* C = p.numEqu x 2 x p.numComp */
31 /* D = p.numEqu x p.numComp */
32 /* X = p.numEqu x 2 */
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_2D(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 2
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, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
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 const double *S = p.shapeFns;
69 dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
70 dim_t len_EM_F = p.numShapes * p.numEqu;
71
72 #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, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)
73 {
74
75 EM_S = new double[len_EM_S];
76 EM_F = new double[len_EM_F];
77 row_index = new index_t[p.numShapes];
78
79 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
80 {
81
82 for (color = elements->minColor; color <= elements->maxColor; color++)
83 {
84 /* open loop over all elements: */
85 #pragma omp for private(e) schedule(static)
86 for (e = 0; e < elements->numElements; e++)
87 {
88 if (elements->Color[e] == color)
89 {
90 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
91
92 A_p = getSampleDataRO(A, e);
93 B_p = getSampleDataRO(B, e);
94 C_p = getSampleDataRO(C, e);
95 D_p = getSampleDataRO(D, e);
96 X_p = getSampleDataRO(X, e);
97 Y_p = getSampleDataRO(Y, e);
98 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 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 A_p = getSampleDataRO(A, e);
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.numQuad)]);
116 for (s = 0; s < p.numShapes; s++)
117 {
118 for (r = 0; r < p.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.numQuad; q++)
126 {
127 rtmp +=
128 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
130 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
131 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
132 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
133 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
134 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
136 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
137 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
138 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
139 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
140 }
141 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
142 }
143 }
144 }
145 }
146 }
147 else
148 {
149 for (s = 0; s < p.numShapes; s++)
150 {
151 for (r = 0; r < p.numShapes; r++)
152 {
153 rtmp00 = 0;
154 rtmp01 = 0;
155 rtmp10 = 0;
156 rtmp11 = 0;
157 for (q = 0; q < p.numQuad; q++)
158 {
159 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
160 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
161 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
162 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
163 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
164 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
165 }
166 for (k = 0; k < p.numEqu; k++)
167 {
168 for (m = 0; m < p.numComp; m++)
169 {
170 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
171 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
172 + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
173 + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
174 + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
175 }
176 }
177 }
178 }
179 }
180 }
181 /************************************************************************************/
182 /* process B: */
183 /************************************************************************************/
184 B_p = getSampleDataRO(B, e);
185 if (NULL != B_p)
186 {
187 add_EM_S = TRUE;
188 if (extendedB)
189 {
190 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]);
191 for (s = 0; s < p.numShapes; s++)
192 {
193 for (r = 0; r < p.numShapes; r++)
194 {
195 for (k = 0; k < p.numEqu; k++)
196 {
197 for (m = 0; m < p.numComp; m++)
198 {
199 rtmp = 0;
200 for (q = 0; q < p.numQuad; q++)
201 {
202 rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
203 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
204 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
205 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
206 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
207 }
208 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
209 }
210 }
211 }
212 }
213 }
214 else
215 {
216 for (s = 0; s < p.numShapes; s++)
217 {
218 for (r = 0; r < p.numShapes; r++)
219 {
220 rtmp0 = 0;
221 rtmp1 = 0;
222 for (q = 0; q < p.numQuad; q++)
223 {
224 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
225 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
226 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
227 }
228 for (k = 0; k < p.numEqu; k++)
229 {
230 for (m = 0; m < p.numComp; m++)
231 {
232 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
233 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
234 rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
235 }
236 }
237 }
238 }
239 }
240 }
241 /************************************************************************************/
242 /* process C: */
243 /************************************************************************************/
244 C_p = getSampleDataRO(C, e);
245 if (NULL != C_p)
246 {
247 add_EM_S = TRUE;
248 if (extendedC)
249 {
250 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]);
251 for (s = 0; s < p.numShapes; s++)
252 {
253 for (r = 0; r < p.numShapes; r++)
254 {
255 for (k = 0; k < p.numEqu; k++)
256 {
257 for (m = 0; m < p.numComp; m++)
258 {
259 rtmp = 0;
260 for (q = 0; q < p.numQuad; q++)
261 {
262 rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
263 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
264 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
265 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
266 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
267 }
268 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
269 }
270 }
271 }
272 }
273 }
274 else
275 {
276 for (s = 0; s < p.numShapes; s++)
277 {
278 for (r = 0; r < p.numShapes; r++)
279 {
280 rtmp0 = 0;
281 rtmp1 = 0;
282 for (q = 0; q < p.numQuad; q++)
283 {
284 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
285 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
286 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
287 }
288 for (k = 0; k < p.numEqu; k++)
289 {
290 for (m = 0; m < p.numComp; m++)
291 {
292 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
293 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
294 rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
295 }
296 }
297 }
298 }
299 }
300 }
301 /*********************************************************************************** */
302 /* process D */
303 /************************************************************************************/
304 D_p = getSampleDataRO(D, e);
305 if (NULL != D_p)
306 {
307 add_EM_S = TRUE;
308 if (extendedD)
309 {
310 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
311 for (s = 0; s < p.numShapes; s++)
312 {
313 for (r = 0; r < p.numShapes; r++)
314 {
315 for (k = 0; k < p.numEqu; k++)
316 {
317 for (m = 0; m < p.numComp; m++)
318 {
319 rtmp = 0;
320 for (q = 0; q < p.numQuad; q++)
321 {
322 rtmp +=
323 vol * S[INDEX2(s, q, p.numShapes)] *
324 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
325 S[INDEX2(r, q, p.numShapes)];
326 }
327 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
328 }
329 }
330 }
331 }
332 }
333 else
334 {
335 for (s = 0; s < p.numShapes; s++)
336 {
337 for (r = 0; r < p.numShapes; r++)
338 {
339 rtmp = 0;
340 for (q = 0; q < p.numQuad; q++)
341 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
342 for (k = 0; k < p.numEqu; k++)
343 {
344 for (m = 0; m < p.numComp; m++)
345 {
346 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
347 rtmp * D_p[INDEX2(k, m, p.numEqu)];
348 }
349 }
350 }
351 }
352 }
353 }
354 /************************************************************************************/
355 /* process X: */
356 /************************************************************************************/
357 X_p = getSampleDataRO(X, e);
358 if (NULL != X_p)
359 {
360 add_EM_F = TRUE;
361 if (extendedX)
362 {
363 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
364 for (s = 0; s < p.numShapes; s++)
365 {
366 for (k = 0; k < p.numEqu; k++)
367 {
368 rtmp = 0;
369 for (q = 0; q < p.numQuad; q++)
370 {
371 rtmp +=
372 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
373 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
375 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376 }
377 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378 }
379 }
380 }
381 else
382 {
383 for (s = 0; s < p.numShapes; s++)
384 {
385 rtmp0 = 0;
386 rtmp1 = 0;
387 for (q = 0; q < p.numQuad; q++)
388 {
389 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
390 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
391 }
392 for (k = 0; k < p.numEqu; k++)
393 EM_F[INDEX2(k, s, p.numEqu)] +=
394 rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
395 }
396 }
397 }
398 /************************************************************************************/
399 /* process Y: */
400 /************************************************************************************/
401 Y_p = getSampleDataRO(Y, e);
402 if (NULL != Y_p)
403 {
404 add_EM_F = TRUE;
405 if (extendedY)
406 {
407 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
408 for (s = 0; s < p.numShapes; s++)
409 {
410 for (k = 0; k < p.numEqu; k++)
411 {
412 rtmp = 0;
413 for (q = 0; q < p.numQuad; q++)
414 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
415 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416 }
417 }
418 }
419 else
420 {
421 for (s = 0; s < p.numShapes; s++)
422 {
423 rtmp = 0;
424 for (q = 0; q < p.numQuad; q++)
425 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
426 for (k = 0; k < p.numEqu; k++)
427 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
428 }
429 }
430 }
431 /*********************************************************************************************************************/
432 /* add the element matrices onto the matrix and right hand side */
433 /*********************************************************************************************************************/
434 for (q = 0; q < p.numShapes; q++)
435 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
436
437 if (add_EM_F)
438 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
439 if (add_EM_S)
440 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
441 p.numShapes, row_index, p.numComp, EM_S);
442
443 } /* end color check */
444 } /* end element loop */
445 } /* end color loop */
446
447 delete[] EM_S;
448 delete[] EM_F;
449 delete[] row_index;
450
451 } /* end of pointer check */
452 } /* end parallel region */
453 }
454
455 /*
456 * $Log$
457 */

  ViewVC Help
Powered by ViewVC 1.1.26