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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4346 - (show annotations)
Tue Apr 2 04:46:45 2013 UTC (6 years, 4 months ago) by jfenwick
File size: 12519 byte(s)
Bringing the changes from doubleplusgood branch.
Can't merge directly because svn doesn't transfer changes to renamed files (mutter grumble).
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 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 since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 /************************************************************************************/
17
18 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
19 /* the shape functions for test and solution must be identical */
20
21 /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
22
23 /* in a 2D domain. The shape functions for test and solution must be identical */
24 /* and row_NS == row_NN */
25
26 /* Shape of the coefficients: */
27
28 /* A = 2 x 2 */
29 /* B = 2 */
30 /* C = 2 */
31 /* D = scalar */
32 /* X = 2 */
33 /* Y = scalar */
34
35 /************************************************************************************/
36
37 #include "Assemble.h"
38 #include "Util.h"
39 #ifdef _OPENMP
40 #include <omp.h>
41 #endif
42
43 #include "ShapeTable.h"
44
45 /************************************************************************************/
46
47 void Dudley_Assemble_PDE_Single2_2D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
48 Paso_SystemMatrix * Mat, escriptDataC * F,
49 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
50 escriptDataC * X, escriptDataC * Y)
51 {
52
53 #define DIM 2
54 index_t color;
55 dim_t e;
56 __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;
57 double *EM_S, *EM_F, *DSDX;
58 index_t *row_index;
59 register dim_t q, s, r;
60 register double rtmp00, rtmp01, rtmp10, rtmp11, rtmp, rtmp0, rtmp1;
61 bool_t add_EM_F, add_EM_S;
62
63 bool_t extendedA = isExpanded(A);
64 bool_t extendedB = isExpanded(B);
65 bool_t extendedC = isExpanded(C);
66 bool_t extendedD = isExpanded(D);
67 bool_t extendedX = isExpanded(X);
68 bool_t extendedY = isExpanded(Y);
69 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
70 const double *S = p.shapeFns;
71 dim_t len_EM_S = p.numShapes * p.numShapes;
72 dim_t len_EM_F = p.numShapes;
73
74 #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)
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
100
101 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]);
102 for (q = 0; q < len_EM_S; ++q)
103 EM_S[q] = 0;
104 for (q = 0; q < len_EM_F; ++q)
105 EM_F[q] = 0;
106 add_EM_F = FALSE;
107 add_EM_S = FALSE;
108 /************************************************************************************/
109 /* process A: */
110 /************************************************************************************/
111 if (NULL != A_p)
112 {
113 add_EM_S = TRUE;
114 if (extendedA)
115 {
116 A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)]);
117 for (s = 0; s < p.numShapes; s++)
118 {
119 for (r = 0; r < p.numShapes; r++)
120 {
121 rtmp = 0;
122 for (q = 0; q < p.numQuad; q++)
123 {
124 rtmp +=
125 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
126 A_q[INDEX3(0, 0, q, DIM, DIM)] *
127 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
128 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
129 A_q[INDEX3(0, 1, q, DIM, DIM)] *
130 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
131 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
132 A_q[INDEX3(1, 0, q, DIM, DIM)] *
133 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
134 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
135 A_q[INDEX3(1, 1, q, DIM, DIM)] *
136 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
137
138 }
139 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
140 }
141 }
142 }
143 else
144 {
145 for (s = 0; s < p.numShapes; s++)
146 {
147 for (r = 0; r < p.numShapes; r++)
148 {
149 rtmp00 = 0;
150 rtmp01 = 0;
151 rtmp10 = 0;
152 rtmp11 = 0;
153 for (q = 0; q < p.numQuad; q++)
154 {
155 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
156 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
157 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
158 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
159 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
160 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
161 }
162 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
163 rtmp00 * A_p[INDEX2(0, 0, DIM)] + rtmp01 * A_p[INDEX2(0, 1, DIM)] +
164 rtmp10 * A_p[INDEX2(1, 0, DIM)] + rtmp11 * A_p[INDEX2(1, 1, DIM)];
165 }
166 }
167 }
168 }
169 /************************************************************************************/
170 /* process B: */
171 /************************************************************************************/
172 if (NULL != B_p)
173 {
174 add_EM_S = TRUE;
175 if (extendedB)
176 {
177 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
178 for (s = 0; s < p.numShapes; s++)
179 {
180 for (r = 0; r < p.numShapes; r++)
181 {
182 rtmp = 0.;
183 for (q = 0; q < p.numQuad; q++)
184 {
185 rtmp +=
186 vol * S[INDEX2(r, q, p.numShapes)] *
187 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
188 B_q[INDEX2(0, q, DIM)] +
189 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * B_q[INDEX2(1, q, DIM)]);
190 }
191 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
192 }
193 }
194 }
195 else
196 {
197 for (s = 0; s < p.numShapes; s++)
198 {
199 for (r = 0; r < p.numShapes; r++)
200 {
201 rtmp0 = 0;
202 rtmp1 = 0;
203 for (q = 0; q < p.numQuad; q++)
204 {
205 rtmp = vol * S[INDEX2(r, q, p.numShapes)];
206 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
207 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
208 }
209 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
210 rtmp0 * B_p[0] + rtmp1 * B_p[1];
211 }
212 }
213 }
214 }
215 /************************************************************************************/
216 /* process C: */
217 /************************************************************************************/
218 if (NULL != C_p)
219 {
220 add_EM_S = TRUE;
221 if (extendedC)
222 {
223 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
224 for (s = 0; s < p.numShapes; s++)
225 {
226 for (r = 0; r < p.numShapes; r++)
227 {
228 rtmp = 0;
229 for (q = 0; q < p.numQuad; q++)
230 {
231 rtmp +=
232 vol * S[INDEX2(s, q, p.numShapes)] * (C_q[INDEX2(0, q, DIM)] *
233 DSDX[INDEX3
234 (r, 0, q,
235 p.numShapes,
236 DIM)] + C_q[INDEX2(1, q,
237 DIM)]
238 *
239 DSDX[INDEX3
240 (r, 1, q,
241 p.numShapes, DIM)]);
242 }
243 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
244 }
245 }
246 }
247 else
248 {
249 for (s = 0; s < p.numShapes; s++)
250 {
251 for (r = 0; r < p.numShapes; r++)
252 {
253 rtmp0 = 0;
254 rtmp1 = 0;
255 for (q = 0; q < p.numQuad; q++)
256 {
257 rtmp = vol * S[INDEX2(s, q, p.numShapes)];
258 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
259 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
260 }
261 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
262 rtmp0 * C_p[0] + rtmp1 * C_p[1];
263 }
264 }
265 }
266 }
267 /*********************************************************************************** */
268 /* process D */
269 /************************************************************************************/
270 if (NULL != D_p)
271 {
272 add_EM_S = TRUE;
273 if (extendedD)
274 {
275 D_q = &(D_p[INDEX2(0, 0, p.numQuad)]);
276 for (s = 0; s < p.numShapes; s++)
277 {
278 for (r = 0; r < p.numShapes; r++)
279 {
280 rtmp = 0;
281 for (q = 0; q < p.numQuad; q++)
282 rtmp +=
283 vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
284 S[INDEX2(r, q, p.numShapes)];
285 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp;
286 }
287 }
288 }
289 else
290 {
291 for (s = 0; s < p.numShapes; s++)
292 {
293 for (r = 0; r < p.numShapes; r++)
294 {
295 rtmp = 0;
296 for (q = 0; q < p.numQuad; q++)
297 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
298 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp * D_p[0];
299 }
300 }
301 }
302 }
303 /************************************************************************************/
304 /* process X: */
305 /************************************************************************************/
306 if (NULL != X_p)
307 {
308 add_EM_F = TRUE;
309 if (extendedX)
310 {
311 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuad)]);
312 for (s = 0; s < p.numShapes; s++)
313 {
314 rtmp = 0.;
315 for (q = 0; q < p.numQuad; q++)
316 {
317 rtmp +=
318 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
319 X_q[INDEX2(0, q, DIM)] +
320 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * X_q[INDEX2(1, q, DIM)]);
321 }
322 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
323 }
324 }
325 else
326 {
327 for (s = 0; s < p.numShapes; s++)
328 {
329 rtmp0 = 0.;
330 rtmp1 = 0.;
331 for (q = 0; q < p.numQuad; q++)
332 {
333 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
334 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
335 }
336 EM_F[INDEX2(0, s, p.numEqu)] += rtmp0 * X_p[0] + rtmp1 * X_p[1];
337 }
338 }
339 }
340 /************************************************************************************/
341 /* process Y: */
342 /************************************************************************************/
343 if (NULL != Y_p)
344 {
345 add_EM_F = TRUE;
346 if (extendedY)
347 {
348 Y_q = &(Y_p[INDEX2(0, 0, p.numQuad)]);
349 for (s = 0; s < p.numShapes; s++)
350 {
351 rtmp = 0;
352 for (q = 0; q < p.numQuad; q++)
353 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
354 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
355 }
356 }
357 else
358 {
359 for (s = 0; s < p.numShapes; s++)
360 {
361 rtmp = 0;
362 for (q = 0; q < p.numQuad; q++)
363 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
364 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
365 }
366 }
367 }
368 /*********************************************************************************************************************/
369 /* add the element matrices onto the matrix and right hand side */
370 /*********************************************************************************************************************/
371
372 for (q = 0; q < p.numShapes; q++)
373 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
374 if (add_EM_F)
375 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound);
376 if (add_EM_S)
377 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu,
378 p.numShapes, row_index, p.numComp, EM_S);
379 } /* end color check */
380 } /* end element loop */
381 } /* end color loop */
382
383 delete[] EM_S;
384 delete[] EM_F;
385 delete[] row_index;
386
387 } /* end of pointer check */
388 } /* end parallel region */
389 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26