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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 1 month ago) by caltinay
File size: 14878 byte(s)
sync and fix.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26