/[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 3223 - (hide annotations)
Wed Sep 29 05:02:52 2010 UTC (8 years, 9 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_Single2_3D.c
File MIME type: text/plain
File size: 14048 byte(s)
Remove some - but not all - nasty hacks

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

  ViewVC Help
Powered by ViewVC 1.1.26