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

Annotation of /trunk/finley/src/Assemble_PDE_System2_2D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 10 months ago) by trankine
File MIME type: text/plain
File size: 17143 byte(s)
And get the *(&(*&(* name right
1 gross 798
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 gross 798 /**************************************************************/
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    
22     /* -(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 */
23    
24     /* u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical */
25     /* and row_NS == row_NN */
26    
27     /* Shape of the coefficients: */
28    
29     /* A = p.numEqu x 2 x p.numComp x 2 */
30     /* B = 2 x p.numEqu x p.numComp */
31     /* C = p.numEqu x 2 x p.numComp */
32     /* D = p.numEqu x p.numComp */
33     /* X = p.numEqu x 2 */
34     /* Y = p.numEqu */
35    
36    
37     /**************************************************************/
38    
39    
40     #include "Assemble.h"
41     #include "Util.h"
42 gross 853 #ifdef _OPENMP
43     #include <omp.h>
44     #endif
45 gross 798
46 gross 853
47 gross 798 /**************************************************************/
48    
49     void Finley_Assemble_PDE_System2_2D(Assemble_Parameters p, Finley_ElementFile* elements,
50     Paso_SystemMatrix* Mat, escriptDataC* F,
51     escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
52    
53     #define DIM 2
54     index_t color;
55     dim_t e;
56 gross 853 double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;
57     index_t *row_index;
58     register dim_t q, s,r,k,m;
59     register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11;
60     bool_t add_EM_F, add_EM_S;
61    
62 gross 798 bool_t extendedA=isExpanded(A);
63     bool_t extendedB=isExpanded(B);
64     bool_t extendedC=isExpanded(C);
65     bool_t extendedD=isExpanded(D);
66     bool_t extendedX=isExpanded(X);
67     bool_t extendedY=isExpanded(Y);
68     double *F_p=getSampleData(F,0);
69     double *S=p.row_jac->ReferenceElement->S;
70     dim_t len_EM_S=p.row_NN*p.col_NN*p.numEqu*p.numComp;
71     dim_t len_EM_F=p.row_NN*p.numEqu;
72    
73    
74 gross 853 #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, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S)
75 gross 798 {
76 gross 853 EM_S=THREAD_MEMALLOC(len_EM_S,double);
77     EM_F=THREAD_MEMALLOC(len_EM_F,double);
78     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
79    
80     if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
81 gross 798
82 gross 853 for (color=elements->minColor;color<=elements->maxColor;color++) {
83     /* open loop over all elements: */
84     #pragma omp for private(e) schedule(static)
85     for(e=0;e<elements->numElements;e++){
86     if (elements->Color[e]==color) {
87     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
88     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
89     for (q=0;q<len_EM_S;++q) EM_S[q]=0;
90     for (q=0;q<len_EM_F;++q) EM_F[q]=0;
91     add_EM_F=FALSE;
92     add_EM_S=FALSE;
93     /**************************************************************/
94     /* process A: */
95     /**************************************************************/
96     A_p=getSampleData(A,e);
97     if (NULL!=A_p) {
98     add_EM_S=TRUE;
99     if (extendedA) {
100     for (s=0;s<p.row_NS;s++) {
101     for (r=0;r<p.col_NS;r++) {
102     for (k=0;k<p.numEqu;k++) {
103     for (m=0;m<p.numComp;m++) {
104     rtmp=0;
105     for (q=0;q<p.numQuad;q++) {
106     rtmp+=Vol[q]* (
107     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)]
108     +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)]
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     }
112     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
113     }
114     }
115     }
116     }
117     } else {
118     for (s=0;s<p.row_NS;s++) {
119     for (r=0;r<p.col_NS;r++) {
120     rtmp00=0;
121     rtmp01=0;
122     rtmp10=0;
123     rtmp11=0;
124 gross 798 for (q=0;q<p.numQuad;q++) {
125 gross 853 rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
126     rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
127     rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
128     rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
129     rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
130     rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
131 gross 798 }
132 gross 853 for (k=0;k<p.numEqu;k++) {
133     for (m=0;m<p.numComp;m++) {
134     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=
135     rtmp00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
136     +rtmp01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
137     +rtmp10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
138     +rtmp11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];
139     }
140     }
141 gross 798 }
142 gross 853 }
143     }
144     }
145     /**************************************************************/
146     /* process B: */
147     /**************************************************************/
148     B_p=getSampleData(B,e);
149     if (NULL!=B_p) {
150     add_EM_S=TRUE;
151     if (extendedB) {
152     for (s=0;s<p.row_NS;s++) {
153     for (r=0;r<p.col_NS;r++) {
154     for (k=0;k<p.numEqu;k++) {
155 gross 798 for (m=0;m<p.numComp;m++) {
156 gross 853 rtmp=0;
157     for (q=0;q<p.numQuad;q++) {
158     rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*
159     ( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
160     + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);
161     }
162     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
163 gross 798 }
164 gross 853 }
165 gross 798 }
166 gross 853 }
167     } else {
168     for (s=0;s<p.row_NS;s++) {
169     for (r=0;r<p.col_NS;r++) {
170     rtmp0=0;
171     rtmp1=0;
172 gross 798 for (q=0;q<p.numQuad;q++) {
173 gross 853 rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];
174     rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
175     rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
176 gross 798 }
177 gross 853 for (k=0;k<p.numEqu;k++) {
178     for (m=0;m<p.numComp;m++) {
179     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
180     + rtmp1*B_p[INDEX3(k,1,m,p.numEqu,DIM)];
181     }
182     }
183     }
184     }
185 gross 798 }
186 gross 853 }
187     /**************************************************************/
188     /* process C: */
189     /**************************************************************/
190     C_p=getSampleData(C,e);
191     if (NULL!=C_p) {
192     add_EM_S=TRUE;
193     if (extendedC) {
194     for (s=0;s<p.row_NS;s++) {
195     for (r=0;r<p.col_NS;r++) {
196     for (k=0;k<p.numEqu;k++) {
197 gross 798 for (m=0;m<p.numComp;m++) {
198 gross 853 rtmp=0;
199     for (q=0;q<p.numQuad;q++) {
200     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*
201     ( C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
202     + C_p[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);
203     }
204     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
205 gross 798 }
206 gross 853 }
207     }
208     }
209     } else {
210     for (s=0;s<p.row_NS;s++) {
211     for (r=0;r<p.col_NS;r++) {
212     rtmp0=0;
213     rtmp1=0;
214     for (q=0;q<p.numQuad;q++) {
215     rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];
216     rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
217     rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
218     }
219     for (k=0;k<p.numEqu;k++) {
220     for (m=0;m<p.numComp;m++) {
221     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)]
222     +rtmp1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];
223     }
224     }
225     }
226     }
227     }
228 gross 798 }
229 gross 853 /************************************************************* */
230     /* process D */
231     /**************************************************************/
232     D_p=getSampleData(D,e);
233     if (NULL!=D_p) {
234     add_EM_S=TRUE;
235     if (extendedD) {
236     for (s=0;s<p.row_NS;s++) {
237     for (r=0;r<p.col_NS;r++) {
238     for (k=0;k<p.numEqu;k++) {
239     for (m=0;m<p.numComp;m++) {
240     rtmp=0;
241     for (q=0;q<p.numQuad;q++) {
242     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)];
243     }
244     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
245     }
246     }
247     }
248     }
249     } else {
250     for (s=0;s<p.row_NS;s++) {
251     for (r=0;r<p.col_NS;r++) {
252 gross 798 rtmp=0;
253 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
254     for (k=0;k<p.numEqu;k++) {
255     for (m=0;m<p.numComp;m++) {
256     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
257     }
258 gross 798 }
259     }
260 gross 853 }
261     }
262     }
263     /**************************************************************/
264     /* process X: */
265     /**************************************************************/
266     X_p=getSampleData(X,e);
267     if (NULL!=X_p) {
268     add_EM_F=TRUE;
269     if (extendedX) {
270     for (s=0;s<p.row_NS;s++) {
271 gross 798 for (k=0;k<p.numEqu;k++) {
272 gross 853 rtmp=0;
273     for (q=0;q<p.numQuad;q++) {
274     rtmp+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)]
275     +DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX3(k,1,q,p.numEqu,DIM)]);
276     }
277     EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
278 gross 798 }
279     }
280 gross 853 } else {
281     for (s=0;s<p.row_NS;s++) {
282     rtmp0=0;
283     rtmp1=0;
284 gross 798 for (q=0;q<p.numQuad;q++) {
285 gross 853 rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
286     rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
287 gross 798 }
288 gross 853 for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp0*X_p[INDEX2(k,0,p.numEqu)]+rtmp1*X_p[INDEX2(k,1,p.numEqu)];
289 gross 798 }
290 gross 853 }
291     }
292     /**************************************************************/
293     /* process Y: */
294     /**************************************************************/
295     Y_p=getSampleData(Y,e);
296     if (NULL!=Y_p) {
297     add_EM_F=TRUE;
298     if (extendedY) {
299     for (s=0;s<p.row_NS;s++) {
300     for (k=0;k<p.numEqu;k++) {
301     rtmp=0;
302     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
303     EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
304     }
305     }
306     } else {
307     for (s=0;s<p.row_NS;s++) {
308 gross 798 rtmp=0;
309 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
310     for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
311 gross 798 }
312     }
313 gross 853 }
314     /***********************************************************************************************/
315     /* add the element matrices onto the matrix and right hand side */
316     /***********************************************************************************************/
317     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
318     if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
319     if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
320    
321     } /* end color check */
322     } /* end element loop */
323     } /* end color loop */
324    
325     THREAD_MEMFREE(EM_S);
326     THREAD_MEMFREE(EM_F);
327     THREAD_MEMFREE(row_index);
328 gross 798
329 gross 853 } /* end of pointer check */
330 gross 798 } /* end parallel region */
331     }
332     /*
333     * $Log$
334     */

  ViewVC Help
Powered by ViewVC 1.1.26