/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp
ViewVC logotype

Annotation of /branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (6 years, 9 months ago) by jfenwick
Original Path: trunk/dudley/src/Assemble_PDE_Single2_1D.c
File MIME type: text/plain
File size: 10956 byte(s)
Copyright changes
1 gross 798
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 3911 * Copyright (c) 2003-2012 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 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 = 1 x 1 */
27     /* B = 1 */
28     /* C = 1 */
29     /* D = scalar */
30     /* X = 1 */
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 caltinay 3247 void Dudley_Assemble_PDE_Single2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
44 jfenwick 3187 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 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;
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_numShapesTotal * p.row_numShapesTotal;
68 jfenwick 3187 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,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 853
76 jfenwick 3187 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
77     {
78 jfenwick 2271
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 853
88 jfenwick 3187 A_p = getSampleDataRO(A, e);
89     C_p = getSampleDataRO(C, e);
90     B_p = getSampleDataRO(B, e);
91     D_p = getSampleDataRO(D, e);
92     X_p = getSampleDataRO(X, e);
93     Y_p = getSampleDataRO(Y, e);
94 gross 2748
95 jfenwick 3187 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     /**************************************************************/
104     /* process A: */
105     /**************************************************************/
106     if (NULL != A_p)
107     {
108     add_EM_S = TRUE;
109     if (extendedA)
110     {
111     A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
112     for (s = 0; s < p.row_numShapes; s++)
113     {
114 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
115 jfenwick 3187 {
116     rtmp = 0;
117     for (q = 0; q < p.numQuadTotal; q++)
118     {
119     rtmp +=
120     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
121     A_q[INDEX3(0, 0, q, DIM, DIM)] *
122     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
123     }
124     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
125     }
126     }
127 jfenwick 3224 }
128     else
129 jfenwick 3187 {
130     for (s = 0; s < p.row_numShapes; s++)
131     {
132 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
133 jfenwick 3187 {
134     rtmp = 0;
135     for (q = 0; q < p.numQuadTotal; q++)
136     rtmp +=
137     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
138     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
139     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
140     rtmp * A_p[INDEX2(0, 0, DIM)];
141     }
142     }
143     }
144     }
145     /**************************************************************/
146     /* process B: */
147     /**************************************************************/
148     if (NULL != B_p)
149     {
150     add_EM_S = TRUE;
151     if (extendedB)
152     {
153     B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
154     for (s = 0; s < p.row_numShapes; s++)
155     {
156 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
157 jfenwick 3187 {
158     rtmp = 0;
159     for (q = 0; q < p.numQuadTotal; q++)
160     {
161     rtmp +=
162     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
163     B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
164     }
165     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
166     }
167     }
168 jfenwick 3224 }
169     else
170 jfenwick 3187 {
171     for (s = 0; s < p.row_numShapes; s++)
172     {
173 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
174 jfenwick 3187 {
175     rtmp = 0;
176     for (q = 0; q < p.numQuadTotal; q++)
177     rtmp +=
178     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
179     S[INDEX2(r, q, p.row_numShapes)];
180     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
181     rtmp * B_p[0];
182     }
183     }
184     }
185     }
186     /**************************************************************/
187     /* process C: */
188     /**************************************************************/
189     if (NULL != C_p)
190     {
191     add_EM_S = TRUE;
192     if (extendedC)
193     {
194     C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
195     for (s = 0; s < p.row_numShapes; s++)
196     {
197 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
198 jfenwick 3187 {
199     rtmp = 0;
200     for (q = 0; q < p.numQuadTotal; q++)
201     {
202     rtmp +=
203     vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
204     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
205     }
206     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
207     }
208     }
209 jfenwick 3224 }
210     else
211 jfenwick 3187 {
212     for (s = 0; s < p.row_numShapes; s++)
213     {
214 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
215 jfenwick 3187 {
216     rtmp = 0;
217     for (q = 0; q < p.numQuadTotal; q++)
218     rtmp +=
219     vol * S[INDEX2(s, q, p.row_numShapes)] *
220     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
221     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
222     rtmp * C_p[0];
223     }
224     }
225     }
226     }
227     /************************************************************* */
228     /* process D */
229     /**************************************************************/
230     if (NULL != D_p)
231     {
232     add_EM_S = TRUE;
233     if (extendedD)
234     {
235     D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
236     for (s = 0; s < p.row_numShapes; s++)
237     {
238 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
239 jfenwick 3187 {
240     rtmp = 0;
241     for (q = 0; q < p.numQuadTotal; q++)
242     {
243     rtmp +=
244     vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
245     S[INDEX2(r, q, p.row_numShapes)];
246     }
247     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
248     }
249     }
250 jfenwick 3224 }
251     else
252 jfenwick 3187 {
253     for (s = 0; s < p.row_numShapes; s++)
254     {
255 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
256 jfenwick 3187 {
257     rtmp = 0;
258     for (q = 0; q < p.numQuadTotal; q++)
259     rtmp +=
260     vol * S[INDEX2(s, q, p.row_numShapes)] *
261     S[INDEX2(r, q, p.row_numShapes)];
262     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
263     rtmp * D_p[0];
264     }
265     }
266     }
267     }
268     /**************************************************************/
269     /* process X: */
270     /**************************************************************/
271     if (NULL != X_p)
272     {
273     add_EM_F = TRUE;
274     if (extendedX)
275     {
276     X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
277     for (s = 0; s < p.row_numShapes; s++)
278     {
279     rtmp = 0;
280     for (q = 0; q < p.numQuadTotal; q++)
281     rtmp +=
282     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
283     X_q[INDEX2(0, q, DIM)];
284     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
285     }
286 jfenwick 3224 }
287     else
288 jfenwick 3187 {
289     for (s = 0; s < p.row_numShapes; s++)
290     {
291     rtmp = 0;
292     for (q = 0; q < p.numQuadTotal; q++)
293     rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
294     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
295     }
296     }
297     }
298     /**************************************************************/
299     /* process Y: */
300     /**************************************************************/
301     if (NULL != Y_p)
302     {
303     add_EM_F = TRUE;
304     if (extendedY)
305     {
306     Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
307     for (s = 0; s < p.row_numShapes; s++)
308     {
309     rtmp = 0;
310     for (q = 0; q < p.numQuadTotal; q++)
311     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
312     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
313     }
314 jfenwick 3224 }
315     else
316 jfenwick 3187 {
317     for (s = 0; s < p.row_numShapes; s++)
318     {
319     rtmp = 0;
320     for (q = 0; q < p.numQuadTotal; q++)
321     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
322     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
323     }
324     }
325     }
326     /***********************************************************************************************/
327     /* add the element matrices onto the matrix and right hand side */
328     /***********************************************************************************************/
329     for (q = 0; q < p.row_numShapesTotal; q++)
330     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
331 gross 798
332 jfenwick 3187 if (add_EM_F)
333     Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
334     p.row_DOF_UpperBound);
335     if (add_EM_S)
336     Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
337 jfenwick 3203 p.row_numShapesTotal, row_index, p.numComp, EM_S);
338 jfenwick 3187
339     } /* end color check */
340     } /* end element loop */
341     } /* end color loop */
342    
343     THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
344     THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
345     THREAD_MEMFREE(row_index);
346    
347     } /* end of pointer check */
348     } /* end parallel region */
349 gross 798 }
350 jfenwick 3187
351 gross 798 /*
352     * $Log$
353     */

  ViewVC Help
Powered by ViewVC 1.1.26