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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 gross 798
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 ksteube 1312
14 ksteube 1811
15 gross 798 /**************************************************************/
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 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 = p.numEqu x 1 x p.numComp x 1 */
29     /* B = 1 x numEqu x p.numComp */
30     /* C = p.numEqu x 1 x p.numComp */
31     /* D = p.numEqu x p.numComp */
32     /* X = p.numEqu x 1 */
33     /* Y = p.numEqu */
34    
35    
36     /**************************************************************/
37    
38    
39     #include "Assemble.h"
40     #include "Util.h"
41 gross 853 #ifdef _OPENMP
42     #include <omp.h>
43     #endif
44 gross 798
45 gross 853
46 gross 798 /**************************************************************/
47    
48     void Finley_Assemble_PDE_System2_1D(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 1
53     index_t color;
54     dim_t e;
55 gross 853 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;
59     bool_t add_EM_F, add_EM_S;
60    
61 gross 798 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 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,add_EM_F, add_EM_S)
74 gross 798 {
75 gross 853 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 gross 798
81 gross 853 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]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*
106     A_p[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
107     }
108     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
109     }
110     }
111     }
112     }
113     } else {
114     for (s=0;s<p.row_NS;s++) {
115     for (r=0;r<p.col_NS;r++) {
116     rtmp=0;
117     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
118     for (k=0;k<p.numEqu;k++) {
119     for (m=0;m<p.numComp;m++) {
120     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)];
121     }
122     }
123     }
124     }
125     }
126     }
127     /**************************************************************/
128     /* process B: */
129     /**************************************************************/
130     B_p=getSampleData(B,e);
131     if (NULL!=B_p) {
132     add_EM_S=TRUE;
133     if (extendedB) {
134     for (s=0;s<p.row_NS;s++) {
135     for (r=0;r<p.col_NS;r++) {
136     for (k=0;k<p.numEqu;k++) {
137     for (m=0;m<p.numComp;m++) {
138     rtmp=0.;
139     for (q=0;q<p.numQuad;q++) {
140     rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*
141     B_p[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]*S[INDEX2(r,q,p.row_NS)];
142     }
143     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
144     }
145     }
146     }
147     }
148     } else {
149     for (s=0;s<p.row_NS;s++) {
150     for (r=0;r<p.col_NS;r++) {
151     rtmp=0;
152     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*S[INDEX2(r,q,p.row_NS)];
153     for (k=0;k<p.numEqu;k++) {
154     for (m=0;m<p.numComp;m++) {
155     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[INDEX3(k,0,m,p.numEqu,DIM)];
156     }
157     }
158     }
159     }
160 gross 798 }
161 gross 853 }
162     /**************************************************************/
163     /* process C: */
164     /**************************************************************/
165     C_p=getSampleData(C,e);
166     if (NULL!=C_p) {
167     add_EM_S=TRUE;
168     if (extendedC) {
169     for (s=0;s<p.row_NS;s++) {
170     for (r=0;r<p.col_NS;r++) {
171     for (k=0;k<p.numEqu;k++) {
172     for (m=0;m<p.numComp;m++) {
173     rtmp=0;
174     for (q=0;q<p.numQuad;q++) {
175     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*C_p[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
176     }
177     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
178     }
179     }
180     }
181     }
182     } else {
183     for (s=0;s<p.row_NS;s++) {
184     for (r=0;r<p.col_NS;r++) {
185     rtmp=0;
186     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
187     for (k=0;k<p.numEqu;k++) {
188 gross 798 for (m=0;m<p.numComp;m++) {
189 gross 853 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)];
190     }
191     }
192     }
193     }
194     }
195     }
196     /************************************************************* */
197     /* process D */
198     /**************************************************************/
199     D_p=getSampleData(D,e);
200     if (NULL!=D_p) {
201     add_EM_S=TRUE;
202     if (extendedD) {
203     for (s=0;s<p.row_NS;s++) {
204     for (r=0;r<p.col_NS;r++) {
205     for (k=0;k<p.numEqu;k++) {
206     for (m=0;m<p.numComp;m++) {
207     rtmp=0;
208     for (q=0;q<p.numQuad;q++) {
209     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)];
210     }
211     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
212    
213 gross 798 }
214 gross 853 }
215     }
216     }
217     } else {
218     for (s=0;s<p.row_NS;s++) {
219     for (r=0;r<p.col_NS;r++) {
220     rtmp=0;
221     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
222     for (k=0;k<p.numEqu;k++) {
223     for (m=0;m<p.numComp;m++) {
224     EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[INDEX2(k,m,p.numEqu)];
225     }
226     }
227     }
228     }
229     }
230 gross 798 }
231 gross 853 /**************************************************************/
232     /* process X: */
233     /**************************************************************/
234     X_p=getSampleData(X,e);
235     if (NULL!=X_p) {
236     add_EM_F=TRUE;
237     if (extendedX) {
238     for (s=0;s<p.row_NS;s++) {
239     for (k=0;k<p.numEqu;k++) {
240 gross 798 rtmp=0;
241 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX3(k,0,q,p.numEqu,DIM)];
242     EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
243     }
244 gross 798 }
245 gross 853 } else {
246     for (s=0;s<p.row_NS;s++) {
247     rtmp=0;
248     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
249     for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*X_p[INDEX2(k,0,p.numEqu)];
250 gross 798 }
251 gross 853 }
252 gross 798 }
253 gross 853 /**************************************************************/
254     /* process Y: */
255     /**************************************************************/
256     Y_p=getSampleData(Y,e);
257     if (NULL!=Y_p) {
258     add_EM_F=TRUE;
259     if (extendedY) {
260     for (s=0;s<p.row_NS;s++) {
261     for (k=0;k<p.numEqu;k++) {
262 gross 798 rtmp=0;
263 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
264     EM_F[INDEX2(k,s,p.numEqu)]+=rtmp;
265     }
266 gross 798 }
267 gross 853 } else {
268     for (s=0;s<p.row_NS;s++) {
269 gross 798 rtmp=0;
270 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
271     for (k=0;k<p.numEqu;k++) EM_F[INDEX2(k,s,p.numEqu)]+=rtmp*Y_p[k];
272 gross 798 }
273     }
274 gross 853 }
275     /***********************************************************************************************/
276     /* add the element matrices onto the matrix and right hand side */
277     /***********************************************************************************************/
278     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
279     if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
280     if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
281    
282     } /* end color check */
283     } /* end element loop */
284     } /* end color loop */
285    
286     THREAD_MEMFREE(EM_S);
287     THREAD_MEMFREE(EM_F);
288     THREAD_MEMFREE(row_index);
289 gross 798
290 gross 853 } /* end of pointer check */
291 gross 798 } /* end parallel region */
292     }
293     /*
294     * $Log$
295     */

  ViewVC Help
Powered by ViewVC 1.1.26