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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26