/[escript]/trunk/dudley/src/Assemble_PDE_Single2_3D.c
ViewVC logotype

Annotation of /trunk/dudley/src/Assemble_PDE_Single2_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3187 - (hide annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 10 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_Single2_3D.c
File MIME type: text/plain
File size: 14943 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_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
20    
21     /* in a 3D 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 = 3 x 3 */
27     /* B = 3 */
28     /* C = 3 */
29     /* D = scalar */
30     /* X = 3 */
31     /* Y = scalar */
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_Single2_3D(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 3
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 3184 double *EM_S, *EM_F, *DSDX;
54 gross 853 index_t *row_index;
55 jfenwick 3187 register dim_t q, s, r;
56 gross 853 register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
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;
68     dim_t len_EM_F = p.row_numShapesTotal;
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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
71 gross 798 {
72 jfenwick 3187 EM_S = THREAD_MEMALLOC(len_EM_S, double);
73     EM_F = THREAD_MEMALLOC(len_EM_F, double);
74     row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
75 gross 798
76 jfenwick 3187 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
77     {
78 gross 2748
79 jfenwick 3187 for (color = elements->minColor; color <= elements->maxColor; color++)
80     {
81     /* open loop over all elements: */
82     #pragma omp for private(e) schedule(static)
83     for (e = 0; e < elements->numElements; e++)
84     {
85     if (elements->Color[e] == color)
86     {
87 gross 2748
88 jfenwick 3187 A_p = getSampleDataRO(A, e);
89     B_p = getSampleDataRO(B, e);
90     C_p = getSampleDataRO(C, e);
91     D_p = getSampleDataRO(D, e);
92     X_p = getSampleDataRO(X, e);
93     Y_p = getSampleDataRO(Y, e);
94    
95 jfenwick 3184 // Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
96 jfenwick 3187 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 gross 2748
105 jfenwick 3187 /**************************************************************/
106     /* process A: */
107     /**************************************************************/
108     if (NULL != A_p)
109     {
110     add_EM_S = TRUE;
111     if (extendedA)
112     {
113     A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
114     for (s = 0; s < p.row_numShapes; s++)
115     {
116     for (r = 0; r < p.col_numShapes; r++)
117     {
118     rtmp = 0;
119     for (q = 0; q < p.numQuadTotal; q++)
120     {
121     rtmp +=
122     vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
123     A_q[INDEX3(0, 0, q, DIM, DIM)] *
124     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
125     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
126     A_q[INDEX3(0, 1, q, DIM, DIM)] *
127     DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
128     DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
129     A_q[INDEX3(0, 2, q, DIM, DIM)] *
130     DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
131     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
132     A_q[INDEX3(1, 0, q, DIM, DIM)] *
133     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
134     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
135     A_q[INDEX3(1, 1, q, DIM, DIM)] *
136     DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
137     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
138     A_q[INDEX3(1, 2, q, DIM, DIM)] *
139     DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
140     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
141     A_q[INDEX3(2, 0, q, DIM, DIM)] *
142     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
143     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
144     A_q[INDEX3(2, 1, q, DIM, DIM)] *
145     DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
146     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
147     A_q[INDEX3(2, 2, q, DIM, DIM)] *
148     DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
149     }
150     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
151     }
152     }
153     } else
154     {
155     for (s = 0; s < p.row_numShapes; s++)
156     {
157     for (r = 0; r < p.col_numShapes; r++)
158     {
159     rtmp00 = 0;
160     rtmp01 = 0;
161     rtmp02 = 0;
162     rtmp10 = 0;
163     rtmp11 = 0;
164     rtmp12 = 0;
165     rtmp20 = 0;
166     rtmp21 = 0;
167     rtmp22 = 0;
168     for (q = 0; q < p.numQuadTotal; q++)
169     {
170 gross 2748
171 jfenwick 3187 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
172     rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
173     rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
174     rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
175 gross 2748
176 jfenwick 3187 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
177     rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
178     rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
179     rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
180 gross 2748
181 jfenwick 3187 rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
182     rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
183     rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
184     rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
185     }
186     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
187     rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
188     rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
189     rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
190     rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
191     rtmp22 * A_p[INDEX2(2, 2, DIM)];
192     }
193     }
194     }
195     }
196     /**************************************************************/
197     /* process B: */
198     /**************************************************************/
199     if (NULL != B_p)
200     {
201     add_EM_S = TRUE;
202     if (extendedB)
203     {
204     B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
205     for (s = 0; s < p.row_numShapes; s++)
206     {
207     for (r = 0; r < p.col_numShapes; r++)
208     {
209     rtmp = 0;
210     for (q = 0; q < p.numQuadTotal; q++)
211     {
212     rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
213     (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
214     B_q[INDEX2(0, q, DIM)] +
215     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
216     B_q[INDEX2(1, q, DIM)] +
217     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
218     B_q[INDEX2(2, q, DIM)]);
219     }
220     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
221     }
222     }
223     } else
224     {
225     for (s = 0; s < p.row_numShapes; s++)
226     {
227     for (r = 0; r < p.col_numShapes; r++)
228     {
229     rtmp0 = 0;
230     rtmp1 = 0;
231     rtmp2 = 0;
232     for (q = 0; q < p.numQuadTotal; q++)
233     {
234     rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
235     rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
236     rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
237     rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
238     }
239     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
240     rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
241     }
242     }
243     }
244     }
245     /**************************************************************/
246     /* process C: */
247     /**************************************************************/
248     if (NULL != C_p)
249     {
250     add_EM_S = TRUE;
251     if (extendedC)
252     {
253     C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
254     for (s = 0; s < p.row_numShapes; s++)
255     {
256     for (r = 0; r < p.col_numShapes; r++)
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[INDEX2(0, q, DIM)] *
263     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264     C_q[INDEX2(1, q, DIM)] *
265     DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
266     C_q[INDEX2(2, q, DIM)] *
267     DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
268     }
269     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
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     rtmp2 = 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     rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
287     }
288     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
289     rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
290     }
291     }
292     }
293     }
294     /************************************************************* */
295     /* process D */
296     /**************************************************************/
297     if (NULL != D_p)
298     {
299     add_EM_S = TRUE;
300     if (extendedD)
301     {
302     D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
303     for (s = 0; s < p.row_numShapes; s++)
304     {
305     for (r = 0; r < p.col_numShapes; r++)
306     {
307     rtmp = 0;
308     for (q = 0; q < p.numQuadTotal; q++)
309     rtmp +=
310     vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
311     S[INDEX2(r, q, p.row_numShapes)];
312     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
313     }
314     }
315     } else
316     {
317     for (s = 0; s < p.row_numShapes; s++)
318     {
319     for (r = 0; r < p.col_numShapes; r++)
320     {
321     rtmp = 0;
322     for (q = 0; q < p.numQuadTotal; q++)
323     rtmp +=
324     vol * S[INDEX2(s, q, p.row_numShapes)] *
325     S[INDEX2(r, q, p.row_numShapes)];
326     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327     rtmp * D_p[0];
328     }
329     }
330     }
331     }
332     /**************************************************************/
333     /* process X: */
334     /**************************************************************/
335     if (NULL != X_p)
336     {
337     add_EM_F = TRUE;
338     if (extendedX)
339     {
340     X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
341     for (s = 0; s < p.row_numShapes; s++)
342     {
343     rtmp = 0;
344     for (q = 0; q < p.numQuadTotal; q++)
345     {
346     rtmp +=
347     vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
348     X_q[INDEX2(0, q, DIM)] +
349     DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
350     X_q[INDEX2(1, q, DIM)] +
351     DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
352     X_q[INDEX2(2, q, DIM)]);
353     }
354     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
355     }
356     } else
357     {
358     for (s = 0; s < p.row_numShapes; s++)
359     {
360     rtmp0 = 0;
361     rtmp1 = 0;
362     rtmp2 = 0;
363     for (q = 0; q < p.numQuadTotal; q++)
364     {
365     rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
366     rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
367     rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
368     }
369     EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
370     }
371     }
372     }
373     /**************************************************************/
374     /* process Y: */
375     /**************************************************************/
376     if (NULL != Y_p)
377     {
378     add_EM_F = TRUE;
379     if (extendedY)
380     {
381     Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
382     for (s = 0; s < p.row_numShapes; s++)
383     {
384     rtmp = 0;
385     for (q = 0; q < p.numQuadTotal; q++)
386     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
387     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
388     }
389     } else
390     {
391     for (s = 0; s < p.row_numShapes; s++)
392     {
393     rtmp = 0;
394     for (q = 0; q < p.numQuadTotal; q++)
395     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
396     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
397     }
398     }
399     }
400     /***********************************************************************************************/
401     /* add the element matrices onto the matrix and right hand side */
402     /***********************************************************************************************/
403 gross 798
404 jfenwick 3187 for (q = 0; q < p.row_numShapesTotal; q++)
405     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
406    
407     if (add_EM_F)
408     Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
409     p.row_DOF_UpperBound);
410     if (add_EM_S)
411     Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
412     p.col_numShapesTotal, row_index, p.numComp, EM_S);
413    
414     } /* end color check */
415     } /* end element loop */
416     } /* end color loop */
417    
418     THREAD_MEMFREE(EM_S);
419     THREAD_MEMFREE(EM_F);
420     THREAD_MEMFREE(row_index);
421    
422     } /* end of pointer check */
423     } /* end parallel region */
424 gross 798 }

  ViewVC Help
Powered by ViewVC 1.1.26