/[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 3203 - (hide annotations)
Thu Sep 23 23:51:26 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: 10990 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 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 jfenwick 3184 // Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
97 jfenwick 3187 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98     for (q = 0; q < len_EM_S; ++q)
99     EM_S[q] = 0;
100     for (q = 0; q < len_EM_F; ++q)
101     EM_F[q] = 0;
102     add_EM_F = FALSE;
103     add_EM_S = FALSE;
104     /**************************************************************/
105     /* process A: */
106     /**************************************************************/
107     if (NULL != A_p)
108     {
109     add_EM_S = TRUE;
110     if (extendedA)
111     {
112     A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
113     for (s = 0; s < p.row_numShapes; s++)
114     {
115 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
116 jfenwick 3187 {
117     rtmp = 0;
118     for (q = 0; q < p.numQuadTotal; q++)
119     {
120     rtmp +=
121     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
122     A_q[INDEX3(0, 0, q, DIM, DIM)] *
123     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
124     }
125     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
126     }
127     }
128     } else
129     {
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     } else
169     {
170     for (s = 0; s < p.row_numShapes; s++)
171     {
172 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
173 jfenwick 3187 {
174     rtmp = 0;
175     for (q = 0; q < p.numQuadTotal; q++)
176     rtmp +=
177     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
178     S[INDEX2(r, q, p.row_numShapes)];
179     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
180     rtmp * B_p[0];
181     }
182     }
183     }
184     }
185     /**************************************************************/
186     /* process C: */
187     /**************************************************************/
188     if (NULL != C_p)
189     {
190     add_EM_S = TRUE;
191     if (extendedC)
192     {
193     C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
194     for (s = 0; s < p.row_numShapes; s++)
195     {
196 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
197 jfenwick 3187 {
198     rtmp = 0;
199     for (q = 0; q < p.numQuadTotal; q++)
200     {
201     rtmp +=
202     vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
203     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
204     }
205     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
206     }
207     }
208     } else
209     {
210     for (s = 0; s < p.row_numShapes; s++)
211     {
212 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
213 jfenwick 3187 {
214     rtmp = 0;
215     for (q = 0; q < p.numQuadTotal; q++)
216     rtmp +=
217     vol * S[INDEX2(s, q, p.row_numShapes)] *
218     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
219     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
220     rtmp * C_p[0];
221     }
222     }
223     }
224     }
225     /************************************************************* */
226     /* process D */
227     /**************************************************************/
228     if (NULL != D_p)
229     {
230     add_EM_S = TRUE;
231     if (extendedD)
232     {
233     D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
234     for (s = 0; s < p.row_numShapes; s++)
235     {
236 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
237 jfenwick 3187 {
238     rtmp = 0;
239     for (q = 0; q < p.numQuadTotal; q++)
240     {
241     rtmp +=
242     vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
243     S[INDEX2(r, q, p.row_numShapes)];
244     }
245     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
246     }
247     }
248     } else
249     {
250     for (s = 0; s < p.row_numShapes; s++)
251     {
252 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
253 jfenwick 3187 {
254     rtmp = 0;
255     for (q = 0; q < p.numQuadTotal; q++)
256     rtmp +=
257     vol * S[INDEX2(s, q, p.row_numShapes)] *
258     S[INDEX2(r, q, p.row_numShapes)];
259     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
260     rtmp * D_p[0];
261     }
262     }
263     }
264     }
265     /**************************************************************/
266     /* process X: */
267     /**************************************************************/
268     if (NULL != X_p)
269     {
270     add_EM_F = TRUE;
271     if (extendedX)
272     {
273     X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
274     for (s = 0; s < p.row_numShapes; s++)
275     {
276     rtmp = 0;
277     for (q = 0; q < p.numQuadTotal; q++)
278     rtmp +=
279     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
280     X_q[INDEX2(0, q, DIM)];
281     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
282     }
283     } else
284     {
285     for (s = 0; s < p.row_numShapes; s++)
286     {
287     rtmp = 0;
288     for (q = 0; q < p.numQuadTotal; q++)
289     rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
290     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
291     }
292     }
293     }
294     /**************************************************************/
295     /* process Y: */
296     /**************************************************************/
297     if (NULL != Y_p)
298     {
299     add_EM_F = TRUE;
300     if (extendedY)
301     {
302     Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
303     for (s = 0; s < p.row_numShapes; s++)
304     {
305     rtmp = 0;
306     for (q = 0; q < p.numQuadTotal; q++)
307     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
308     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
309     }
310     } else
311     {
312     for (s = 0; s < p.row_numShapes; s++)
313     {
314     rtmp = 0;
315     for (q = 0; q < p.numQuadTotal; q++)
316     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
317     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
318     }
319     }
320     }
321     /***********************************************************************************************/
322     /* add the element matrices onto the matrix and right hand side */
323     /***********************************************************************************************/
324     for (q = 0; q < p.row_numShapesTotal; q++)
325     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
326 gross 798
327 jfenwick 3187 if (add_EM_F)
328     Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
329     p.row_DOF_UpperBound);
330     if (add_EM_S)
331     Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
332 jfenwick 3203 p.row_numShapesTotal, row_index, p.numComp, EM_S);
333 jfenwick 3187
334     } /* end color check */
335     } /* end element loop */
336     } /* end color loop */
337    
338     THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
339     THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
340     THREAD_MEMFREE(row_index);
341    
342     } /* end of pointer check */
343     } /* end parallel region */
344 gross 798 }
345 jfenwick 3187
346 gross 798 /*
347     * $Log$
348     */

  ViewVC Help
Powered by ViewVC 1.1.26