/[escript]/branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp
ViewVC logotype

Annotation of /branches/doubleplusgood/dudley/src/Assemble_PDE_Single2_1D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (hide annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 1 month ago) by ksteube
Original Path: trunk/finley/src/Assemble_PDE_Single2_1D.c
File MIME type: text/plain
File size: 12039 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


1 gross 798 /*
2     ************************************************************
3     * Copyright 2006 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 1.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11     */
12    
13     /**************************************************************/
14    
15     /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
16     /* the shape functions for test and solution must be identical */
17    
18    
19     /* -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y */
20    
21     /* 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 = 1 x 1 */
27     /* B = 1 */
28     /* C = 1 */
29     /* D = scalar */
30     /* X = 1 */
31     /* Y = scalar */
32    
33    
34     /**************************************************************/
35    
36     /* Author: gross@access.edu.au */
37     /* Version: $Id:$ */
38    
39     /**************************************************************/
40    
41    
42     #include "Assemble.h"
43     #include "Util.h"
44 gross 853 #ifdef _OPENMP
45     #include <omp.h>
46     #endif
47 gross 798
48 gross 853
49 gross 798 /**************************************************************/
50    
51     void Finley_Assemble_PDE_Single2_1D(Assemble_Parameters p, Finley_ElementFile* elements,
52     Paso_SystemMatrix* Mat, escriptDataC* F,
53     escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
54    
55     #define DIM 1
56     index_t color;
57     dim_t e;
58 gross 853 double *EM_S, *EM_F, *Vol, *DSDX, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p;
59     index_t *row_index;
60     register dim_t q, s,r;
61     register double rtmp;
62     bool_t add_EM_F, add_EM_S;
63    
64 gross 798 bool_t extendedA=isExpanded(A);
65     bool_t extendedB=isExpanded(B);
66     bool_t extendedC=isExpanded(C);
67     bool_t extendedD=isExpanded(D);
68     bool_t extendedX=isExpanded(X);
69     bool_t extendedY=isExpanded(Y);
70     double *F_p=getSampleData(F,0);
71     double *S=p.row_jac->ReferenceElement->S;
72     dim_t len_EM_S=p.row_NN*p.col_NN;
73     dim_t len_EM_F=p.row_NN;
74    
75    
76 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,rtmp,add_EM_F, add_EM_S)
77 gross 798 {
78 gross 853 EM_S=THREAD_MEMALLOC(len_EM_S,double);
79     EM_F=THREAD_MEMALLOC(len_EM_F,double);
80     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
81    
82     if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
83    
84     #ifndef PASO_MPI
85     for (color=elements->minColor;color<=elements->maxColor;color++) {
86     /* open loop over all elements: */
87     #pragma omp for private(e) schedule(static)
88     for(e=0;e<elements->numElements;e++){
89     if (elements->Color[e]==color) {
90     #else
91 ksteube 971 {
92 gross 853 for(e=0;e<elements->numElements;e++) {
93     {
94     #endif
95     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
96     DSDX=&(p.row_jac->DSDX[INDEX4(0,0,0,e,p.row_NN,DIM,p.numQuad)]);
97     for (q=0;q<len_EM_S;++q) EM_S[q]=0;
98     for (q=0;q<len_EM_F;++q) EM_F[q]=0;
99     add_EM_F=FALSE;
100     add_EM_S=FALSE;
101     /**************************************************************/
102     /* process A: */
103     /**************************************************************/
104     A_p=getSampleData(A,e);
105     if (NULL!=A_p) {
106     add_EM_S=TRUE;
107     if (extendedA) {
108     for (s=0;s<p.row_NS;s++) {
109     for (r=0;r<p.col_NS;r++) {
110     rtmp=0;
111     for (q=0;q<p.numQuad;q++) {
112     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)];
113     }
114     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
115     }
116     }
117     } else {
118     for (s=0;s<p.row_NS;s++) {
119     for (r=0;r<p.col_NS;r++) {
120     rtmp=0;
121     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)];
122     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*A_p[INDEX2(0,0,DIM)];
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     rtmp=0;
137     for (q=0;q<p.numQuad;q++) {
138     rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*B_p[INDEX2(0,q,DIM)]*S[INDEX2(r,q,p.row_NS)];
139     }
140     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
141     }
142     }
143     } else {
144     for (s=0;s<p.row_NS;s++) {
145     for (r=0;r<p.col_NS;r++) {
146     rtmp=0;
147     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)];
148     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*B_p[0];
149     }
150     }
151 gross 798 }
152 gross 853 }
153     /**************************************************************/
154     /* process C: */
155     /**************************************************************/
156     C_p=getSampleData(C,e);
157     if (NULL!=C_p) {
158     add_EM_S=TRUE;
159     if (extendedC) {
160     for (s=0;s<p.row_NS;s++) {
161     for (r=0;r<p.col_NS;r++) {
162     rtmp=0;
163     for (q=0;q<p.numQuad;q++) {
164     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)];
165     }
166     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
167     }
168     }
169     } else {
170     for (s=0;s<p.row_NS;s++) {
171     for (r=0;r<p.col_NS;r++) {
172     rtmp=0;
173     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)];
174     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*C_p[0];
175     }
176     }
177     }
178     }
179     /************************************************************* */
180     /* process D */
181     /**************************************************************/
182     D_p=getSampleData(D,e);
183     if (NULL!=D_p) {
184     add_EM_S=TRUE;
185     if (extendedD) {
186     for (s=0;s<p.row_NS;s++) {
187     for (r=0;r<p.col_NS;r++) {
188     rtmp=0;
189     for (q=0;q<p.numQuad;q++) {
190     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q]*S[INDEX2(r,q,p.row_NS)];
191     }
192     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp;
193     }
194     }
195     } else {
196     for (s=0;s<p.row_NS;s++) {
197     for (r=0;r<p.col_NS;r++) {
198     rtmp=0;
199     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
200     EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_NN)]+=rtmp*D_p[0];
201     }
202     }
203     }
204     }
205     /**************************************************************/
206     /* process X: */
207     /**************************************************************/
208     X_p=getSampleData(X,e);
209     if (NULL!=X_p) {
210     add_EM_F=TRUE;
211     if (extendedX) {
212     for (s=0;s<p.row_NS;s++) {
213 gross 798 rtmp=0;
214 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)]*X_p[INDEX2(0,q,DIM)];
215     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp;
216 gross 798 }
217 gross 853 } else {
218     for (s=0;s<p.row_NS;s++) {
219 gross 798 rtmp=0;
220 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_NN,DIM)];
221     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*X_p[0];
222 gross 798 }
223 gross 853 }
224 gross 798 }
225 gross 853 /**************************************************************/
226     /* process Y: */
227     /**************************************************************/
228     Y_p=getSampleData(Y,e);
229     if (NULL!=Y_p) {
230     add_EM_F=TRUE;
231     if (extendedY) {
232     for (s=0;s<p.row_NS;s++) {
233 gross 798 rtmp=0;
234 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[q];
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 gross 798 rtmp=0;
240 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
241     EM_F[INDEX2(0,s,p.numEqu)]+=rtmp*Y_p[0];
242 gross 798 }
243     }
244 gross 853 }
245     /***********************************************************************************************/
246     /* add the element matrices onto the matrix and right hand side */
247     /***********************************************************************************************/
248     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
249     if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
250     if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
251    
252     } /* end color check */
253     } /* end element loop */
254     } /* end color loop */
255    
256     THREAD_MEMFREE(EM_S);
257     THREAD_MEMFREE(EM_F);
258     THREAD_MEMFREE(row_index);
259 gross 798
260 gross 853 } /* end of pointer check */
261 gross 798 } /* end parallel region */
262     }
263     /*
264     * $Log$
265     */

  ViewVC Help
Powered by ViewVC 1.1.26