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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3203 - (hide annotations)
Thu Sep 23 23:51:26 2010 UTC (8 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 14664 byte(s)
Collapsing some col_params

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

  ViewVC Help
Powered by ViewVC 1.1.26