/[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 3187 - (show annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 8 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_System2_2D.c
File MIME type: text/plain
File size: 14971 byte(s)
Enforcing indent -linux -i4 -bl -bli0 -l120

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
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
11 *
12 *******************************************************/
13
14 /**************************************************************/
15
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 */
18
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 */
20
21 /* u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical */
22 /* and row_NS == row_NN */
23
24 /* Shape of the coefficients: */
25
26 /* A = p.numEqu x 2 x p.numComp x 2 */
27 /* B = 2 x p.numEqu x p.numComp */
28 /* C = p.numEqu x 2 x p.numComp */
29 /* D = p.numEqu x p.numComp */
30 /* X = p.numEqu x 2 */
31 /* Y = p.numEqu */
32
33 /**************************************************************/
34
35 #include "Assemble.h"
36 #include "Util.h"
37 #ifdef _OPENMP
38 #include <omp.h>
39 #endif
40
41 /**************************************************************/
42
43 void Dudley_Assemble_PDE_System2_2D(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)
47 {
48
49 #define DIM 2
50 index_t color;
51 dim_t e;
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;
54 index_t *row_index;
55 register dim_t q, s, r, k, m;
56 register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
57 bool_t add_EM_F, add_EM_S;
58
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.col_numShapesTotal * p.numEqu * p.numComp;
68 dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;
69
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, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)
71 {
72
73 EM_S = THREAD_MEMALLOC(len_EM_S, double);
74 EM_F = THREAD_MEMALLOC(len_EM_F, double);
75 row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
76
77 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
78 {
79
80 for (color = elements->minColor; color <= elements->maxColor; color++)
81 {
82 /* open loop over all elements: */
83 #pragma omp for private(e) schedule(static)
84 for (e = 0; e < elements->numElements; e++)
85 {
86 if (elements->Color[e] == color)
87 {
88
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);
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)
98 EM_S[q] = 0;
99 for (q = 0; q < len_EM_F; ++q)
100 EM_F[q] = 0;
101 add_EM_F = FALSE;
102 add_EM_S = FALSE;
103
104 /**************************************************************/
105 /* process A: */
106 /**************************************************************/
107 A_p = getSampleDataRO(A, e);
108 if (NULL != A_p)
109 {
110 add_EM_S = TRUE;
111 if (extendedA)
112 {
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++)
115 {
116 for (r = 0; r < p.col_numShapes; r++)
117 {
118 for (k = 0; k < p.numEqu; k++)
119 {
120 for (m = 0; m < p.numComp; m++)
121 {
122 rtmp = 0;
123 for (q = 0; q < p.numQuadTotal; q++)
124 {
125 rtmp +=
126 vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
127 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
128 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
129 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
130 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
131 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
132 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
133 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
134 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
135 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
136 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
137 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
138 }
139 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
140 rtmp;
141 }
142 }
143 }
144 }
145 } else
146 {
147 for (s = 0; s < p.row_numShapes; s++)
148 {
149 for (r = 0; r < p.col_numShapes; r++)
150 {
151 rtmp00 = 0;
152 rtmp01 = 0;
153 rtmp10 = 0;
154 rtmp11 = 0;
155 for (q = 0; q < p.numQuadTotal; q++)
156 {
157 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
158 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
159 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
160 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
161 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
162 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
163 }
164 for (k = 0; k < p.numEqu; k++)
165 {
166 for (m = 0; m < p.numComp; m++)
167 {
168 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
169 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
170 + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
171 + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
172 + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
173 }
174 }
175 }
176 }
177 }
178 }
179 /**************************************************************/
180 /* process B: */
181 /**************************************************************/
182 B_p = getSampleDataRO(B, e);
183 if (NULL != B_p)
184 {
185 add_EM_S = TRUE;
186 if (extendedB)
187 {
188 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
189 for (s = 0; s < p.row_numShapes; s++)
190 {
191 for (r = 0; r < p.col_numShapes; r++)
192 {
193 for (k = 0; k < p.numEqu; k++)
194 {
195 for (m = 0; m < p.numComp; m++)
196 {
197 rtmp = 0;
198 for (q = 0; q < p.numQuadTotal; q++)
199 {
200 rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
201 (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
202 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
203 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
204 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
205 }
206 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
207 rtmp;
208 }
209 }
210 }
211 }
212 } else
213 {
214 for (s = 0; s < p.row_numShapes; s++)
215 {
216 for (r = 0; r < p.col_numShapes; r++)
217 {
218 rtmp0 = 0;
219 rtmp1 = 0;
220 for (q = 0; q < p.numQuadTotal; q++)
221 {
222 rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
223 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
224 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
225 }
226 for (k = 0; k < p.numEqu; k++)
227 {
228 for (m = 0; m < p.numComp; m++)
229 {
230 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
231 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
232 rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
233 }
234 }
235 }
236 }
237 }
238 }
239 /**************************************************************/
240 /* process C: */
241 /**************************************************************/
242 C_p = getSampleDataRO(C, e);
243 if (NULL != C_p)
244 {
245 add_EM_S = TRUE;
246 if (extendedC)
247 {
248 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
249 for (s = 0; s < p.row_numShapes; s++)
250 {
251 for (r = 0; r < p.col_numShapes; r++)
252 {
253 for (k = 0; k < p.numEqu; k++)
254 {
255 for (m = 0; m < p.numComp; m++)
256 {
257 rtmp = 0;
258 for (q = 0; q < p.numQuadTotal; q++)
259 {
260 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
261 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
262 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
263 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
264 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
265 }
266 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
267 rtmp;
268 }
269 }
270 }
271 }
272 } else
273 {
274 for (s = 0; s < p.row_numShapes; s++)
275 {
276 for (r = 0; r < p.col_numShapes; r++)
277 {
278 rtmp0 = 0;
279 rtmp1 = 0;
280 for (q = 0; q < p.numQuadTotal; q++)
281 {
282 rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
283 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
284 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
285 }
286 for (k = 0; k < p.numEqu; k++)
287 {
288 for (m = 0; m < p.numComp; m++)
289 {
290 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
291 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
292 rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
293 }
294 }
295 }
296 }
297 }
298 }
299 /************************************************************* */
300 /* process D */
301 /**************************************************************/
302 D_p = getSampleDataRO(D, e);
303 if (NULL != D_p)
304 {
305 add_EM_S = TRUE;
306 if (extendedD)
307 {
308 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
309 for (s = 0; s < p.row_numShapes; s++)
310 {
311 for (r = 0; r < p.col_numShapes; r++)
312 {
313 for (k = 0; k < p.numEqu; k++)
314 {
315 for (m = 0; m < p.numComp; m++)
316 {
317 rtmp = 0;
318 for (q = 0; q < p.numQuadTotal; q++)
319 {
320 rtmp +=
321 vol * S[INDEX2(s, q, p.row_numShapes)] *
322 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
323 S[INDEX2(r, q, p.row_numShapes)];
324 }
325 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
326 rtmp;
327 }
328 }
329 }
330 }
331 } else
332 {
333 for (s = 0; s < p.row_numShapes; s++)
334 {
335 for (r = 0; r < p.col_numShapes; r++)
336 {
337 rtmp = 0;
338 for (q = 0; q < p.numQuadTotal; q++)
339 rtmp +=
340 vol * S[INDEX2(s, q, p.row_numShapes)] *
341 S[INDEX2(r, q, p.row_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.row_numShapesTotal)] +=
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.numQuadTotal)]);
364 for (s = 0; s < p.row_numShapes; s++)
365 {
366 for (k = 0; k < p.numEqu; k++)
367 {
368 rtmp = 0;
369 for (q = 0; q < p.numQuadTotal; q++)
370 {
371 rtmp +=
372 vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
373 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
375 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376 }
377 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378 }
379 }
380 } else
381 {
382 for (s = 0; s < p.row_numShapes; s++)
383 {
384 rtmp0 = 0;
385 rtmp1 = 0;
386 for (q = 0; q < p.numQuadTotal; q++)
387 {
388 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
389 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
390 }
391 for (k = 0; k < p.numEqu; k++)
392 EM_F[INDEX2(k, s, p.numEqu)] +=
393 rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
394 }
395 }
396 }
397 /**************************************************************/
398 /* process Y: */
399 /**************************************************************/
400 Y_p = getSampleDataRO(Y, e);
401 if (NULL != Y_p)
402 {
403 add_EM_F = TRUE;
404 if (extendedY)
405 {
406 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
407 for (s = 0; s < p.row_numShapes; s++)
408 {
409 for (k = 0; k < p.numEqu; k++)
410 {
411 rtmp = 0;
412 for (q = 0; q < p.numQuadTotal; q++)
413 rtmp +=
414 vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
415 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416 }
417 }
418 } else
419 {
420 for (s = 0; s < p.row_numShapes; s++)
421 {
422 rtmp = 0;
423 for (q = 0; q < p.numQuadTotal; q++)
424 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
425 for (k = 0; k < p.numEqu; k++)
426 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
427 }
428 }
429 }
430 /***********************************************************************************************/
431 /* add the element matrices onto the matrix and right hand side */
432 /***********************************************************************************************/
433 for (q = 0; q < p.row_numShapesTotal; q++)
434 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
435
436 if (add_EM_F)
437 Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
438 p.row_DOF_UpperBound);
439 if (add_EM_S)
440 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
441 p.col_numShapesTotal, row_index, p.numComp, EM_S);
442
443 } /* end color check */
444 } /* end element loop */
445 } /* end color loop */
446
447 THREAD_MEMFREE(EM_S);
448 THREAD_MEMFREE(EM_F);
449 THREAD_MEMFREE(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