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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5933 - (show annotations)
Wed Feb 17 23:53:30 2016 UTC (21 months, 4 weeks ago) by caltinay
File size: 14617 byte(s)
sync with trunk.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 /************************************************************************************/
18
19 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
20 /* the shape functions for test and solution must be identical */
21
22 /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
23
24 /* in a 3D domain. The shape functions for test and solution must be identical */
25 /* and row_NS == row_NN */
26
27 /* Shape of the coefficients: */
28
29 /* A = 3 x 3 */
30 /* B = 3 */
31 /* C = 3 */
32 /* D = scalar */
33 /* X = 3 */
34 /* Y = scalar */
35
36 /************************************************************************************/
37
38 #define ESNEEDPYTHON
39 #include "esysUtils/first.h"
40
41 #include "Assemble.h"
42 #include "Util.h"
43 #ifdef _OPENMP
44 #include <omp.h>
45 #endif
46
47 /************************************************************************************/
48
49 void Dudley_Assemble_PDE_Single2_3D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
50 escript::ASM_ptr mat, escript::Data* F,
51 const escript::Data* A, const escript::Data* B, const escript::Data* C, const escript::Data* D,
52 const escript::Data* X, const escript::Data* Y)
53 {
54
55 #define DIM 3
56 index_t color;
57 dim_t e;
58 __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;
59 double *EM_S, *EM_F, *DSDX;
60 index_t *row_index;
61 register dim_t q, s, r;
62 register double rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2;
63 bool add_EM_F, add_EM_S;
64
65 bool extendedA = isExpanded(A);
66 bool extendedB = isExpanded(B);
67 bool extendedC = isExpanded(C);
68 bool extendedD = isExpanded(D);
69 bool extendedX = isExpanded(X);
70 bool extendedY = isExpanded(Y);
71 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
72 const double *S = p.shapeFns;
73 dim_t len_EM_S = p.numShapes * p.numShapes;
74 dim_t len_EM_F = p.numShapes;
75
76 #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,rtmp, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22, rtmp0, rtmp1, rtmp2,add_EM_F, add_EM_S)
77 {
78 EM_S = new double[len_EM_S];
79 EM_F = new double[len_EM_F];
80 row_index = new index_t[p.numShapes];
81
82 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
83 {
84
85 for (color = elements->minColor; color <= elements->maxColor; color++)
86 {
87 /* open loop over all elements: */
88 #pragma omp for private(e) schedule(static)
89 for (e = 0; e < elements->numElements; e++)
90 {
91 if (elements->Color[e] == color)
92 {
93 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
94
95 A_p = getSampleDataRO(A, e);
96 B_p = getSampleDataRO(B, e);
97 C_p = getSampleDataRO(C, e);
98 D_p = getSampleDataRO(D, e);
99 X_p = getSampleDataRO(X, e);
100 Y_p = getSampleDataRO(Y, e);
101
102 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
103 for (q = 0; q < len_EM_S; ++q)
104 EM_S[q] = 0;
105 for (q = 0; q < len_EM_F; ++q)
106 EM_F[q] = 0;
107 add_EM_F = FALSE;
108 add_EM_S = FALSE;
109
110 /************************************************************************************/
111 /* process A: */
112 /************************************************************************************/
113 if (NULL != A_p)
114 {
115 add_EM_S = TRUE;
116 if (extendedA)
117 {
118 A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)]);
119 for (s = 0; s < p.numShapes; s++)
120 {
121 for (r = 0; r < p.numShapes; r++)
122 {
123 rtmp = 0;
124 for (q = 0; q < p.numQuad; q++)
125 {
126 rtmp +=
127 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
128 A_q[INDEX3(0, 0, q, DIM, DIM)] *
129 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
130 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
131 A_q[INDEX3(0, 1, q, DIM, DIM)] *
132 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
133 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
134 A_q[INDEX3(0, 2, q, DIM, DIM)] *
135 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
136 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
137 A_q[INDEX3(1, 0, q, DIM, DIM)] *
138 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
139 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
140 A_q[INDEX3(1, 1, q, DIM, DIM)] *
141 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
142 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
143 A_q[INDEX3(1, 2, q, DIM, DIM)] *
144 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
145 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
146 A_q[INDEX3(2, 0, q, DIM, DIM)] *
147 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
148 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
149 A_q[INDEX3(2, 1, q, DIM, DIM)] *
150 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
151 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
152 A_q[INDEX3(2, 2, q, DIM, DIM)] *
153 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
154 }
155 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
156 }
157 }
158 }
159 else
160 {
161 for (s = 0; s < p.numShapes; s++)
162 {
163 for (r = 0; r < p.numShapes; r++)
164 {
165 rtmp00 = 0;
166 rtmp01 = 0;
167 rtmp02 = 0;
168 rtmp10 = 0;
169 rtmp11 = 0;
170 rtmp12 = 0;
171 rtmp20 = 0;
172 rtmp21 = 0;
173 rtmp22 = 0;
174 for (q = 0; q < p.numQuad; q++)
175 {
176
177 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
178 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
179 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
180 rtmp02 += rtmp0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
181
182 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
183 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
184 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
185 rtmp12 += rtmp1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
186
187 rtmp2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
188 rtmp20 += rtmp2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
189 rtmp21 += rtmp2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
190 rtmp22 += rtmp2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
191 }
192 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
193 rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
194 rtmp02 * A_p[INDEX2(0, 2, DIM)] + rtmp10 * A_p[INDEX2(1, 0, DIM)] +
195 rtmp11 * A_p[INDEX2(1, 1, DIM)] + rtmp12 * A_p[INDEX2(1, 2, DIM)] +
196 rtmp20 * A_p[INDEX2(2, 0, DIM)] + rtmp21 * A_p[INDEX2(2, 1, DIM)] +
197 rtmp22 * A_p[INDEX2(2, 2, DIM)];
198 }
199 }
200 }
201 }
202 /************************************************************************************/
203 /* process B: */
204 /************************************************************************************/
205 if (NULL != B_p)
206 {
207 add_EM_S = TRUE;
208 if (extendedB)
209 {
210 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
211 for (s = 0; s < p.numShapes; s++)
212 {
213 for (r = 0; r < p.numShapes; r++)
214 {
215 rtmp = 0;
216 for (q = 0; q < p.numQuad; q++)
217 {
218 rtmp += vol * S[INDEX2(r, q, p.numShapes)] *
219 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
220 B_q[INDEX2(0, q, DIM)] +
221 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
222 B_q[INDEX2(1, q, DIM)] +
223 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * B_q[INDEX2(2, q, DIM)]);
224 }
225 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
226 }
227 }
228 }
229 else
230 {
231 for (s = 0; s < p.numShapes; s++)
232 {
233 for (r = 0; r < p.numShapes; r++)
234 {
235 rtmp0 = 0;
236 rtmp1 = 0;
237 rtmp2 = 0;
238 for (q = 0; q < p.numQuad; q++)
239 {
240 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
241 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
242 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
243 rtmp2 += rtmp * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
244 }
245 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
246 rtmp0 * B_p[0] + rtmp1 * B_p[1] + rtmp2 * B_p[2];
247 }
248 }
249 }
250 }
251 /************************************************************************************/
252 /* process C: */
253 /************************************************************************************/
254 if (NULL != C_p)
255 {
256 add_EM_S = TRUE;
257 if (extendedC)
258 {
259 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
260 for (s = 0; s < p.numShapes; s++)
261 {
262 for (r = 0; r < p.numShapes; r++)
263 {
264 rtmp = 0;
265 for (q = 0; q < p.numQuad; q++)
266 {
267 rtmp += vol * S[INDEX2(s, q, p.numShapes)] *
268 (C_q[INDEX2(0, q, DIM)] *
269 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
270 C_q[INDEX2(1, q, DIM)] *
271 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
272 C_q[INDEX2(2, q, DIM)] * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
273 }
274 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
275 }
276 }
277 }
278 else
279 {
280 for (s = 0; s < p.numShapes; s++)
281 {
282 for (r = 0; r < p.numShapes; r++)
283 {
284 rtmp0 = 0;
285 rtmp1 = 0;
286 rtmp2 = 0;
287 for (q = 0; q < p.numQuad; q++)
288 {
289 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
290 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
291 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
292 rtmp2 += rtmp * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
293 }
294 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
295 rtmp0 * C_p[0] + rtmp1 * C_p[1] + rtmp2 * C_p[2];
296 }
297 }
298 }
299 }
300 /*********************************************************************************** */
301 /* process D */
302 /************************************************************************************/
303 if (NULL != D_p)
304 {
305 add_EM_S = TRUE;
306 if (extendedD)
307 {
308 D_q = &(D_p[INDEX2(0, 0, p.numQuad)]);
309 for (s = 0; s < p.numShapes; s++)
310 {
311 for (r = 0; r < p.numShapes; r++)
312 {
313 rtmp = 0;
314 for (q = 0; q < p.numQuad; q++)
315 rtmp +=
316 vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
317 S[INDEX2(r, q, p.numShapes)];
318 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
319 }
320 }
321 }
322 else
323 {
324 for (s = 0; s < p.numShapes; s++)
325 {
326 for (r = 0; r < p.numShapes; r++)
327 {
328 rtmp = 0;
329 for (q = 0; q < p.numQuad; q++)
330 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
331 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp * D_p[0];
332 }
333 }
334 }
335 }
336 /************************************************************************************/
337 /* process X: */
338 /************************************************************************************/
339 if (NULL != X_p)
340 {
341 add_EM_F = TRUE;
342 if (extendedX)
343 {
344 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
345 for (s = 0; s < p.numShapes; s++)
346 {
347 rtmp = 0;
348 for (q = 0; q < p.numQuad; q++)
349 {
350 rtmp +=
351 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
352 X_q[INDEX2(0, q, DIM)] +
353 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
354 X_q[INDEX2(1, q, DIM)] +
355 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * X_q[INDEX2(2, q, DIM)]);
356 }
357 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
358 }
359 }
360 else
361 {
362 for (s = 0; s < p.numShapes; s++)
363 {
364 rtmp0 = 0;
365 rtmp1 = 0;
366 rtmp2 = 0;
367 for (q = 0; q < p.numQuad; q++)
368 {
369 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
370 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
371 rtmp2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
372 }
373 EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1] + rtmp2 * X_p[2];
374 }
375 }
376 }
377 /************************************************************************************/
378 /* process Y: */
379 /************************************************************************************/
380 if (NULL != Y_p)
381 {
382 add_EM_F = TRUE;
383 if (extendedY)
384 {
385 Y_q = &(Y_p[INDEX2(0, 0, p.numQuad)]);
386 for (s = 0; s < p.numShapes; s++)
387 {
388 rtmp = 0;
389 for (q = 0; q < p.numQuad; q++)
390 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
391 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
392 }
393 }
394 else
395 {
396 for (s = 0; s < p.numShapes; s++)
397 {
398 rtmp = 0;
399 for (q = 0; q < p.numQuad; q++)
400 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
401 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
402 }
403 }
404 }
405 /*********************************************************************************************************************/
406 /* add the element matrices onto the matrix and right hand side */
407 /*********************************************************************************************************************/
408
409 for (q = 0; q < p.numShapes; q++)
410 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
411
412 if (add_EM_F)
413 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
414 if (add_EM_S)
415 Dudley_Assemble_addToSystemMatrix(mat, p.numShapes, row_index, p.numEqu,
416 p.numShapes, row_index, p.numComp, EM_S);
417
418 } /* end color check */
419 } /* end element loop */
420 } /* end color loop */
421
422 delete[] EM_S;
423 delete[] EM_F;
424 delete[] row_index;
425
426 } /* end of pointer check */
427 } /* end parallel region */
428 }

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/Assemble_PDE_Single2_3D.cpp:5567-5588 /branches/lapack2681/finley/src/Assemble_PDE_Single2_3D.cpp:2682-2741 /branches/pasowrap/dudley/src/Assemble_PDE_Single2_3D.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Assemble_PDE_Single2_3D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_Single2_3D.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Assemble_PDE_Single2_3D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_Single2_3D.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Assemble_PDE_Single2_3D.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_PDE_Single2_3D.cpp:3517-3974 /release/3.0/finley/src/Assemble_PDE_Single2_3D.cpp:2591-2601 /release/4.0/dudley/src/Assemble_PDE_Single2_3D.cpp:5380-5406 /trunk/dudley/src/Assemble_PDE_Single2_3D.cpp:4257-4344,5898-5932 /trunk/ripley/test/python/dudley/src/Assemble_PDE_Single2_3D.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26