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

Annotation of /branches/domexper/dudley/src/Assemble_PDE_Single2_2D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 4 months ago) by ksteube
Original Path: trunk/finley/src/Assemble_PDE_Single2_2D.c
File MIME type: text/plain
File size: 14346 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
23    
24     /* 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 = 2 x 2 */
30     /* B = 2 */
31     /* C = 2 */
32     /* D = scalar */
33     /* X = 2 */
34     /* Y = scalar */
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     /**************************************************************/
47    
48     void Finley_Assemble_PDE_Single2_2D(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 2
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;
58     register double rtmp00, rtmp01, rtmp10, rtmp11, rtmp, rtmp0, rtmp1;
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;
70     dim_t len_EM_F=p.row_NN;
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,rtmp00, rtmp01, rtmp10, rtmp11, rtmp, rtmp0, rtmp1,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    
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     rtmp=0;
102     for (q=0;q<p.numQuad;q++) {
103     rtmp+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
104     + DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*A_p[INDEX3(0,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]
105     + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX3(1,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
106     + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*A_p[INDEX3(1,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)] );
107     }
108     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
109     }
110     }
111     } else {
112     for (s=0;s<p.row_NS;s++) {
113     for (r=0;r<p.col_NS;r++) {
114     rtmp00=0;
115     rtmp01=0;
116     rtmp10=0;
117     rtmp11=0;
118     for (q=0;q<p.numQuad;q++) {
119     rtmp0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
120     rtmp1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
121     rtmp00+=rtmp0*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
122     rtmp01+=rtmp0*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
123     rtmp10+=rtmp1*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
124     rtmp11+=rtmp1*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
125     }
126     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+= rtmp00*A_p[INDEX2(0,0,DIM)]
127     + rtmp01*A_p[INDEX2(0,1,DIM)]
128     + rtmp10*A_p[INDEX2(1,0,DIM)]
129     + rtmp11*A_p[INDEX2(1,1,DIM)];
130     }
131     }
132     }
133     }
134     /**************************************************************/
135     /* process B: */
136     /**************************************************************/
137     B_p=getSampleData(B,e);
138     if (NULL!=B_p) {
139     add_EM_S=TRUE;
140     if (extendedB) {
141     for (s=0;s<p.row_NS;s++) {
142     for (r=0;r<p.col_NS;r++) {
143     rtmp=0.;
144     for (q=0;q<p.numQuad;q++) {
145     rtmp+=Vol[q]*S[INDEX2(r,q,p.row_NS)]*(DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX2(0,q,DIM)]
146     + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*B_p[INDEX2(1,q,DIM)]);
147     }
148     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
149     }
150     }
151     } else {
152     for (s=0;s<p.row_NS;s++) {
153     for (r=0;r<p.col_NS;r++) {
154     rtmp0=0;
155     rtmp1=0;
156     for (q=0;q<p.numQuad;q++) {
157     rtmp=Vol[q]*S[INDEX2(r,q,p.row_NS)];
158     rtmp0+=rtmp*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
159     rtmp1+=rtmp*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
160     }
161     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*B_p[0]+rtmp1*B_p[1];
162     }
163     }
164 gross 798 }
165 gross 853 }
166     /**************************************************************/
167     /* process C: */
168     /**************************************************************/
169     C_p=getSampleData(C,e);
170     if (NULL!=C_p) {
171     add_EM_S=TRUE;
172     if (extendedC) {
173     for (s=0;s<p.row_NS;s++) {
174     for (r=0;r<p.col_NS;r++) {
175     rtmp=0;
176     for (q=0;q<p.numQuad;q++) {
177     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*(C_p[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_NN,DIM)]
178     + C_p[INDEX2(1,q,DIM)]*DSDX[INDEX3(r,1,q,p.row_NN,DIM)]);
179     }
180     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
181     }
182     }
183     } else {
184     for (s=0;s<p.row_NS;s++) {
185     for (r=0;r<p.col_NS;r++) {
186     rtmp0=0;
187     rtmp1=0;
188     for (q=0;q<p.numQuad;q++) {
189     rtmp=Vol[q]*S[INDEX2(s,q,p.row_NS)];
190     rtmp0+=rtmp*DSDX[INDEX3(r,0,q,p.row_NN,DIM)];
191     rtmp1+=rtmp*DSDX[INDEX3(r,1,q,p.row_NN,DIM)];
192     }
193     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp0*C_p[0]+rtmp1*C_p[1];
194     }
195     }
196     }
197     }
198     /************************************************************* */
199     /* process D */
200     /**************************************************************/
201     D_p=getSampleData(D,e);
202     if (NULL!=D_p) {
203     add_EM_S=TRUE;
204     if (extendedD) {
205     for (s=0;s<p.row_NS;s++) {
206     for (r=0;r<p.col_NS;r++) {
207     rtmp=0;
208     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q]*S[INDEX2(r,q,p.row_NS)];
209     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
210     }
211     }
212     } else {
213     for (s=0;s<p.row_NS;s++) {
214     for (r=0;r<p.col_NS;r++) {
215     rtmp=0;
216     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
217     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[0];
218     }
219     }
220     }
221     }
222     /**************************************************************/
223     /* process X: */
224     /**************************************************************/
225     X_p=getSampleData(X,e);
226     if (NULL!=X_p) {
227     add_EM_F=TRUE;
228     if (extendedX) {
229     for (s=0;s<p.row_NS;s++) {
230     rtmp=0.;
231 gross 798 for (q=0;q<p.numQuad;q++) {
232 gross 853 rtmp+=Vol[q]*( DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX2(0,q,DIM)]
233     + DSDX[INDEX3(s,1,q,p.row_NN,DIM)]*X_p[INDEX2(1,q,DIM)]);
234     }
235     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;
236 gross 798 }
237 gross 853 } else {
238     for (s=0;s<p.row_NS;s++) {
239     rtmp0=0.;
240     rtmp1=0.;
241 gross 798 for (q=0;q<p.numQuad;q++) {
242 gross 853 rtmp0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
243     rtmp1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_NN,DIM)];
244 gross 798 }
245 gross 853 EM_F[INDEX2(0,s,p.numEqu)]+=rtmp0*X_p[0]+rtmp1*X_p[1];
246 gross 798 }
247 gross 853 }
248 gross 798 }
249 gross 853 /**************************************************************/
250     /* process Y: */
251     /**************************************************************/
252     Y_p=getSampleData(Y,e);
253     if (NULL!=Y_p) {
254     add_EM_F=TRUE;
255     if (extendedY) {
256     for (s=0;s<p.row_NS;s++) {
257     rtmp=0;
258     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[q];
259     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;
260 gross 798 }
261 gross 853 } else {
262     for (s=0;s<p.row_NS;s++) {
263 gross 798 rtmp=0;
264 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
265     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];
266 gross 798 }
267     }
268 gross 853 }
269     /***********************************************************************************************/
270     /* add the element matrices onto the matrix and right hand side */
271     /***********************************************************************************************/
272     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
273     if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
274     if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
275    
276     } /* end color check */
277     } /* end element loop */
278     } /* end color loop */
279    
280     THREAD_MEMFREE(EM_S);
281     THREAD_MEMFREE(EM_F);
282     THREAD_MEMFREE(row_index);
283 gross 798
284 gross 853 } /* end of pointer check */
285 gross 798 } /* end parallel region */
286     }
287     /*
288     * $Log$
289     */

  ViewVC Help
Powered by ViewVC 1.1.26