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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 3 months ago) by caltinay
File size: 12987 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 1D 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 1 x p.numComp x 1 */
30 /* B = 1 x numEqu x p.numComp */
31 /* C = p.numEqu x 1 x p.numComp */
32 /* D = p.numEqu x p.numComp */
33 /* X = p.numEqu x 1 */
34 /* Y = p.numEqu */
35
36 /************************************************************************************/
37
38 #include "Assemble.h"
39 #include "Util.h"
40 #ifdef _OPENMP
41 #include <omp.h>
42 #endif
43
44 /************************************************************************************/
45
46 void Dudley_Assemble_PDE_System2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
47 escript::ASM_ptr mat, escriptDataC * F,
48 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
49 escriptDataC * X, escriptDataC * Y)
50 {
51
52 #define DIM 1
53 index_t color;
54 dim_t e;
55 __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;
56 double *EM_S, *EM_F, *DSDX;
57 index_t *row_index;
58 dim_t q, s, r, k, m;
59 double rtmp;
60 bool add_EM_F, add_EM_S;
61
62 bool extendedA = isExpanded(A);
63 bool extendedB = isExpanded(B);
64 bool extendedC = isExpanded(C);
65 bool extendedD = isExpanded(D);
66 bool extendedX = isExpanded(X);
67 bool extendedY = isExpanded(Y);
68 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
69 double *S = p.row_jac->BasisFunctions->S;
70 dim_t len_EM_S = p.row_numShapes * p.row_numShapes * p.numEqu * p.numComp;
71 dim_t len_EM_F = p.row_numShapes * p.numEqu;
72
73 #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,k,m,rtmp,add_EM_F, add_EM_S)
74 {
75 EM_S = new double[len_EM_S];
76 EM_F = new double[len_EM_F];
77 row_index = new index_t[p.row_numShapes];
78
79 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
80 {
81
82 for (color = elements->minColor; color <= elements->maxColor; color++)
83 {
84 /* open loop over all elements: */
85 #pragma omp for private(e) schedule(static)
86 for (e = 0; e < elements->numElements; e++)
87 {
88 if (elements->Color[e] == color)
89 {
90
91 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 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
98 Vol = &(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal, 1)]);
99 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapes, DIM, p.numQuadTotal, 1)]);
100 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
107 /************************************************************************************/
108 /* process A: */
109 /************************************************************************************/
110
111 if (NULL != A_p)
112 {
113 add_EM_S = TRUE;
114 if (extendedA)
115 {
116 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
117 for (s = 0; s < p.row_numShapes; s++)
118 {
119 for (r = 0; r < p.row_numShapes; r++)
120 {
121 for (k = 0; k < p.numEqu; k++)
122 {
123 for (m = 0; m < p.numComp; m++)
124 {
125 rtmp = 0.;
126 for (q = 0; q < p.numQuadTotal; q++)
127 {
128 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
129 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] *
130 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
131 }
132 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
133 }
134 }
135 }
136 }
137 }
138 else
139 {
140 for (s = 0; s < p.row_numShapes; s++)
141 {
142 for (r = 0; r < p.row_numShapes; r++)
143 {
144 rtmp = 0;
145 for (q = 0; q < p.numQuadTotal; q++)
146 rtmp +=
147 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
148 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
149 for (k = 0; k < p.numEqu; k++)
150 {
151 for (m = 0; m < p.numComp; m++)
152 {
153 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
154 rtmp * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)];
155 }
156 }
157 }
158 }
159 }
160 }
161 /************************************************************************************/
162 /* process B: */
163 /************************************************************************************/
164
165 if (NULL != B_p)
166 {
167 add_EM_S = TRUE;
168 if (extendedB)
169 {
170 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
171 for (s = 0; s < p.row_numShapes; s++)
172 {
173 for (r = 0; r < p.row_numShapes; r++)
174 {
175 for (k = 0; k < p.numEqu; k++)
176 {
177 for (m = 0; m < p.numComp; m++)
178 {
179 rtmp = 0.;
180 for (q = 0; q < p.numQuadTotal; q++)
181 {
182 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
183 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] *
184 S[INDEX2(r, q, p.row_numShapes)];
185 }
186 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
187 }
188 }
189 }
190 }
191 }
192 else
193 {
194 for (s = 0; s < p.row_numShapes; s++)
195 {
196 for (r = 0; r < p.row_numShapes; r++)
197 {
198 rtmp = 0;
199 for (q = 0; q < p.numQuadTotal; q++)
200 rtmp +=
201 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
202 S[INDEX2(r, q, p.row_numShapes)];
203 for (k = 0; k < p.numEqu; k++)
204 {
205 for (m = 0; m < p.numComp; m++)
206 {
207 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
208 rtmp * B_p[INDEX3(k, 0, m, p.numEqu, DIM)];
209 }
210 }
211 }
212 }
213 }
214 }
215 /************************************************************************************/
216 /* process C: */
217 /************************************************************************************/
218
219 if (NULL != C_p)
220 {
221 add_EM_S = TRUE;
222 if (extendedC)
223 {
224 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
225 for (s = 0; s < p.row_numShapes; s++)
226 {
227 for (r = 0; r < p.row_numShapes; r++)
228 {
229 for (k = 0; k < p.numEqu; k++)
230 {
231 for (m = 0; m < p.numComp; m++)
232 {
233 rtmp = 0;
234 for (q = 0; q < p.numQuadTotal; q++)
235 {
236 rtmp +=
237 vol * S[INDEX2(s, q, p.row_numShapes)] *
238 C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
239 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
240 }
241 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
242 }
243 }
244 }
245 }
246 }
247 else
248 {
249 for (s = 0; s < p.row_numShapes; s++)
250 {
251 for (r = 0; r < p.row_numShapes; r++)
252 {
253 rtmp = 0;
254 for (q = 0; q < p.numQuadTotal; q++)
255 rtmp +=
256 vol * S[INDEX2(s, q, p.row_numShapes)] *
257 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
258 for (k = 0; k < p.numEqu; k++)
259 {
260 for (m = 0; m < p.numComp; m++)
261 {
262 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
263 rtmp * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)];
264 }
265 }
266 }
267 }
268 }
269 }
270 /*********************************************************************************** */
271 /* process D */
272 /************************************************************************************/
273
274 if (NULL != D_p)
275 {
276 add_EM_S = TRUE;
277 if (extendedD)
278 {
279 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
280 for (s = 0; s < p.row_numShapes; s++)
281 {
282 for (r = 0; r < p.row_numShapes; r++)
283 {
284 for (k = 0; k < p.numEqu; k++)
285 {
286 for (m = 0; m < p.numComp; m++)
287 {
288 rtmp = 0;
289 for (q = 0; q < p.numQuadTotal; q++)
290 {
291 rtmp +=
292 vol * S[INDEX2(s, q, p.row_numShapes)] *
293 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
294 S[INDEX2(r, q, p.row_numShapes)];
295 }
296 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
297
298 }
299 }
300 }
301 }
302 }
303 else
304 {
305 for (s = 0; s < p.row_numShapes; s++)
306 {
307 for (r = 0; r < p.row_numShapes; r++)
308 {
309 rtmp = 0;
310 for (q = 0; q < p.numQuadTotal; q++)
311 rtmp +=
312 vol * S[INDEX2(s, q, p.row_numShapes)] *
313 S[INDEX2(r, q, p.row_numShapes)];
314 for (k = 0; k < p.numEqu; k++)
315 {
316 for (m = 0; m < p.numComp; m++)
317 {
318 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
319 rtmp * D_p[INDEX2(k, m, p.numEqu)];
320 }
321 }
322 }
323 }
324 }
325 }
326 /************************************************************************************/
327 /* process X: */
328 /************************************************************************************/
329
330 if (NULL != X_p)
331 {
332 add_EM_F = TRUE;
333 if (extendedX)
334 {
335 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
336 for (s = 0; s < p.row_numShapes; s++)
337 {
338 for (k = 0; k < p.numEqu; k++)
339 {
340 rtmp = 0;
341 for (q = 0; q < p.numQuadTotal; q++)
342 rtmp +=
343 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
344 X_q[INDEX3(k, 0, q, p.numEqu, DIM)];
345 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
346 }
347 }
348 }
349 else
350 {
351 for (s = 0; s < p.row_numShapes; s++)
352 {
353 rtmp = 0;
354 for (q = 0; q < p.numQuadTotal; q++)
355 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];
356 for (k = 0; k < p.numEqu; k++)
357 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * X_p[INDEX2(k, 0, p.numEqu)];
358 }
359 }
360 }
361 /************************************************************************************/
362 /* process Y: */
363 /************************************************************************************/
364
365 if (NULL != Y_p)
366 {
367 add_EM_F = TRUE;
368 if (extendedY)
369 {
370 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
371 for (s = 0; s < p.row_numShapes; s++)
372 {
373 for (k = 0; k < p.numEqu; k++)
374 {
375 rtmp = 0;
376 for (q = 0; q < p.numQuadTotal; q++)
377 rtmp +=
378 vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
379 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
380 }
381 }
382 }
383 else
384 {
385 for (s = 0; s < p.row_numShapes; s++)
386 {
387 rtmp = 0;
388 for (q = 0; q < p.numQuadTotal; q++)
389 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
390 for (k = 0; k < p.numEqu; k++)
391 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
392 }
393 }
394 }
395 /*********************************************************************************************************************/
396 /* add the element matrices onto the matrix and right hand side */
397 /*********************************************************************************************************************/
398
399 for (q = 0; q < p.row_numShapes; q++)
400 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
401
402 if (add_EM_F)
403 Dudley_Util_AddScatter(p.row_numShapes, row_index, p.numEqu, EM_F, F_p,
404 p.row_DOF_UpperBound);
405 if (add_EM_S)
406 Dudley_Assemble_addToSystemMatrix(mat, p.row_numShapes, row_index, p.numEqu,
407 p.row_numShapes, row_index, p.numComp, EM_S);
408
409 } /* end color check */
410 } /* end element loop */
411 } /* end color loop */
412
413 delete[] EM_S;
414 delete[] EM_F;
415 delete[] row_index;
416
417 } /* end of pointer check */
418 } /* end parallel region */
419 }
420
421 /*
422 * $Log$
423 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26