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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26