/[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 3221 - (hide annotations)
Wed Sep 29 01:00:21 2010 UTC (8 years, 6 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_Single2_1D.c
File MIME type: text/plain
File size: 10907 byte(s)
Comment stripping

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 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 jfenwick 3187 void Dudley_Assemble_PDE_Single2_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 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     } else
128     {
129     for (s = 0; s < p.row_numShapes; s++)
130     {
131 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
132 jfenwick 3187 {
133     rtmp = 0;
134     for (q = 0; q < p.numQuadTotal; q++)
135     rtmp +=
136     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
137     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
138     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
139     rtmp * A_p[INDEX2(0, 0, DIM)];
140     }
141     }
142     }
143     }
144     /**************************************************************/
145     /* process B: */
146     /**************************************************************/
147     if (NULL != B_p)
148     {
149     add_EM_S = TRUE;
150     if (extendedB)
151     {
152     B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
153     for (s = 0; s < p.row_numShapes; s++)
154     {
155 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
156 jfenwick 3187 {
157     rtmp = 0;
158     for (q = 0; q < p.numQuadTotal; q++)
159     {
160     rtmp +=
161     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
162     B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
163     }
164     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
165     }
166     }
167     } else
168     {
169     for (s = 0; s < p.row_numShapes; s++)
170     {
171 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
172 jfenwick 3187 {
173     rtmp = 0;
174     for (q = 0; q < p.numQuadTotal; q++)
175     rtmp +=
176     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
177     S[INDEX2(r, q, p.row_numShapes)];
178     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
179     rtmp * B_p[0];
180     }
181     }
182     }
183     }
184     /**************************************************************/
185     /* process C: */
186     /**************************************************************/
187     if (NULL != C_p)
188     {
189     add_EM_S = TRUE;
190     if (extendedC)
191     {
192     C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
193     for (s = 0; s < p.row_numShapes; s++)
194     {
195 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
196 jfenwick 3187 {
197     rtmp = 0;
198     for (q = 0; q < p.numQuadTotal; q++)
199     {
200     rtmp +=
201     vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
202     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
203     }
204     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
205     }
206     }
207     } else
208     {
209     for (s = 0; s < p.row_numShapes; s++)
210     {
211 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
212 jfenwick 3187 {
213     rtmp = 0;
214     for (q = 0; q < p.numQuadTotal; q++)
215     rtmp +=
216     vol * S[INDEX2(s, q, p.row_numShapes)] *
217     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
218     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
219     rtmp * C_p[0];
220     }
221     }
222     }
223     }
224     /************************************************************* */
225     /* process D */
226     /**************************************************************/
227     if (NULL != D_p)
228     {
229     add_EM_S = TRUE;
230     if (extendedD)
231     {
232     D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
233     for (s = 0; s < p.row_numShapes; s++)
234     {
235 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
236 jfenwick 3187 {
237     rtmp = 0;
238     for (q = 0; q < p.numQuadTotal; q++)
239     {
240     rtmp +=
241     vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
242     S[INDEX2(r, q, p.row_numShapes)];
243     }
244     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
245     }
246     }
247     } else
248     {
249     for (s = 0; s < p.row_numShapes; s++)
250     {
251 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
252 jfenwick 3187 {
253     rtmp = 0;
254     for (q = 0; q < p.numQuadTotal; q++)
255     rtmp +=
256     vol * S[INDEX2(s, q, p.row_numShapes)] *
257     S[INDEX2(r, q, p.row_numShapes)];
258     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
259     rtmp * D_p[0];
260     }
261     }
262     }
263     }
264     /**************************************************************/
265     /* process X: */
266     /**************************************************************/
267     if (NULL != X_p)
268     {
269     add_EM_F = TRUE;
270     if (extendedX)
271     {
272     X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
273     for (s = 0; s < p.row_numShapes; s++)
274     {
275     rtmp = 0;
276     for (q = 0; q < p.numQuadTotal; q++)
277     rtmp +=
278     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
279     X_q[INDEX2(0, q, DIM)];
280     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
281     }
282     } else
283     {
284     for (s = 0; s < p.row_numShapes; s++)
285     {
286     rtmp = 0;
287     for (q = 0; q < p.numQuadTotal; q++)
288     rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
289     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
290     }
291     }
292     }
293     /**************************************************************/
294     /* process Y: */
295     /**************************************************************/
296     if (NULL != Y_p)
297     {
298     add_EM_F = TRUE;
299     if (extendedY)
300     {
301     Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
302     for (s = 0; s < p.row_numShapes; s++)
303     {
304     rtmp = 0;
305     for (q = 0; q < p.numQuadTotal; q++)
306     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
307     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
308     }
309     } else
310     {
311     for (s = 0; s < p.row_numShapes; s++)
312     {
313     rtmp = 0;
314     for (q = 0; q < p.numQuadTotal; q++)
315     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
316     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
317     }
318     }
319     }
320     /***********************************************************************************************/
321     /* add the element matrices onto the matrix and right hand side */
322     /***********************************************************************************************/
323     for (q = 0; q < p.row_numShapesTotal; q++)
324     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
325 gross 798
326 jfenwick 3187 if (add_EM_F)
327     Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
328     p.row_DOF_UpperBound);
329     if (add_EM_S)
330     Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
331 jfenwick 3203 p.row_numShapesTotal, row_index, p.numComp, EM_S);
332 jfenwick 3187
333     } /* end color check */
334     } /* end element loop */
335     } /* end color loop */
336    
337     THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
338     THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
339     THREAD_MEMFREE(row_index);
340    
341     } /* end of pointer check */
342     } /* end parallel region */
343 gross 798 }
344 jfenwick 3187
345 gross 798 /*
346     * $Log$
347     */

  ViewVC Help
Powered by ViewVC 1.1.26