/[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 3187 - (show annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 7 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_Single2_1D.c
File MIME type: text/plain
File size: 10990 byte(s)
Enforcing indent -linux -i4 -bl -bli0 -l120

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 /**************************************************************/
15
16 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
17 /* the shape functions for test and solution must be identical */
18
19 /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
20
21 /* in a 1D domain. The shape functions for test and solution must be identical */
22 /* and row_NS == row_NN */
23
24 /* Shape of the coefficients: */
25
26 /* A = 1 x 1 */
27 /* B = 1 */
28 /* C = 1 */
29 /* D = scalar */
30 /* X = 1 */
31 /* Y = scalar */
32
33 /**************************************************************/
34
35 #include "Assemble.h"
36 #include "Util.h"
37 #ifdef _OPENMP
38 #include <omp.h>
39 #endif
40
41 /**************************************************************/
42
43 void Dudley_Assemble_PDE_Single2_1D(Assemble_Parameters p, Dudley_ElementFile * elements,
44 Paso_SystemMatrix * Mat, escriptDataC * F,
45 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
46 escriptDataC * X, escriptDataC * Y)
47 {
48
49 #define DIM 1
50 index_t color;
51 dim_t e;
52 __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;
53 double *EM_S, *EM_F, *DSDX;
54 index_t *row_index;
55 register dim_t q, s, r;
56 register double rtmp;
57 bool_t add_EM_F, add_EM_S;
58
59 bool_t extendedA = isExpanded(A);
60 bool_t extendedB = isExpanded(B);
61 bool_t extendedC = isExpanded(C);
62 bool_t extendedD = isExpanded(D);
63 bool_t extendedX = isExpanded(X);
64 bool_t extendedY = isExpanded(Y);
65 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */
66 double *S = p.row_jac->BasisFunctions->S;
67 dim_t len_EM_S = p.row_numShapesTotal * p.col_numShapesTotal;
68 dim_t len_EM_F = p.row_numShapesTotal;
69
70 #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)
71 {
72 EM_S = THREAD_MEMALLOC(len_EM_S, double);
73 EM_F = THREAD_MEMALLOC(len_EM_F, double);
74 row_index = THREAD_MEMALLOC(p.row_numShapesTotal, index_t);
75
76 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index))
77 {
78
79 for (color = elements->minColor; color <= elements->maxColor; color++)
80 {
81 /* open loop over all elements: */
82 #pragma omp for private(e) schedule(static)
83 for (e = 0; e < elements->numElements; e++)
84 {
85 if (elements->Color[e] == color)
86 {
87
88 A_p = getSampleDataRO(A, e);
89 C_p = getSampleDataRO(C, e);
90 B_p = getSampleDataRO(B, e);
91 D_p = getSampleDataRO(D, e);
92 X_p = getSampleDataRO(X, e);
93 Y_p = getSampleDataRO(Y, e);
94
95 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
96 // Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
97 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapesTotal, DIM, p.numQuadTotal, 1)]);
98 for (q = 0; q < len_EM_S; ++q)
99 EM_S[q] = 0;
100 for (q = 0; q < len_EM_F; ++q)
101 EM_F[q] = 0;
102 add_EM_F = FALSE;
103 add_EM_S = FALSE;
104 /**************************************************************/
105 /* process A: */
106 /**************************************************************/
107 if (NULL != A_p)
108 {
109 add_EM_S = TRUE;
110 if (extendedA)
111 {
112 A_q = &(A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuadTotal)]);
113 for (s = 0; s < p.row_numShapes; s++)
114 {
115 for (r = 0; r < p.col_numShapes; r++)
116 {
117 rtmp = 0;
118 for (q = 0; q < p.numQuadTotal; q++)
119 {
120 rtmp +=
121 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
122 A_q[INDEX3(0, 0, q, DIM, DIM)] *
123 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
124 }
125 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
126 }
127 }
128 } else
129 {
130 for (s = 0; s < p.row_numShapes; s++)
131 {
132 for (r = 0; r < p.col_numShapes; r++)
133 {
134 rtmp = 0;
135 for (q = 0; q < p.numQuadTotal; q++)
136 rtmp +=
137 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
138 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
139 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
140 rtmp * A_p[INDEX2(0, 0, DIM)];
141 }
142 }
143 }
144 }
145 /**************************************************************/
146 /* process B: */
147 /**************************************************************/
148 if (NULL != B_p)
149 {
150 add_EM_S = TRUE;
151 if (extendedB)
152 {
153 B_q = &(B_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
154 for (s = 0; s < p.row_numShapes; s++)
155 {
156 for (r = 0; r < p.col_numShapes; r++)
157 {
158 rtmp = 0;
159 for (q = 0; q < p.numQuadTotal; q++)
160 {
161 rtmp +=
162 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
163 B_q[INDEX2(0, q, DIM)] * S[INDEX2(r, q, p.row_numShapes)];
164 }
165 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
166 }
167 }
168 } else
169 {
170 for (s = 0; s < p.row_numShapes; s++)
171 {
172 for (r = 0; r < p.col_numShapes; r++)
173 {
174 rtmp = 0;
175 for (q = 0; q < p.numQuadTotal; q++)
176 rtmp +=
177 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
178 S[INDEX2(r, q, p.row_numShapes)];
179 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
180 rtmp * B_p[0];
181 }
182 }
183 }
184 }
185 /**************************************************************/
186 /* process C: */
187 /**************************************************************/
188 if (NULL != C_p)
189 {
190 add_EM_S = TRUE;
191 if (extendedC)
192 {
193 C_q = &(C_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
194 for (s = 0; s < p.row_numShapes; s++)
195 {
196 for (r = 0; r < p.col_numShapes; r++)
197 {
198 rtmp = 0;
199 for (q = 0; q < p.numQuadTotal; q++)
200 {
201 rtmp +=
202 vol * S[INDEX2(s, q, p.row_numShapes)] * C_q[INDEX2(0, q, DIM)] *
203 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
204 }
205 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
206 }
207 }
208 } else
209 {
210 for (s = 0; s < p.row_numShapes; s++)
211 {
212 for (r = 0; r < p.col_numShapes; r++)
213 {
214 rtmp = 0;
215 for (q = 0; q < p.numQuadTotal; q++)
216 rtmp +=
217 vol * S[INDEX2(s, q, p.row_numShapes)] *
218 DSDX[INDEX3(r, 0, q, p.row_numShapesTotal, DIM)];
219 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
220 rtmp * C_p[0];
221 }
222 }
223 }
224 }
225 /************************************************************* */
226 /* process D */
227 /**************************************************************/
228 if (NULL != D_p)
229 {
230 add_EM_S = TRUE;
231 if (extendedD)
232 {
233 D_q = &(D_p[INDEX2(0, 0, p.numQuadTotal)]);
234 for (s = 0; s < p.row_numShapes; s++)
235 {
236 for (r = 0; r < p.col_numShapes; r++)
237 {
238 rtmp = 0;
239 for (q = 0; q < p.numQuadTotal; q++)
240 {
241 rtmp +=
242 vol * S[INDEX2(s, q, p.row_numShapes)] * D_q[q] *
243 S[INDEX2(r, q, p.row_numShapes)];
244 }
245 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] += rtmp;
246 }
247 }
248 } else
249 {
250 for (s = 0; s < p.row_numShapes; s++)
251 {
252 for (r = 0; r < p.col_numShapes; r++)
253 {
254 rtmp = 0;
255 for (q = 0; q < p.numQuadTotal; q++)
256 rtmp +=
257 vol * S[INDEX2(s, q, p.row_numShapes)] *
258 S[INDEX2(r, q, p.row_numShapes)];
259 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.row_numShapesTotal)] +=
260 rtmp * D_p[0];
261 }
262 }
263 }
264 }
265 /**************************************************************/
266 /* process X: */
267 /**************************************************************/
268 if (NULL != X_p)
269 {
270 add_EM_F = TRUE;
271 if (extendedX)
272 {
273 X_q = &(X_p[INDEX3(0, 0, 0, DIM, p.numQuadTotal)]);
274 for (s = 0; s < p.row_numShapes; s++)
275 {
276 rtmp = 0;
277 for (q = 0; q < p.numQuadTotal; q++)
278 rtmp +=
279 vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)] *
280 X_q[INDEX2(0, q, DIM)];
281 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
282 }
283 } else
284 {
285 for (s = 0; s < p.row_numShapes; s++)
286 {
287 rtmp = 0;
288 for (q = 0; q < p.numQuadTotal; q++)
289 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapesTotal, DIM)];
290 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * X_p[0];
291 }
292 }
293 }
294 /**************************************************************/
295 /* process Y: */
296 /**************************************************************/
297 if (NULL != Y_p)
298 {
299 add_EM_F = TRUE;
300 if (extendedY)
301 {
302 Y_q = &(Y_p[INDEX2(0, 0, p.numQuadTotal)]);
303 for (s = 0; s < p.row_numShapes; s++)
304 {
305 rtmp = 0;
306 for (q = 0; q < p.numQuadTotal; q++)
307 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[q];
308 EM_F[INDEX2(0, s, p.numEqu)] += rtmp;
309 }
310 } else
311 {
312 for (s = 0; s < p.row_numShapes; s++)
313 {
314 rtmp = 0;
315 for (q = 0; q < p.numQuadTotal; q++)
316 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
317 EM_F[INDEX2(0, s, p.numEqu)] += rtmp * Y_p[0];
318 }
319 }
320 }
321 /***********************************************************************************************/
322 /* add the element matrices onto the matrix and right hand side */
323 /***********************************************************************************************/
324 for (q = 0; q < p.row_numShapesTotal; q++)
325 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
326
327 if (add_EM_F)
328 Dudley_Util_AddScatter(p.row_numShapesTotal, row_index, p.numEqu, EM_F, F_p,
329 p.row_DOF_UpperBound);
330 if (add_EM_S)
331 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapesTotal, row_index, p.numEqu,
332 p.col_numShapesTotal, row_index, p.numComp, EM_S);
333
334 } /* end color check */
335 } /* end element loop */
336 } /* end color loop */
337
338 THREAD_MEMFREE(EM_S); /* these FREEs appear to be inside the if because if any of the allocs */
339 THREAD_MEMFREE(EM_F); /* failed it means an out of memory (which is not recoverable anyway) */
340 THREAD_MEMFREE(row_index);
341
342 } /* end of pointer check */
343 } /* end parallel region */
344 }
345
346 /*
347 * $Log$
348 */

  ViewVC Help
Powered by ViewVC 1.1.26