/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_PDE_Single_3D.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_Single_3D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3187 - (show annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 8 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_Single2_3D.c
File MIME type: text/plain
File size: 14943 byte(s)
Enforcing indent -linux -i4 -bl -bli0 -l120

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * 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
14 /**************************************************************/
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 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 = 3 x 3 */
27 /* B = 3 */
28 /* C = 3 */
29 /* D = scalar */
30 /* X = 3 */
31 /* Y = scalar */
32
33 /**************************************************************/
34
35 #include "Assemble.h"
36 #include "Util.h"
37 #ifdef _OPENMP
38 #include <omp.h>
39 #endif
40
41 /**************************************************************/
42
43 void Dudley_Assemble_PDE_Single2_3D(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
49 #define DIM 3
50 index_t color;
51 dim_t e;
52 __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 double *EM_S, *EM_F, *DSDX;
54 index_t *row_index;
55 register dim_t q, s, r;
56 register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
57 bool_t add_EM_F, add_EM_S;
58
59 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 dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal;
68 dim_t len_EM_F = p.row_numShapesTotal;
69
70 #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, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
71 {
72 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
76 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
77 {
78
79 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
88 A_p = getSampleDataRO(A, e);
89 B_p = getSampleDataRO(B, e);
90 C_p = getSampleDataRO(C, e);
91 D_p = getSampleDataRO(D, e);
92 X_p = getSampleDataRO(X, e);
93 Y_p = getSampleDataRO(Y, e);
94
95 // Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
96 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
97 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 /**************************************************************/
106 /* process A: */
107 /**************************************************************/
108 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 for (r = 0; r < p.col_numShapes; r++)
117 {
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 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
126 A_q[INDEX3(0, 1, q, DIM, DIM)] *
127 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
128 DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
129 A_q[INDEX3(0, 2, q, DIM, DIM)] *
130 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
131 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
132 A_q[INDEX3(1, 0, q, DIM, DIM)] *
133 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
134 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
135 A_q[INDEX3(1, 1, q, DIM, DIM)] *
136 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
137 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
138 A_q[INDEX3(1, 2, q, DIM, DIM)] *
139 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)] +
140 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
141 A_q[INDEX3(2, 0, q, DIM, DIM)] *
142 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
143 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
144 A_q[INDEX3(2, 1, q, DIM, DIM)] *
145 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
146 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
147 A_q[INDEX3(2, 2, q, DIM, DIM)] *
148 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
149 }
150 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
151 }
152 }
153 } else
154 {
155 for (s = 0; s < p.row_numShapes; s++)
156 {
157 for (r = 0; r < p.col_numShapes; r++)
158 {
159 rtmp00 = 0;
160 rtmp01 = 0;
161 rtmp02 = 0;
162 rtmp10 = 0;
163 rtmp11 = 0;
164 rtmp12 = 0;
165 rtmp20 = 0;
166 rtmp21 = 0;
167 rtmp22 = 0;
168 for (q = 0; q < p.numQuadTotal; q++)
169 {
170
171 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
172 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
173 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
174 rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
175
176 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
177 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
178 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
179 rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
180
181 rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
182 rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
183 rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
184 rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
185 }
186 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
187 rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
188 rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
189 rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
190 rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
191 rtmp22 * A_p[INDEX2(2, 2, DIM)];
192 }
193 }
194 }
195 }
196 /**************************************************************/
197 /* process B: */
198 /**************************************************************/
199 if (NULL != B_p)
200 {
201 add_EM_S = TRUE;
202 if (extendedB)
203 {
204 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
205 for (s = 0; s < p.row_numShapes; s++)
206 {
207 for (r = 0; r < p.col_numShapes; r++)
208 {
209 rtmp = 0;
210 for (q = 0; q < p.numQuadTotal; q++)
211 {
212 rtmp += vol * S[INDEX2(r, q, p.row_numShapes)] *
213 (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
214 B_q[INDEX2(0, q, DIM)] +
215 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
216 B_q[INDEX2(1, q, DIM)] +
217 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
218 B_q[INDEX2(2, q, DIM)]);
219 }
220 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
221 }
222 }
223 } else
224 {
225 for (s = 0; s < p.row_numShapes; s++)
226 {
227 for (r = 0; r < p.col_numShapes; r++)
228 {
229 rtmp0 = 0;
230 rtmp1 = 0;
231 rtmp2 = 0;
232 for (q = 0; q < p.numQuadTotal; q++)
233 {
234 rtmp = vol * S[INDEX2(r, q, p.row_numShapes)];
235 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
236 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
237 rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
238 }
239 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
240 rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
241 }
242 }
243 }
244 }
245 /**************************************************************/
246 /* process C: */
247 /**************************************************************/
248 if (NULL != C_p)
249 {
250 add_EM_S = TRUE;
251 if (extendedC)
252 {
253 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
254 for (s = 0; s < p.row_numShapes; s++)
255 {
256 for (r = 0; r < p.col_numShapes; r++)
257 {
258 rtmp = 0;
259 for (q = 0; q < p.numQuadTotal; q++)
260 {
261 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] *
262 (C_q[INDEX2(0, q, DIM)] *
263 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)] +
264 C_q[INDEX2(1, q, DIM)] *
265 DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)] +
266 C_q[INDEX2(2, q, DIM)] *
267 DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)]);
268 }
269 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
270 }
271 }
272 } else
273 {
274 for (s = 0; s < p.row_numShapes; s++)
275 {
276 for (r = 0; r < p.col_numShapes; r++)
277 {
278 rtmp0 = 0;
279 rtmp1 = 0;
280 rtmp2 = 0;
281 for (q = 0; q < p.numQuadTotal; q++)
282 {
283 rtmp = vol * S[INDEX2(s, q, p.row_numShapes)];
284 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
285 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.row_numShapesTotal, DIM)];
286 rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.row_numShapesTotal, DIM)];
287 }
288 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
289 rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
290 }
291 }
292 }
293 }
294 /************************************************************* */
295 /* process D */
296 /**************************************************************/
297 if (NULL != D_p)
298 {
299 add_EM_S = TRUE;
300 if (extendedD)
301 {
302 D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
303 for (s = 0; s < p.row_numShapes; s++)
304 {
305 for (r = 0; r < p.col_numShapes; r++)
306 {
307 rtmp = 0;
308 for (q = 0; q < p.numQuadTotal; q++)
309 rtmp +=
310 vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
311 S[INDEX2(r, q, p.row_numShapes)];
312 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
313 }
314 }
315 } else
316 {
317 for (s = 0; s < p.row_numShapes; s++)
318 {
319 for (r = 0; r < p.col_numShapes; r++)
320 {
321 rtmp = 0;
322 for (q = 0; q < p.numQuadTotal; q++)
323 rtmp +=
324 vol * S[INDEX2(s, q, p.row_numShapes)] *
325 S[INDEX2(r, q, p.row_numShapes)];
326 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
327 rtmp * D_p[0];
328 }
329 }
330 }
331 }
332 /**************************************************************/
333 /* process X: */
334 /**************************************************************/
335 if (NULL != X_p)
336 {
337 add_EM_F = TRUE;
338 if (extendedX)
339 {
340 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
341 for (s = 0; s < p.row_numShapes; s++)
342 {
343 rtmp = 0;
344 for (q = 0; q < p.numQuadTotal; q++)
345 {
346 rtmp +=
347 vol * (DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
348 X_q[INDEX2(0, q, DIM)] +
349 DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)] *
350 X_q[INDEX2(1, q, DIM)] +
351 DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)] *
352 X_q[INDEX2(2, q, DIM)]);
353 }
354 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
355 }
356 } else
357 {
358 for (s = 0; s < p.row_numShapes; s++)
359 {
360 rtmp0 = 0;
361 rtmp1 = 0;
362 rtmp2 = 0;
363 for (q = 0; q < p.numQuadTotal; q++)
364 {
365 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
366 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.row_numShapesTotal, DIM)];
367 rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.row_numShapesTotal, DIM)];
368 }
369 EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
370 }
371 }
372 }
373 /**************************************************************/
374 /* process Y: */
375 /**************************************************************/
376 if (NULL != Y_p)
377 {
378 add_EM_F = TRUE;
379 if (extendedY)
380 {
381 Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
382 for (s = 0; s < p.row_numShapes; s++)
383 {
384 rtmp = 0;
385 for (q = 0; q < p.numQuadTotal; q++)
386 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
387 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
388 }
389 } else
390 {
391 for (s = 0; s < p.row_numShapes; s++)
392 {
393 rtmp = 0;
394 for (q = 0; q < p.numQuadTotal; q++)
395 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
396 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
397 }
398 }
399 }
400 /***********************************************************************************************/
401 /* add the element matrices onto the matrix and right hand side */
402 /***********************************************************************************************/
403
404 for (q = 0; q < p.row_numShapesTotal; q++)
405 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
406
407 if (add_EM_F)
408 Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
409 p.row_DOF_UpperBound);
410 if (add_EM_S)
411 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
412 p.col_numShapesTotal, row_index, p.numComp, EM_S);
413
414 } /* end color check */
415 } /* end element loop */
416 } /* end color loop */
417
418 THREAD_MEMFREE(EM_S);
419 THREAD_MEMFREE(EM_F);
420 THREAD_MEMFREE(row_index);
421
422 } /* end of pointer check */
423 } /* end parallel region */
424 }

  ViewVC Help
Powered by ViewVC 1.1.26