/[escript]/branches/domexper/dudley/src/Assemble_PDE_System2_2D.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_PDE_System2_2D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3202 - (show annotations)
Thu Sep 23 23:24:09 2010 UTC (8 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 15008 byte(s)


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 const double* S = p.shapeFns;
68 dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal * p.numEqu * p.numComp;
69 dim_t len_EM_F = p.row_numShapesTotal * p.numEqu;
70
71 #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)
72 {
73
74 EM_S = THREAD_MEMALLOC(len_EM_S, double);
75 EM_F = THREAD_MEMALLOC(len_EM_F, double);
76 row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
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 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98 for (q = 0; q < len_EM_S; ++q)
99 EM_S[q] = 0;
100 for (q = 0; q < len_EM_F; ++q)
101 EM_F[q] = 0;
102 add_EM_F = FALSE;
103 add_EM_S = FALSE;
104
105 /**************************************************************/
106 /* process A: */
107 /**************************************************************/
108 A_p = getSampleDataRO(A, e);
109 if (NULL != A_p)
110 {
111 add_EM_S = TRUE;
112 if (extendedA)
113 {
114 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
115 for (s = 0; s < p.row_numShapes; s++)
116 {
117 for (r = 0; r < p.col_numShapes; r++)
118 {
119 for (k = 0; k < p.numEqu; k++)
120 {
121 for (m = 0; m < p.numComp; m++)
122 {
123 rtmp = 0;
124 for (q = 0; q < p.numQuadTotal; q++)
125 {
126 rtmp +=
127 vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
128 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
129 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
130 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
131 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
132 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
133 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
134 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
135 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
136 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
137 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
138 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
139 }
140 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
141 rtmp;
142 }
143 }
144 }
145 }
146 } else
147 {
148 for (s = 0; s < p.row_numShapes; s++)
149 {
150 for (r = 0; r < p.col_numShapes; r++)
151 {
152 rtmp00 = 0;
153 rtmp01 = 0;
154 rtmp10 = 0;
155 rtmp11 = 0;
156 for (q = 0; q < p.numQuadTotal; q++)
157 {
158 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
159 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
160 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
161 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
162 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
163 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
164 }
165 for (k = 0; k < p.numEqu; k++)
166 {
167 for (m = 0; m < p.numComp; m++)
168 {
169 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
170 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
171 + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
172 + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
173 + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
174 }
175 }
176 }
177 }
178 }
179 }
180 /**************************************************************/
181 /* process B: */
182 /**************************************************************/
183 B_p = getSampleDataRO(B, e);
184 if (NULL != B_p)
185 {
186 add_EM_S = TRUE;
187 if (extendedB)
188 {
189 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
190 for (s = 0; s < p.row_numShapes; s++)
191 {
192 for (r = 0; r < p.col_numShapes; r++)
193 {
194 for (k = 0; k < p.numEqu; k++)
195 {
196 for (m = 0; m < p.numComp; m++)
197 {
198 rtmp = 0;
199 for (q = 0; q < p.numQuadTotal; q++)
200 {
201 rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
202 (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
203 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
204 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
205 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
206 }
207 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
208 rtmp;
209 }
210 }
211 }
212 }
213 } else
214 {
215 for (s = 0; s < p.row_numShapes; s++)
216 {
217 for (r = 0; r < p.col_numShapes; r++)
218 {
219 rtmp0 = 0;
220 rtmp1 = 0;
221 for (q = 0; q < p.numQuadTotal; q++)
222 {
223 rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
224 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
225 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
226 }
227 for (k = 0; k < p.numEqu; k++)
228 {
229 for (m = 0; m < p.numComp; m++)
230 {
231 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
232 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
233 rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
234 }
235 }
236 }
237 }
238 }
239 }
240 /**************************************************************/
241 /* process C: */
242 /**************************************************************/
243 C_p = getSampleDataRO(C, e);
244 if (NULL != C_p)
245 {
246 add_EM_S = TRUE;
247 if (extendedC)
248 {
249 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
250 for (s = 0; s < p.row_numShapes; s++)
251 {
252 for (r = 0; r < p.col_numShapes; r++)
253 {
254 for (k = 0; k < p.numEqu; k++)
255 {
256 for (m = 0; m < p.numComp; m++)
257 {
258 rtmp = 0;
259 for (q = 0; q < p.numQuadTotal; q++)
260 {
261 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
262 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
263 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
265 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)]);
266 }
267 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
268 rtmp;
269 }
270 }
271 }
272 }
273 } else
274 {
275 for (s = 0; s < p.row_numShapes; s++)
276 {
277 for (r = 0; r < p.col_numShapes; r++)
278 {
279 rtmp0 = 0;
280 rtmp1 = 0;
281 for (q = 0; q < p.numQuadTotal; q++)
282 {
283 rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
284 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
285 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
286 }
287 for (k = 0; k < p.numEqu; k++)
288 {
289 for (m = 0; m < p.numComp; m++)
290 {
291 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
292 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
293 rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
294 }
295 }
296 }
297 }
298 }
299 }
300 /************************************************************* */
301 /* process D */
302 /**************************************************************/
303 D_p = getSampleDataRO(D, e);
304 if (NULL != D_p)
305 {
306 add_EM_S = TRUE;
307 if (extendedD)
308 {
309 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
310 for (s = 0; s < p.row_numShapes; s++)
311 {
312 for (r = 0; r < p.col_numShapes; r++)
313 {
314 for (k = 0; k < p.numEqu; k++)
315 {
316 for (m = 0; m < p.numComp; m++)
317 {
318 rtmp = 0;
319 for (q = 0; q < p.numQuadTotal; q++)
320 {
321 rtmp +=
322 vol * S[INDEX2(s, q, p.row_numShapes)] *
323 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
324 S[INDEX2(r, q, p.row_numShapes)];
325 }
326 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327 rtmp;
328 }
329 }
330 }
331 }
332 } else
333 {
334 for (s = 0; s < p.row_numShapes; s++)
335 {
336 for (r = 0; r < p.col_numShapes; r++)
337 {
338 rtmp = 0;
339 for (q = 0; q < p.numQuadTotal; q++)
340 rtmp +=
341 vol * S[INDEX2(s, q, p.row_numShapes)] *
342 S[INDEX2(r, q, p.row_numShapes)];
343 for (k = 0; k < p.numEqu; k++)
344 {
345 for (m = 0; m < p.numComp; m++)
346 {
347 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
348 rtmp * D_p[INDEX2(k, m, p.numEqu)];
349 }
350 }
351 }
352 }
353 }
354 }
355 /**************************************************************/
356 /* process X: */
357 /**************************************************************/
358 X_p = getSampleDataRO(X, e);
359 if (NULL != X_p)
360 {
361 add_EM_F = TRUE;
362 if (extendedX)
363 {
364 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
365 for (s = 0; s < p.row_numShapes; s++)
366 {
367 for (k = 0; k < p.numEqu; k++)
368 {
369 rtmp = 0;
370 for (q = 0; q < p.numQuadTotal; q++)
371 {
372 rtmp +=
373 vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
374 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
375 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
376 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
377 }
378 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
379 }
380 }
381 } else
382 {
383 for (s = 0; s < p.row_numShapes; s++)
384 {
385 rtmp0 = 0;
386 rtmp1 = 0;
387 for (q = 0; q < p.numQuadTotal; q++)
388 {
389 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
390 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, 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.numQuadTotal)]);
408 for (s = 0; s < p.row_numShapes; s++)
409 {
410 for (k = 0; k < p.numEqu; k++)
411 {
412 rtmp = 0;
413 for (q = 0; q < p.numQuadTotal; q++)
414 rtmp +=
415 vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
416 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
417 }
418 }
419 } else
420 {
421 for (s = 0; s < p.row_numShapes; s++)
422 {
423 rtmp = 0;
424 for (q = 0; q < p.numQuadTotal; q++)
425 rtmp += vol * S[INDEX2(s, q, p.row_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.row_numShapesTotal; 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.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
439 p.row_DOF_UpperBound);
440 if (add_EM_S)
441 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
442 p.col_numShapesTotal, row_index, p.numComp, EM_S);
443
444 } /* end color check */
445 } /* end element loop */
446 } /* end color loop */
447
448 THREAD_MEMFREE(EM_S);
449 THREAD_MEMFREE(EM_F);
450 THREAD_MEMFREE(row_index);
451
452 } /* end of pointer check */
453 } /* end parallel region */
454 }
455
456 /*
457 * $Log$
458 */

  ViewVC Help
Powered by ViewVC 1.1.26