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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4257 - (hide annotations)
Wed Feb 27 03:42:40 2013 UTC (6 years, 2 months ago) by jfenwick
Original Path: branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.c
File MIME type: text/plain
File size: 14827 byte(s)
Some simple experiments for c++ Finley

1 gross 798
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4154 * Copyright (c) 2003-2013 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_{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 */
22    
23     /* u has p.numComp components in a 2D 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 = p.numEqu x 2 x p.numComp x 2 */
29     /* B = 2 x p.numEqu x p.numComp */
30     /* C = p.numEqu x 2 x p.numComp */
31     /* D = p.numEqu x p.numComp */
32     /* X = p.numEqu x 2 */
33     /* Y = p.numEqu */
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_System2_2D(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 2
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 3187 double *EM_S, *EM_F, *DSDX;
56 gross 853 index_t *row_index;
57 jfenwick 3187 register dim_t q, s, r, k, m;
58 gross 853 register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
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 jfenwick 3224 const double *S = p.shapeFns;
69 jfenwick 3204 dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
70     dim_t len_EM_F = p.numShapes * p.numEqu;
71 gross 798
72 jfenwick 3231 #pragma omp parallel private(color,EM_S, EM_F, 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, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)
73 gross 798 {
74 jfenwick 2271
75 jfenwick 3187 EM_S = THREAD_MEMALLOC(len_EM_S, double);
76     EM_F = THREAD_MEMALLOC(len_EM_F, double);
77 jfenwick 3204 row_index = THREAD_MEMALLOC(p.numShapes, index_t);
78 gross 798
79 jfenwick 3187 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
80     {
81 gross 2748
82 jfenwick 3187 for (color = elements->minColor; color <= elements->maxColor; color++)
83     {
84     /* open loop over all elements: */
85     #pragma omp for private(e) schedule(static)
86     for (e = 0; e < elements->numElements; e++)
87     {
88     if (elements->Color[e] == color)
89     {
90 jfenwick 3251 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
91 gross 2748
92 jfenwick 3187 A_p = getSampleDataRO(A, e);
93     B_p = getSampleDataRO(B, e);
94     C_p = getSampleDataRO(C, e);
95     D_p = getSampleDataRO(D, e);
96     X_p = getSampleDataRO(X, e);
97     Y_p = getSampleDataRO(Y, e);
98 jfenwick 3205 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
99 jfenwick 3187 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 gross 2748
106 jfenwick 3981 /************************************************************************************/
107 jfenwick 3187 /* process A: */
108 jfenwick 3981 /************************************************************************************/
109 jfenwick 3187 A_p = getSampleDataRO(A, e);
110     if (NULL != A_p)
111     {
112     add_EM_S = TRUE;
113     if (extendedA)
114     {
115 jfenwick 3205 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)]);
116 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
117 jfenwick 3187 {
118 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
119 jfenwick 3187 {
120     for (k = 0; k < p.numEqu; k++)
121     {
122     for (m = 0; m < p.numComp; m++)
123     {
124     rtmp = 0;
125 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
126 jfenwick 3187 {
127     rtmp +=
128 jfenwick 3204 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129 jfenwick 3187 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
130 jfenwick 3204 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
131     DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
132 jfenwick 3187 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
133 jfenwick 3204 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
134     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135 jfenwick 3187 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
136 jfenwick 3204 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
137     DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
138 jfenwick 3187 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
139 jfenwick 3204 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
140 jfenwick 3187 }
141 jfenwick 3224 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
142 jfenwick 3187 }
143     }
144     }
145     }
146 jfenwick 3224 }
147     else
148 jfenwick 3187 {
149 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
150 jfenwick 3187 {
151 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
152 jfenwick 3187 {
153     rtmp00 = 0;
154     rtmp01 = 0;
155     rtmp10 = 0;
156     rtmp11 = 0;
157 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
158 jfenwick 3187 {
159 jfenwick 3204 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
160     rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
161     rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
162     rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
163     rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
164     rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
165 jfenwick 3187 }
166     for (k = 0; k < p.numEqu; k++)
167     {
168     for (m = 0; m < p.numComp; m++)
169     {
170 jfenwick 3204 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
171 jfenwick 3187 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
172     + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
173     + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
174     + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
175     }
176     }
177     }
178     }
179     }
180     }
181 jfenwick 3981 /************************************************************************************/
182 jfenwick 3187 /* process B: */
183 jfenwick 3981 /************************************************************************************/
184 jfenwick 3187 B_p = getSampleDataRO(B, e);
185     if (NULL != B_p)
186     {
187     add_EM_S = TRUE;
188     if (extendedB)
189     {
190 jfenwick 3205 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]);
191 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
192 jfenwick 3187 {
193 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
194 jfenwick 3187 {
195     for (k = 0; k < p.numEqu; k++)
196     {
197     for (m = 0; m < p.numComp; m++)
198     {
199     rtmp = 0;
200 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
201 jfenwick 3187 {
202 jfenwick 3204 rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
203     (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
204 jfenwick 3187 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
205 jfenwick 3204 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
206 jfenwick 3187 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
207     }
208 jfenwick 3224 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
209 jfenwick 3187 }
210     }
211     }
212     }
213 jfenwick 3224 }
214     else
215 jfenwick 3187 {
216 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
217 jfenwick 3187 {
218 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
219 jfenwick 3187 {
220     rtmp0 = 0;
221     rtmp1 = 0;
222 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
223 jfenwick 3187 {
224 jfenwick 3204 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
225     rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
226     rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
227 jfenwick 3187 }
228     for (k = 0; k < p.numEqu; k++)
229     {
230     for (m = 0; m < p.numComp; m++)
231     {
232 jfenwick 3204 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
233 jfenwick 3187 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
234     rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
235     }
236     }
237     }
238     }
239     }
240     }
241 jfenwick 3981 /************************************************************************************/
242 jfenwick 3187 /* process C: */
243 jfenwick 3981 /************************************************************************************/
244 jfenwick 3187 C_p = getSampleDataRO(C, e);
245     if (NULL != C_p)
246     {
247     add_EM_S = TRUE;
248     if (extendedC)
249     {
250 jfenwick 3205 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]);
251 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
252 jfenwick 3187 {
253 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
254 jfenwick 3187 {
255     for (k = 0; k < p.numEqu; k++)
256     {
257     for (m = 0; m < p.numComp; m++)
258     {
259     rtmp = 0;
260 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
261 jfenwick 3187 {
262 jfenwick 3204 rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
263 jfenwick 3187 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
264 jfenwick 3204 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
265 jfenwick 3187 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
266 jfenwick 3204 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
267 jfenwick 3187 }
268 jfenwick 3224 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
269 jfenwick 3187 }
270     }
271     }
272     }
273 jfenwick 3224 }
274     else
275 jfenwick 3187 {
276 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
277 jfenwick 3187 {
278 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
279 jfenwick 3187 {
280     rtmp0 = 0;
281     rtmp1 = 0;
282 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
283 jfenwick 3187 {
284 jfenwick 3204 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
285     rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
286     rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
287 jfenwick 3187 }
288     for (k = 0; k < p.numEqu; k++)
289     {
290     for (m = 0; m < p.numComp; m++)
291     {
292 jfenwick 3204 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
293 jfenwick 3187 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
294     rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
295     }
296     }
297     }
298     }
299     }
300     }
301 jfenwick 3981 /*********************************************************************************** */
302 jfenwick 3187 /* process D */
303 jfenwick 3981 /************************************************************************************/
304 jfenwick 3187 D_p = getSampleDataRO(D, e);
305     if (NULL != D_p)
306     {
307     add_EM_S = TRUE;
308     if (extendedD)
309     {
310 jfenwick 3205 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]);
311 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
312 jfenwick 3187 {
313 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
314 jfenwick 3187 {
315     for (k = 0; k < p.numEqu; k++)
316     {
317     for (m = 0; m < p.numComp; m++)
318     {
319     rtmp = 0;
320 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
321 jfenwick 3187 {
322     rtmp +=
323 jfenwick 3204 vol * S[INDEX2(s, q, p.numShapes)] *
324 jfenwick 3187 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
325 jfenwick 3204 S[INDEX2(r, q, p.numShapes)];
326 jfenwick 3187 }
327 jfenwick 3224 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
328 jfenwick 3187 }
329     }
330     }
331     }
332 jfenwick 3224 }
333     else
334 jfenwick 3187 {
335 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
336 jfenwick 3187 {
337 jfenwick 3204 for (r = 0; r < p.numShapes; r++)
338 jfenwick 3187 {
339     rtmp = 0;
340 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
341 jfenwick 3224 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
342 jfenwick 3187 for (k = 0; k < p.numEqu; k++)
343     {
344     for (m = 0; m < p.numComp; m++)
345     {
346 jfenwick 3204 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
347 jfenwick 3187 rtmp * D_p[INDEX2(k, m, p.numEqu)];
348     }
349     }
350     }
351     }
352     }
353     }
354 jfenwick 3981 /************************************************************************************/
355 jfenwick 3187 /* process X: */
356 jfenwick 3981 /************************************************************************************/
357 jfenwick 3187 X_p = getSampleDataRO(X, e);
358     if (NULL != X_p)
359     {
360     add_EM_F = TRUE;
361     if (extendedX)
362     {
363 jfenwick 3205 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]);
364 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
365 jfenwick 3187 {
366     for (k = 0; k < p.numEqu; k++)
367     {
368     rtmp = 0;
369 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
370 jfenwick 3187 {
371     rtmp +=
372 jfenwick 3204 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
373 jfenwick 3187 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
374 jfenwick 3204 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
375 jfenwick 3187 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
376     }
377     EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
378     }
379     }
380 jfenwick 3224 }
381     else
382 jfenwick 3187 {
383 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
384 jfenwick 3187 {
385     rtmp0 = 0;
386     rtmp1 = 0;
387 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
388 jfenwick 3187 {
389 jfenwick 3204 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
390     rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
391 jfenwick 3187 }
392     for (k = 0; k < p.numEqu; k++)
393     EM_F[INDEX2(k, s, p.numEqu)] +=
394     rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)];
395     }
396     }
397     }
398 jfenwick 3981 /************************************************************************************/
399 jfenwick 3187 /* process Y: */
400 jfenwick 3981 /************************************************************************************/
401 jfenwick 3187 Y_p = getSampleDataRO(Y, e);
402     if (NULL != Y_p)
403     {
404     add_EM_F = TRUE;
405     if (extendedY)
406     {
407 jfenwick 3205 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]);
408 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
409 jfenwick 3187 {
410     for (k = 0; k < p.numEqu; k++)
411     {
412     rtmp = 0;
413 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
414 jfenwick 3224 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
415 jfenwick 3187 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
416     }
417     }
418 jfenwick 3224 }
419     else
420 jfenwick 3187 {
421 jfenwick 3204 for (s = 0; s < p.numShapes; s++)
422 jfenwick 3187 {
423     rtmp = 0;
424 jfenwick 3205 for (q = 0; q < p.numQuad; q++)
425 jfenwick 3204 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
426 jfenwick 3187 for (k = 0; k < p.numEqu; k++)
427     EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
428     }
429     }
430     }
431 jfenwick 3981 /*********************************************************************************************************************/
432 jfenwick 3187 /* add the element matrices onto the matrix and right hand side */
433 jfenwick 3981 /*********************************************************************************************************************/
434 jfenwick 3204 for (q = 0; q < p.numShapes; q++)
435 jfenwick 3187 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
436 gross 798
437 jfenwick 3187 if (add_EM_F)
438 jfenwick 3224 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
439 jfenwick 3187 if (add_EM_S)
440 jfenwick 3204 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
441     p.numShapes, row_index, p.numComp, EM_S);
442 jfenwick 3187
443     } /* end color check */
444     } /* end element loop */
445     } /* end color loop */
446    
447     THREAD_MEMFREE(EM_S);
448     THREAD_MEMFREE(EM_F);
449     THREAD_MEMFREE(row_index);
450    
451     } /* end of pointer check */
452     } /* end parallel region */
453 gross 798 }
454 jfenwick 3187
455 gross 798 /*
456     * $Log$
457     */

  ViewVC Help
Powered by ViewVC 1.1.26