/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 11392 byte(s)
like sand though the hourglass
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 1D 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 = 1 x 1 */
29 /* B = 1 */
30 /* C = 1 */
31 /* D = scalar */
32 /* X = 1 */
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 /************************************************************************************/
44
45 void Dudley_Assemble_PDE_Single2_1D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements,
46 Paso_SystemMatrix * Mat, escriptDataC * F,
47 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
48 escriptDataC * X, escriptDataC * Y)
49 {
50
51 #define DIM 1
52 index_t color;
53 dim_t e;
54 __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;
55 double *EM_S, *EM_F, *DSDX;
56 index_t *row_index;
57 register dim_t q, s, r;
58 register double rtmp;
59 bool_t add_EM_F, add_EM_S;
60
61 bool_t extendedA = isExpanded(A);
62 bool_t extendedB = isExpanded(B);
63 bool_t extendedC = isExpanded(C);
64 bool_t extendedD = isExpanded(D);
65 bool_t extendedX = isExpanded(X);
66 bool_t extendedY = isExpanded(Y);
67 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
68 double *S = p.row_jac->BasisFunctions->S;
69 dim_t len_EM_S = p.row_numShapesTotal * p.row_numShapesTotal;
70 dim_t len_EM_F = p.row_numShapesTotal;
71
72 #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,rtmp,add_EM_F, add_EM_S)
73 {
74 EM_S = new double[len_EM_S];
75 EM_F = new double[len_EM_F];
76 row_index = new index_t[p.row_numShapesTotal];
77
78 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
79 {
80
81 for (color = elements->minColor; color <= elements->maxColor; color++)
82 {
83 /* open loop over all elements: */
84 #pragma omp for private(e) schedule(static)
85 for (e = 0; e < elements->numElements; e++)
86 {
87 if (elements->Color[e] == color)
88 {
89
90 A_p = getSampleDataRO(A, e);
91 C_p = getSampleDataRO(C, e);
92 B_p = getSampleDataRO(B, e);
93 D_p = getSampleDataRO(D, e);
94 X_p = getSampleDataRO(X, e);
95 Y_p = getSampleDataRO(Y, e);
96
97 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
98 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
99 for (q = 0; q < len_EM_S; ++q)
100 EM_S[q] = 0;
101 for (q = 0; q < len_EM_F; ++q)
102 EM_F[q] = 0;
103 add_EM_F = FALSE;
104 add_EM_S = FALSE;
105 /************************************************************************************/
106 /* process A: */
107 /************************************************************************************/
108 if (NULL != A_p)
109 {
110 add_EM_S = TRUE;
111 if (extendedA)
112 {
113 A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
114 for (s = 0; s < p.row_numShapes; s++)
115 {
116 for (r = 0; r < p.row_numShapes; r++)
117 {
118 rtmp = 0;
119 for (q = 0; q < p.numQuadTotal; q++)
120 {
121 rtmp +=
122 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
123 A_q[INDEX3(0, 0, q, DIM, DIM)] *
124 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
125 }
126 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
127 }
128 }
129 }
130 else
131 {
132 for (s = 0; s < p.row_numShapes; s++)
133 {
134 for (r = 0; r < p.row_numShapes; r++)
135 {
136 rtmp = 0;
137 for (q = 0; q < p.numQuadTotal; q++)
138 rtmp +=
139 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
140 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
141 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
142 rtmp * A_p[INDEX2(0, 0, DIM)];
143 }
144 }
145 }
146 }
147 /************************************************************************************/
148 /* process B: */
149 /************************************************************************************/
150 if (NULL != B_p)
151 {
152 add_EM_S = TRUE;
153 if (extendedB)
154 {
155 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
156 for (s = 0; s < p.row_numShapes; s++)
157 {
158 for (r = 0; r < p.row_numShapes; r++)
159 {
160 rtmp = 0;
161 for (q = 0; q < p.numQuadTotal; q++)
162 {
163 rtmp +=
164 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
165 B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
166 }
167 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
168 }
169 }
170 }
171 else
172 {
173 for (s = 0; s < p.row_numShapes; s++)
174 {
175 for (r = 0; r < p.row_numShapes; r++)
176 {
177 rtmp = 0;
178 for (q = 0; q < p.numQuadTotal; q++)
179 rtmp +=
180 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
181 S[INDEX2(r, q, p.row_numShapes)];
182 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
183 rtmp * B_p[0];
184 }
185 }
186 }
187 }
188 /************************************************************************************/
189 /* process C: */
190 /************************************************************************************/
191 if (NULL != C_p)
192 {
193 add_EM_S = TRUE;
194 if (extendedC)
195 {
196 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
197 for (s = 0; s < p.row_numShapes; s++)
198 {
199 for (r = 0; r < p.row_numShapes; r++)
200 {
201 rtmp = 0;
202 for (q = 0; q < p.numQuadTotal; q++)
203 {
204 rtmp +=
205 vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
206 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
207 }
208 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
209 }
210 }
211 }
212 else
213 {
214 for (s = 0; s < p.row_numShapes; s++)
215 {
216 for (r = 0; r < p.row_numShapes; r++)
217 {
218 rtmp = 0;
219 for (q = 0; q < p.numQuadTotal; q++)
220 rtmp +=
221 vol * S[INDEX2(s, q, p.row_numShapes)] *
222 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
223 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
224 rtmp * C_p[0];
225 }
226 }
227 }
228 }
229 /*********************************************************************************** */
230 /* process D */
231 /************************************************************************************/
232 if (NULL != D_p)
233 {
234 add_EM_S = TRUE;
235 if (extendedD)
236 {
237 D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
238 for (s = 0; s < p.row_numShapes; s++)
239 {
240 for (r = 0; r < p.row_numShapes; r++)
241 {
242 rtmp = 0;
243 for (q = 0; q < p.numQuadTotal; q++)
244 {
245 rtmp +=
246 vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
247 S[INDEX2(r, q, p.row_numShapes)];
248 }
249 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
250 }
251 }
252 }
253 else
254 {
255 for (s = 0; s < p.row_numShapes; s++)
256 {
257 for (r = 0; r < p.row_numShapes; r++)
258 {
259 rtmp = 0;
260 for (q = 0; q < p.numQuadTotal; q++)
261 rtmp +=
262 vol * S[INDEX2(s, q, p.row_numShapes)] *
263 S[INDEX2(r, q, p.row_numShapes)];
264 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
265 rtmp * D_p[0];
266 }
267 }
268 }
269 }
270 /************************************************************************************/
271 /* process X: */
272 /************************************************************************************/
273 if (NULL != X_p)
274 {
275 add_EM_F = TRUE;
276 if (extendedX)
277 {
278 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
279 for (s = 0; s < p.row_numShapes; s++)
280 {
281 rtmp = 0;
282 for (q = 0; q < p.numQuadTotal; q++)
283 rtmp +=
284 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
285 X_q[INDEX2(0, q, DIM)];
286 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
287 }
288 }
289 else
290 {
291 for (s = 0; s < p.row_numShapes; s++)
292 {
293 rtmp = 0;
294 for (q = 0; q < p.numQuadTotal; q++)
295 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
296 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
297 }
298 }
299 }
300 /************************************************************************************/
301 /* process Y: */
302 /************************************************************************************/
303 if (NULL != Y_p)
304 {
305 add_EM_F = TRUE;
306 if (extendedY)
307 {
308 Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
309 for (s = 0; s < p.row_numShapes; s++)
310 {
311 rtmp = 0;
312 for (q = 0; q < p.numQuadTotal; q++)
313 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
314 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
315 }
316 }
317 else
318 {
319 for (s = 0; s < p.row_numShapes; s++)
320 {
321 rtmp = 0;
322 for (q = 0; q < p.numQuadTotal; q++)
323 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
324 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
325 }
326 }
327 }
328 /*********************************************************************************************************************/
329 /* add the element matrices onto the matrix and right hand side */
330 /*********************************************************************************************************************/
331 for (q = 0; q < p.row_numShapesTotal; q++)
332 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
333
334 if (add_EM_F)
335 Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
336 p.row_DOF_UpperBound);
337 if (add_EM_S)
338 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
339 p.row_numShapesTotal, row_index, p.numComp, EM_S);
340
341 } /* end color check */
342 } /* end element loop */
343 } /* end color loop */
344
345 delete[] EM_S; /* these FREEs appear to be inside the if because if any of the allocs */
346 delete[] EM_F; /* failed it means an out of memory (which is not recoverable anyway) */
347 delete[] row_index;
348
349 } /* end of pointer check */
350 } /* end parallel region */
351 }
352
353 /*
354 * $Log$
355 */

  ViewVC Help
Powered by ViewVC 1.1.26