/[escript]/trunk/dudley/src/Assemble_PDE_System2_2D.cpp
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_PDE_System2_2D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:mergeinfo /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 /trunk/ripley/test/python/dudley/src/Assemble_PDE_System2_2D.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26