/[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 3981 - (hide annotations)
Fri Sep 21 02:47:54 2012 UTC (6 years, 6 months ago) by jfenwick
Original Path: trunk/dudley/src/Assemble_PDE_Single2_1D.c
File MIME type: text/plain
File size: 11449 byte(s)
First pass of updating copyright notices
1 gross 798
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 3911 * Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
7     * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10     *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 ksteube 1312
16 jfenwick 3981 /************************************************************************************/
17 gross 798
18     /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
19     /* the shape functions for test and solution must be identical */
20    
21     /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
22    
23     /* in a 1D domain. The shape functions for test and solution must be identical */
24     /* and row_NS == row_NN */
25    
26     /* Shape of the coefficients: */
27    
28     /* A = 1 x 1 */
29     /* B = 1 */
30     /* C = 1 */
31     /* D = scalar */
32     /* X = 1 */
33     /* Y = scalar */
34    
35 jfenwick 3981 /************************************************************************************/
36 gross 798
37     #include "Assemble.h"
38     #include "Util.h"
39 gross 853 #ifdef _OPENMP
40     #include <omp.h>
41     #endif
42 gross 798
43 jfenwick 3981 /************************************************************************************/
44 gross 798
45 caltinay 3247 void Dudley_Assemble_PDE_Single2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46 jfenwick 3187 Paso_SystemMatrix * Mat, escriptDataC * F,
47     escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48     escriptDataC * X, escriptDataC * Y)
49     {
50 gross 798
51 jfenwick 3187 #define DIM 1
52 gross 798 index_t color;
53 jfenwick 3136 dim_t e;
54 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;
55 jfenwick 3184 double *EM_S, *EM_F, *DSDX;
56 gross 853 index_t *row_index;
57 jfenwick 3187 register dim_t q, s, r;
58 gross 853 register double rtmp;
59     bool_t add_EM_F, add_EM_S;
60    
61 jfenwick 3187 bool_t extendedA = isExpanded(A);
62     bool_t extendedB = isExpanded(B);
63     bool_t extendedC = isExpanded(C);
64     bool_t extendedD = isExpanded(D);
65     bool_t extendedX = isExpanded(X);
66     bool_t extendedY = isExpanded(Y);
67     double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68     double *S = p.row_jac->BasisFunctions->S;
69 jfenwick 3203 dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal;
70 jfenwick 3187 dim_t len_EM_F = p.row_numShapesTotal;
71 gross 798
72 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)
73 gross 798 {
74 jfenwick 3187 EM_S = THREAD_MEMALLOC(len_EM_S, double);
75     EM_F = THREAD_MEMALLOC(len_EM_F, double);
76     row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
77 gross 853
78 jfenwick 3187 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79     {
80 jfenwick 2271
81 jfenwick 3187 for (color = elements->minColor; color <= elements->maxColor; color++)
82     {
83     /* open loop over all elements: */
84     #pragma omp for private(e) schedule(static)
85     for (e = 0; e < elements->numElements; e++)
86     {
87     if (elements->Color[e] == color)
88     {
89 gross 853
90 jfenwick 3187 A_p = getSampleDataRO(A, e);
91     C_p = getSampleDataRO(C, e);
92     B_p = getSampleDataRO(B, e);
93     D_p = getSampleDataRO(D, e);
94     X_p = getSampleDataRO(X, e);
95     Y_p = getSampleDataRO(Y, e);
96 gross 2748
97 jfenwick 3187 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
98     DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
99     for (q = 0; q < len_EM_S; ++q)
100     EM_S[q] = 0;
101     for (q = 0; q < len_EM_F; ++q)
102     EM_F[q] = 0;
103     add_EM_F = FALSE;
104     add_EM_S = FALSE;
105 jfenwick 3981 /************************************************************************************/
106 jfenwick 3187 /* process A: */
107 jfenwick 3981 /************************************************************************************/
108 jfenwick 3187 if (NULL != A_p)
109     {
110     add_EM_S = TRUE;
111     if (extendedA)
112     {
113     A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, 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     rtmp = 0;
119     for (q = 0; q < p.numQuadTotal; q++)
120     {
121     rtmp +=
122     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
123     A_q[INDEX3(0, 0, q, DIM, DIM)] *
124     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
125     }
126     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
127     }
128     }
129 jfenwick 3224 }
130     else
131 jfenwick 3187 {
132     for (s = 0; s < p.row_numShapes; s++)
133     {
134 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
135 jfenwick 3187 {
136     rtmp = 0;
137     for (q = 0; q < p.numQuadTotal; q++)
138     rtmp +=
139     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
140     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
141     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
142     rtmp * A_p[INDEX2(0, 0, DIM)];
143     }
144     }
145     }
146     }
147 jfenwick 3981 /************************************************************************************/
148 jfenwick 3187 /* process B: */
149 jfenwick 3981 /************************************************************************************/
150 jfenwick 3187 if (NULL != B_p)
151     {
152     add_EM_S = TRUE;
153     if (extendedB)
154     {
155     B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
156     for (s = 0; s < p.row_numShapes; s++)
157     {
158 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
159 jfenwick 3187 {
160     rtmp = 0;
161     for (q = 0; q < p.numQuadTotal; q++)
162     {
163     rtmp +=
164     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
165     B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
166     }
167     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
168     }
169     }
170 jfenwick 3224 }
171     else
172 jfenwick 3187 {
173     for (s = 0; s < p.row_numShapes; s++)
174     {
175 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
176 jfenwick 3187 {
177     rtmp = 0;
178     for (q = 0; q < p.numQuadTotal; q++)
179     rtmp +=
180     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
181     S[INDEX2(r, q, p.row_numShapes)];
182     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
183     rtmp * B_p[0];
184     }
185     }
186     }
187     }
188 jfenwick 3981 /************************************************************************************/
189 jfenwick 3187 /* process C: */
190 jfenwick 3981 /************************************************************************************/
191 jfenwick 3187 if (NULL != C_p)
192     {
193     add_EM_S = TRUE;
194     if (extendedC)
195     {
196     C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
197     for (s = 0; s < p.row_numShapes; s++)
198     {
199 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
200 jfenwick 3187 {
201     rtmp = 0;
202     for (q = 0; q < p.numQuadTotal; q++)
203     {
204     rtmp +=
205     vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
206     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
207     }
208     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
209     }
210     }
211 jfenwick 3224 }
212     else
213 jfenwick 3187 {
214     for (s = 0; s < p.row_numShapes; s++)
215     {
216 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
217 jfenwick 3187 {
218     rtmp = 0;
219     for (q = 0; q < p.numQuadTotal; q++)
220     rtmp +=
221     vol * S[INDEX2(s, q, p.row_numShapes)] *
222     DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
223     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
224     rtmp * C_p[0];
225     }
226     }
227     }
228     }
229 jfenwick 3981 /*********************************************************************************** */
230 jfenwick 3187 /* process D */
231 jfenwick 3981 /************************************************************************************/
232 jfenwick 3187 if (NULL != D_p)
233     {
234     add_EM_S = TRUE;
235     if (extendedD)
236     {
237     D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
238     for (s = 0; s < p.row_numShapes; s++)
239     {
240 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
241 jfenwick 3187 {
242     rtmp = 0;
243     for (q = 0; q < p.numQuadTotal; q++)
244     {
245     rtmp +=
246     vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
247     S[INDEX2(r, q, p.row_numShapes)];
248     }
249     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
250     }
251     }
252 jfenwick 3224 }
253     else
254 jfenwick 3187 {
255     for (s = 0; s < p.row_numShapes; s++)
256     {
257 jfenwick 3203 for (r = 0; r < p.row_numShapes; r++)
258 jfenwick 3187 {
259     rtmp = 0;
260     for (q = 0; q < p.numQuadTotal; q++)
261     rtmp +=
262     vol * S[INDEX2(s, q, p.row_numShapes)] *
263     S[INDEX2(r, q, p.row_numShapes)];
264     EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
265     rtmp * D_p[0];
266     }
267     }
268     }
269     }
270 jfenwick 3981 /************************************************************************************/
271 jfenwick 3187 /* process X: */
272 jfenwick 3981 /************************************************************************************/
273 jfenwick 3187 if (NULL != X_p)
274     {
275     add_EM_F = TRUE;
276     if (extendedX)
277     {
278     X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
279     for (s = 0; s < p.row_numShapes; s++)
280     {
281     rtmp = 0;
282     for (q = 0; q < p.numQuadTotal; q++)
283     rtmp +=
284     vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
285     X_q[INDEX2(0, q, DIM)];
286     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
287     }
288 jfenwick 3224 }
289     else
290 jfenwick 3187 {
291     for (s = 0; s < p.row_numShapes; s++)
292     {
293     rtmp = 0;
294     for (q = 0; q < p.numQuadTotal; q++)
295     rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
296     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
297     }
298     }
299     }
300 jfenwick 3981 /************************************************************************************/
301 jfenwick 3187 /* process Y: */
302 jfenwick 3981 /************************************************************************************/
303 jfenwick 3187 if (NULL != Y_p)
304     {
305     add_EM_F = TRUE;
306     if (extendedY)
307     {
308     Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
309     for (s = 0; s < p.row_numShapes; s++)
310     {
311     rtmp = 0;
312     for (q = 0; q < p.numQuadTotal; q++)
313     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
314     EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
315     }
316 jfenwick 3224 }
317     else
318 jfenwick 3187 {
319     for (s = 0; s < p.row_numShapes; s++)
320     {
321     rtmp = 0;
322     for (q = 0; q < p.numQuadTotal; q++)
323     rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
324     EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
325     }
326     }
327     }
328 jfenwick 3981 /*********************************************************************************************************************/
329 jfenwick 3187 /* add the element matrices onto the matrix and right hand side */
330 jfenwick 3981 /*********************************************************************************************************************/
331 jfenwick 3187 for (q = 0; q < p.row_numShapesTotal; q++)
332     row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
333 gross 798
334 jfenwick 3187 if (add_EM_F)
335     Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
336     p.row_DOF_UpperBound);
337     if (add_EM_S)
338     Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
339 jfenwick 3203 p.row_numShapesTotal, row_index, p.numComp, EM_S);
340 jfenwick 3187
341     } /* end color check */
342     } /* end element loop */
343     } /* end color loop */
344    
345     THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
346     THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
347     THREAD_MEMFREE(row_index);
348    
349     } /* end of pointer check */
350     } /* end parallel region */
351 gross 798 }
352 jfenwick 3187
353 gross 798 /*
354     * $Log$
355     */

  ViewVC Help
Powered by ViewVC 1.1.26