/[escript]/trunk/dudley/src/Assemble_PDE_System2_3D.c
ViewVC logotype

Annotation of /trunk/dudley/src/Assemble_PDE_System2_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 16758 byte(s)
Merging dudley and scons updates from branches

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

  ViewVC Help
Powered by ViewVC 1.1.26