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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3187 - (hide annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 14971 byte(s)
Enforcing indent -linux -i4 -bl -bli0 -l120

1 gross 798
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 gross 798 /**************************************************************/
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 gross 853 #ifdef _OPENMP
38     #include <omp.h>
39     #endif
40 gross 798
41     /**************************************************************/
42    
43 jfenwick 3187 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 gross 798
49 jfenwick 3187 #define DIM 2
50 gross 798 index_t color;
51 jfenwick 3136 dim_t e;
52 gross 2748 __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 jfenwick 3187 double *EM_S, *EM_F, *DSDX;
54 gross 853 index_t *row_index;
55 jfenwick 3187 register dim_t q, s, r, k, m;
56 gross 853 register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
57     bool_t add_EM_F, add_EM_S;
58    
59 jfenwick 3187 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 gross 798
70 jfenwick 3187 #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 gross 798 {
72 jfenwick 2271
73 jfenwick 3187 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 gross 798
77 jfenwick 3187 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
78     {
79 gross 2748
80 jfenwick 3187 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 gross 2748
89 jfenwick 3187 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 gross 2748
104 jfenwick 3187 /**************************************************************/
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 gross 798
436 jfenwick 3187 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 gross 798 }
454 jfenwick 3187
455 gross 798 /*
456     * $Log$
457     */

  ViewVC Help
Powered by ViewVC 1.1.26