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

  ViewVC Help
Powered by ViewVC 1.1.26