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

Contents of /branches/domexper/dudley/src/Assemble_PDE_Single2_1D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (show annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 3 months 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 /*
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 #ifdef _OPENMP
45 #include <omp.h>
46 #endif
47
48
49 /**************************************************************/
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 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 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 #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 {
78 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 {
92 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 }
152 }
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 rtmp=0;
214 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 }
217 } else {
218 for (s=0;s<p.row_NS;s++) {
219 rtmp=0;
220 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 }
223 }
224 }
225 /**************************************************************/
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 rtmp=0;
234 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 }
237 } else {
238 for (s=0;s<p.row_NS;s++) {
239 rtmp=0;
240 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 }
243 }
244 }
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
260 } /* end of pointer check */
261 } /* end parallel region */
262 }
263 /*
264 * $Log$
265 */

  ViewVC Help
Powered by ViewVC 1.1.26