/[escript]/branches/domexper/dudley/src/Assemble_PDE_System2_1D.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_PDE_System2_1D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 12514 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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_{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 */
20
21 /* u has p.numComp components 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 = p.numEqu x 1 x p.numComp x 1 */
27 /* B = 1 x numEqu x p.numComp */
28 /* C = p.numEqu x 1 x p.numComp */
29 /* D = p.numEqu x p.numComp */
30 /* X = p.numEqu x 1 */
31 /* Y = p.numEqu */
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_System2_1D(Dudley_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, k, m;
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_numShapes * p.row_numShapes * p.numEqu * p.numComp;
68 dim_t len_EM_F = p.row_numShapes * p.numEqu;
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,k,m,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_numShapes, 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 B_p = getSampleDataRO(B, e);
90 C_p = getSampleDataRO(C, e);
91 D_p = getSampleDataRO(D, e);
92 X_p = getSampleDataRO(X, e);
93 Y_p = getSampleDataRO(Y, e);
94 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
95 Vol = &(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal, 1)]);
96 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.row_numShapes, DIM, p.numQuadTotal, 1)]);
97 for (q = 0; q < len_EM_S; ++q)
98 EM_S[q] = 0;
99 for (q = 0; q < len_EM_F; ++q)
100 EM_F[q] = 0;
101 add_EM_F = FALSE;
102 add_EM_S = FALSE;
103
104 /**************************************************************/
105 /* process A: */
106 /**************************************************************/
107
108 if (NULL != A_p)
109 {
110 add_EM_S = TRUE;
111 if (extendedA)
112 {
113 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuadTotal)]);
114 for (s = 0; s < p.row_numShapes; s++)
115 {
116 for (r = 0; r < p.row_numShapes; r++)
117 {
118 for (k = 0; k < p.numEqu; k++)
119 {
120 for (m = 0; m < p.numComp; m++)
121 {
122 rtmp = 0.;
123 for (q = 0; q < p.numQuadTotal; q++)
124 {
125 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
126 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] *
127 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
128 }
129 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
130 }
131 }
132 }
133 }
134 }
135 else
136 {
137 for (s = 0; s < p.row_numShapes; s++)
138 {
139 for (r = 0; r < p.row_numShapes; r++)
140 {
141 rtmp = 0;
142 for (q = 0; q < p.numQuadTotal; q++)
143 rtmp +=
144 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
145 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
146 for (k = 0; k < p.numEqu; k++)
147 {
148 for (m = 0; m < p.numComp; m++)
149 {
150 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
151 rtmp * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)];
152 }
153 }
154 }
155 }
156 }
157 }
158 /**************************************************************/
159 /* process B: */
160 /**************************************************************/
161
162 if (NULL != B_p)
163 {
164 add_EM_S = TRUE;
165 if (extendedB)
166 {
167 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuadTotal)]);
168 for (s = 0; s < p.row_numShapes; s++)
169 {
170 for (r = 0; r < p.row_numShapes; r++)
171 {
172 for (k = 0; k < p.numEqu; k++)
173 {
174 for (m = 0; m < p.numComp; m++)
175 {
176 rtmp = 0.;
177 for (q = 0; q < p.numQuadTotal; q++)
178 {
179 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
180 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] *
181 S[INDEX2(r, q, p.row_numShapes)];
182 }
183 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
184 }
185 }
186 }
187 }
188 }
189 else
190 {
191 for (s = 0; s < p.row_numShapes; s++)
192 {
193 for (r = 0; r < p.row_numShapes; r++)
194 {
195 rtmp = 0;
196 for (q = 0; q < p.numQuadTotal; q++)
197 rtmp +=
198 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
199 S[INDEX2(r, q, p.row_numShapes)];
200 for (k = 0; k < p.numEqu; k++)
201 {
202 for (m = 0; m < p.numComp; m++)
203 {
204 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
205 rtmp * B_p[INDEX3(k, 0, m, p.numEqu, DIM)];
206 }
207 }
208 }
209 }
210 }
211 }
212 /**************************************************************/
213 /* process C: */
214 /**************************************************************/
215
216 if (NULL != C_p)
217 {
218 add_EM_S = TRUE;
219 if (extendedC)
220 {
221 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuadTotal)]);
222 for (s = 0; s < p.row_numShapes; s++)
223 {
224 for (r = 0; r < p.row_numShapes; r++)
225 {
226 for (k = 0; k < p.numEqu; k++)
227 {
228 for (m = 0; m < p.numComp; m++)
229 {
230 rtmp = 0;
231 for (q = 0; q < p.numQuadTotal; q++)
232 {
233 rtmp +=
234 vol * S[INDEX2(s, q, p.row_numShapes)] *
235 C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
236 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
237 }
238 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
239 }
240 }
241 }
242 }
243 }
244 else
245 {
246 for (s = 0; s < p.row_numShapes; s++)
247 {
248 for (r = 0; r < p.row_numShapes; r++)
249 {
250 rtmp = 0;
251 for (q = 0; q < p.numQuadTotal; q++)
252 rtmp +=
253 vol * S[INDEX2(s, q, p.row_numShapes)] *
254 DSDX[INDEX3(r, 0, q, p.row_numShapes, DIM)];
255 for (k = 0; k < p.numEqu; k++)
256 {
257 for (m = 0; m < p.numComp; m++)
258 {
259 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
260 rtmp * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)];
261 }
262 }
263 }
264 }
265 }
266 }
267 /************************************************************* */
268 /* process D */
269 /**************************************************************/
270
271 if (NULL != D_p)
272 {
273 add_EM_S = TRUE;
274 if (extendedD)
275 {
276 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuadTotal)]);
277 for (s = 0; s < p.row_numShapes; s++)
278 {
279 for (r = 0; r < p.row_numShapes; r++)
280 {
281 for (k = 0; k < p.numEqu; k++)
282 {
283 for (m = 0; m < p.numComp; m++)
284 {
285 rtmp = 0;
286 for (q = 0; q < p.numQuadTotal; q++)
287 {
288 rtmp +=
289 vol * S[INDEX2(s, q, p.row_numShapes)] *
290 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
291 S[INDEX2(r, q, p.row_numShapes)];
292 }
293 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] += rtmp;
294
295 }
296 }
297 }
298 }
299 }
300 else
301 {
302 for (s = 0; s < p.row_numShapes; s++)
303 {
304 for (r = 0; r < p.row_numShapes; r++)
305 {
306 rtmp = 0;
307 for (q = 0; q < p.numQuadTotal; q++)
308 rtmp +=
309 vol * S[INDEX2(s, q, p.row_numShapes)] *
310 S[INDEX2(r, q, p.row_numShapes)];
311 for (k = 0; k < p.numEqu; k++)
312 {
313 for (m = 0; m < p.numComp; m++)
314 {
315 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.row_numShapes)] +=
316 rtmp * D_p[INDEX2(k, m, p.numEqu)];
317 }
318 }
319 }
320 }
321 }
322 }
323 /**************************************************************/
324 /* process X: */
325 /**************************************************************/
326
327 if (NULL != X_p)
328 {
329 add_EM_F = TRUE;
330 if (extendedX)
331 {
332 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuadTotal)]);
333 for (s = 0; s < p.row_numShapes; s++)
334 {
335 for (k = 0; k < p.numEqu; k++)
336 {
337 rtmp = 0;
338 for (q = 0; q < p.numQuadTotal; q++)
339 rtmp +=
340 vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)] *
341 X_q[INDEX3(k, 0, q, p.numEqu, DIM)];
342 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
343 }
344 }
345 }
346 else
347 {
348 for (s = 0; s < p.row_numShapes; s++)
349 {
350 rtmp = 0;
351 for (q = 0; q < p.numQuadTotal; q++)
352 rtmp += vol * DSDX[INDEX3(s, 0, q, p.row_numShapes, DIM)];
353 for (k = 0; k < p.numEqu; k++)
354 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * X_p[INDEX2(k, 0, p.numEqu)];
355 }
356 }
357 }
358 /**************************************************************/
359 /* process Y: */
360 /**************************************************************/
361
362 if (NULL != Y_p)
363 {
364 add_EM_F = TRUE;
365 if (extendedY)
366 {
367 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuadTotal)]);
368 for (s = 0; s < p.row_numShapes; s++)
369 {
370 for (k = 0; k < p.numEqu; k++)
371 {
372 rtmp = 0;
373 for (q = 0; q < p.numQuadTotal; q++)
374 rtmp +=
375 vol * S[INDEX2(s, q, p.row_numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
376 EM_F[INDEX2(k, s, p.numEqu)] += rtmp;
377 }
378 }
379 }
380 else
381 {
382 for (s = 0; s < p.row_numShapes; s++)
383 {
384 rtmp = 0;
385 for (q = 0; q < p.numQuadTotal; q++)
386 rtmp += vol * S[INDEX2(s, q, p.row_numShapes)];
387 for (k = 0; k < p.numEqu; k++)
388 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k];
389 }
390 }
391 }
392 /***********************************************************************************************/
393 /* add the element matrices onto the matrix and right hand side */
394 /***********************************************************************************************/
395
396 for (q = 0; q < p.row_numShapes; q++)
397 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
398
399 if (add_EM_F)
400 Dudley_Util_AddScatter(p.row_numShapes, row_index, p.numEqu, EM_F, F_p,
401 p.row_DOF_UpperBound);
402 if (add_EM_S)
403 Dudley_Assemble_addToSystemMatrix(Mat, p.row_numShapes, row_index, p.numEqu,
404 p.row_numShapes, row_index, p.numComp, EM_S);
405
406 } /* end color check */
407 } /* end element loop */
408 } /* end color loop */
409
410 THREAD_MEMFREE(EM_S);
411 THREAD_MEMFREE(EM_F);
412 THREAD_MEMFREE(row_index);
413
414 } /* end of pointer check */
415 } /* end parallel region */
416 }
417
418 /*
419 * $Log$
420 */

  ViewVC Help
Powered by ViewVC 1.1.26