/[escript]/trunk/finley/src/Assemble_PDE_System2_3D.c
ViewVC logotype

Contents of /trunk/finley/src/Assemble_PDE_System2_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 20813 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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
17 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
18 /* the shape functions for test and solution must be identical */
19
20
21 /* -(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 */
22
23 /* u has p.numComp components in a 3D 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 = p.numEqu x 3 x p.numComp x 3 */
29 /* B = 3 x p.numEqu x p.numComp */
30 /* C = p.numEqu x 3 x p.numComp */
31 /* D = p.numEqu x p.numComp */
32 /* X = p.numEqu x 3 */
33 /* Y = p.numEqu */
34
35
36 /**************************************************************/
37
38
39 #include "Assemble.h"
40 #include "Util.h"
41 #ifdef _OPENMP
42 #include <omp.h>
43 #endif
44
45
46 /**************************************************************/
47
48 void Finley_Assemble_PDE_System2_3D(Assemble_Parameters p, Finley_ElementFile* elements,
49 Paso_SystemMatrix* Mat, escriptDataC* F,
50 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
51
52 #define DIM 3
53 index_t color;
54 dim_t e;
55 double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;
56 index_t *row_index;
57 register dim_t q, s,r,k,m;
58 register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
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=getSampleData(F,0);
68 double *S=p.row_jac->ReferenceElement->S;
69 dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;
70 dim_t len_EM_F=p.row_NN*p.numEqu;
71
72
73 #pragma omp parallel private(color,EM_S, EM_F, Vol, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p,row_index,q, s,r,k,m,rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22,add_EM_F, add_EM_S)
74 {
75 EM_S=THREAD_MEMALLOC(len_EM_S,double);
76 EM_F=THREAD_MEMALLOC(len_EM_F,double);
77 row_index=THREAD_MEMALLOC(p.row_NN,index_t);
78
79 if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
80
81 for (color=elements->minColor;color<=elements->maxColor;color++) {
82 /* open loop over all elements: */
83 #pragma omp for private(e) schedule(static)
84 for(e=0;e<elements->numElements;e++){
85 if (elements->Color[e]==color) {
86 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
87 DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
88 for (q=0;q<len_EM_S;++q) EM_S[q]=0;
89 for (q=0;q<len_EM_F;++q) EM_F[q]=0;
90 add_EM_F=FALSE;
91 add_EM_S=FALSE;
92 /**************************************************************/
93 /* process A: */
94 /**************************************************************/
95 A_p=getSampleData(A,e);
96 if (NULL!=A_p) {
97 add_EM_S=TRUE;
98 if (extendedA) {
99 for (s=0;s<p.row_NS;s++) {
100 for (r=0;r<p.col_NS;r++) {
101 for (k=0;k<p.numEqu;k++) {
102 for (m=0;m<p.numComp;m++) {
103 rtmp=0;
104 for (q=0;q<p.numQuad;q++) {
105 rtmp+=Vol[q]*(
106 DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
107 +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
108 +DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX5(k,0,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]
109 +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
110 +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
111 +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX5(k,1,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]
112 +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
113 +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
114 +DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*A_p[INDEX5(k,2,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);
115
116 }
117 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
118 }
119 }
120 }
121 }
122 } else {
123 for (s=0;s<p.row_NS;s++) {
124 for (r=0;r<p.col_NS;r++) {
125 rtmp00=0;
126 rtmp01=0;
127 rtmp02=0;
128 rtmp10=0;
129 rtmp11=0;
130 rtmp12=0;
131 rtmp20=0;
132 rtmp21=0;
133 rtmp22=0;
134 for (q=0;q<p.numQuad;q++) {
135 rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
136 rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
137 rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
138 rtmp02+=rtmp0*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
139
140 rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
141 rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
142 rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
143 rtmp12+=rtmp1*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
144
145 rtmp2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
146 rtmp20+=rtmp2*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
147 rtmp21+=rtmp2*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
148 rtmp22+=rtmp2*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
149 }
150 for (k=0;k<p.numEqu;k++) {
151 for (m=0;m<p.numComp;m++) {
152 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
153 +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
154 +rtmp02*A_p[INDEX4(k,0,m,2,p.numEqu,DIM,p.numComp)]
155 +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
156 +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)]
157 +rtmp12*A_p[INDEX4(k,1,m,2,p.numEqu,DIM,p.numComp)]
158 +rtmp20*A_p[INDEX4(k,2,m,0,p.numEqu,DIM,p.numComp)]
159 +rtmp21*A_p[INDEX4(k,2,m,1,p.numEqu,DIM,p.numComp)]
160 +rtmp22*A_p[INDEX4(k,2,m,2,p.numEqu,DIM,p.numComp)];
161 }
162 }
163 }
164 }
165 }
166 }
167 /**************************************************************/
168 /* process B: */
169 /**************************************************************/
170 B_p=getSampleData(B,e);
171 if (NULL!=B_p) {
172 add_EM_S=TRUE;
173 if (extendedB) {
174 for (s=0;s<p.row_NS;s++) {
175 for (r=0;r<p.col_NS;r++) {
176 for (k=0;k<p.numEqu;k++) {
177 for (m=0;m<p.numComp;m++) {
178 rtmp=0;
179 for (q=0;q<p.numQuad;q++) {
180 rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*(
181 DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
182 + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]
183 + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*B_p[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)] );
184 }
185 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
186 }
187 }
188 }
189 }
190 } else {
191 for (s=0;s<p.row_NS;s++) {
192 for (r=0;r<p.col_NS;r++) {
193 rtmp0=0;
194 rtmp1=0;
195 rtmp2=0;
196 for (q=0;q<p.numQuad;q++) {
197 rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];
198 rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
199 rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
200 rtmp2+=rtmp*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
201 }
202 for (k=0;k<p.numEqu;k++) {
203 for (m=0;m<p.numComp;m++) {
204 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
205 +rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]
206 +rtmp2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];
207 }
208 }
209 }
210 }
211 }
212 }
213 /**************************************************************/
214 /* process C: */
215 /**************************************************************/
216 C_p=getSampleData(C,e);
217 if (NULL!=C_p) {
218 add_EM_S=TRUE;
219 if (extendedC) {
220 for (s=0;s<p.row_NS;s++) {
221 for (r=0;r<p.col_NS;r++) {
222 for (k=0;k<p.numEqu;k++) {
223 for (m=0;m<p.numComp;m++) {
224 rtmp=0;
225 for (q=0;q<p.numQuad;q++) {
226 rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*(
227 C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
228 +C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
229 +C_p[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_NN,DIM)]);
230 }
231 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
232 }
233 }
234 }
235 }
236 } else {
237 for (s=0;s<p.row_NS;s++) {
238 for (r=0;r<p.col_NS;r++) {
239 rtmp0=0;
240 rtmp1=0;
241 rtmp2=0;
242 for (q=0;q<p.numQuad;q++) {
243 rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];
244 rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
245 rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
246 rtmp2+=rtmp*DSDX[INDEX3(r,2,q,p.row_NN,DIM)];
247 }
248 for (k=0;k<p.numEqu;k++) {
249 for (m=0;m<p.numComp;m++) {
250 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]
251 +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]
252 +rtmp2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];
253 }
254 }
255 }
256 }
257 }
258 }
259 /************************************************************* */
260 /* process D */
261 /**************************************************************/
262 D_p=getSampleData(D,e);
263 if (NULL!=D_p) {
264 add_EM_S=TRUE;
265 if (extendedD) {
266 for (s=0;s<p.row_NS;s++) {
267 for (r=0;r<p.col_NS;r++) {
268 for (k=0;k<p.numEqu;k++) {
269 for (m=0;m<p.numComp;m++) {
270 rtmp=0;
271 for (q=0;q<p.numQuad;q++) {
272 rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_NS)];
273 }
274 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
275 }
276 }
277 }
278 }
279 } else {
280 for (s=0;s<p.row_NS;s++) {
281 for (r=0;r<p.col_NS;r++) {
282 rtmp=0;
283 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
284 for (k=0;k<p.numEqu;k++) {
285 for (m=0;m<p.numComp;m++) {
286 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
287 }
288 }
289 }
290 }
291 }
292 }
293 /**************************************************************/
294 /* process X: */
295 /**************************************************************/
296 X_p=getSampleData(X,e);
297 if (NULL!=X_p) {
298 add_EM_F=TRUE;
299 if (extendedX) {
300 for (s=0;s<p.row_NS;s++) {
301 for (k=0;k<p.numEqu;k++) {
302 rtmp=0;
303 for (q=0;q<p.numQuad;q++) {
304 rtmp+=Vol[q]* ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]
305 + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]
306 + DSDX[INDEX3(s,2,q,p.row_NN,DIM)]*X_p[INDEX3(k,2,q,p.numEqu,DIM)]);
307 }
308 EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
309 }
310 }
311 } else {
312 for (s=0;s<p.row_NS;s++) {
313 rtmp0=0;
314 rtmp1=0;
315 rtmp2=0;
316 for (q=0;q<p.numQuad;q++) {
317 rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
318 rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
319 rtmp2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_NN,DIM)];
320 }
321 for (k=0;k<p.numEqu;k++) {
322 EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]
323 +rtmp1*X_p[INDEX2(k,1,p.numEqu)]
324 +rtmp2*X_p[INDEX2(k,2,p.numEqu)];
325 }
326 }
327 }
328 }
329 /**************************************************************/
330 /* process Y: */
331 /**************************************************************/
332 Y_p=getSampleData(Y,e);
333 if (NULL!=Y_p) {
334 add_EM_F=TRUE;
335 if (extendedY) {
336 for (s=0;s<p.row_NS;s++) {
337 for (k=0;k<p.numEqu;k++) {
338 rtmp=0.;
339 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
340 EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
341 }
342 }
343 } else {
344 for (s=0;s<p.row_NS;s++) {
345 rtmp=0;
346 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
347 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
348 }
349 }
350 }
351 /***********************************************************************************************/
352 /* add the element matrices onto the matrix and right hand side */
353 /***********************************************************************************************/
354 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
355 if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
356 if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
357
358 } /* end color check */
359 } /* end element loop */
360 } /* end color loop */
361
362 THREAD_MEMFREE(EM_S);
363 THREAD_MEMFREE(EM_F);
364 THREAD_MEMFREE(row_index);
365
366 } /* end of pointer check */
367 } /* end parallel region */
368 }
369 /*
370 * $Log$
371 */

  ViewVC Help
Powered by ViewVC 1.1.26