/[escript]/trunk/ripley/src/Assemble_PDE_System_3D_reduced.c
ViewVC logotype

Contents of /trunk/ripley/src/Assemble_PDE_System_3D_reduced.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3641 - (show annotations)
Thu Oct 27 02:16:12 2011 UTC (7 years, 4 months ago) by gross
File MIME type: text/plain
File size: 18199 byte(s)
more work on the assemblage
1 /*******************************************************
2 *
3 * Copyright (c) 2003-2010 by University of Queensland
4 * Earth Systems Science Computational Center (ESSCC)
5 * http://www.uq.edu.au/esscc
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 *******************************************************/
12 /**************************************************************/
13 /* */
14 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
15 /* for coefficients on the reduced integration scheme */
16 /* */
17 /* -(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 */
18 /* */
19 /* u has p.numComp components in a 2D domain. The shape functions for test and solution must be identical */
20 /* */
21 /* A = p.numEqu x DIM x p.numComp x DIM */
22 /* B = p.numEqu x p.numComp x DIM */
23 /* C = p.numEqu x DIM x p.numComp */
24 /* D = p.numEqu x p.numComp */
25 /* X = p.numEqu x DIM */
26 /* Y = p.numEqu */
27 /* */
28 /**************************************************************/
29 #include "Assemble.h"
30 #include "Util.h"
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34 /**************************************************************/
35 void Finley_Assemble_PDE_System_3D_reduced(Finley_Assemble_Parameters p,
36 Finley_ElementFile* elements,
37 Paso_SystemMatrix* Mat, escriptDataC* F,
38 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
39 #define DIM 3
40 index_t color;
41 dim_t e, isub;
42 __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
43 double *EM_S, *EM_F, *Vol, *DSDX;
44 index_t *row_index;
45 register dim_t q, s,r,k,m;
46 register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
47 bool_t add_EM_F, add_EM_S;
48 bool_t extendedA=isExpanded(A);
49 bool_t extendedB=isExpanded(B);
50 bool_t extendedC=isExpanded(C);
51 bool_t extendedD=isExpanded(D);
52 bool_t extendedX=isExpanded(X);
53 bool_t extendedY=isExpanded(Y);
54 double *F_p=(requireWrite(F), getSampleDataRW(F,0)); /* use comma, to get around the mixed code and declarations thing */
55 dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
56 dim_t len_EM_F=p.row_numShapesTotal*p.numEqu;
57 {
58 /* GENERATOR SNIP_PRE TOP */
59 const double w12 = -0.5*h0;
60 const double w2 = -0.25000000000000000000;
61 const double w1 = 0.25000000000000000000;
62 const double w0 = -0.25*h1/h0;
63 const double w15 = 0.25*h0*h1;
64 const double w7 = -0.125*h0;
65 const double w9 = 0.125*h0;
66 const double w14 = 0.5*h0;
67 const double w11 = -0.5*h1;
68 const double w4 = -0.25*h0/h1;
69 const double w10 = 0.0625*h0*h1;
70 const double w13 = 0.5*h1;
71 const double w8 = 0.125*h1;
72 const double w6 = -0.125*h1;
73 const double w3 = 0.25*h0/h1;
74 const double w5 = 0.25*h1/h0;
75 /* GENERATOR SNIP_PRE BOTTOM */
76 #pragma omp parallel private(EM_S, EM_F,k2_0, k0, k1, k2, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p)
77 {
78 EM_S=THREAD_MEMALLOC(len_EM_S,double);
79 EM_F=THREAD_MEMALLOC(len_EM_F,double);
80 if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
81 for (k2_0 = 0; k2_0 <2; k2_0++) { /* coloring */
82 #pragma omp parallel for private(i2, i1,i0)
83 for (k2 = k2_0; k2< N0; k2=k2+2) {
84 for (k1 = 0; k1< N1; ++k1) {
85 for (k0 = 0; k0< N0; ++k0) {
86 bool_t add_EM_F=FALSE;
87 bool_t add_EM_S=FALSE;
88 index_t e = k0 + M0 * k1 + M0*M1 * k2;
89 A_p=getSampleDataRO(A,e);
90 B_p=getSampleDataRO(B,e);
91 C_p=getSampleDataRO(C,e);
92 D_p=getSampleDataRO(D,e);
93 X_p=getSampleDataRO(X,e);
94 Y_p=getSampleDataRO(Y,e);
95 for (q=0;q<len_EM_S;++q) EM_S[q]=0;
96 for (q=0;q<len_EM_F;++q) EM_F[q]=0;
97 /* GENERATOR SNIP TOP */
98 /**************************************************************/
99 /* process A: */
100 /**************************************************************/
101 if (NULL!=A_p) {
102 add_EM_S=TRUE;
103 for (k=0;k<p.numEqu;k++) {
104 for (m=0;m<p.numComp;m++) {
105 const register double A_00 = A_p[INDEX4(k,0,m,0 p.numEqu,2, p.numComp)];
106 const register double A_01 = A_p[INDEX4(k,0,m,1 p.numEqu,2, p.numComp)];
107 const register double A_10 = A_p[INDEX4(k,1,m,0 p.numEqu,2, p.numComp)];
108 const register double A_11 = A_p[INDEX4(k,1,m,1 p.numEqu,2, p.numComp)];
109 const register double tmp0_0 = A_01 + A_10;
110 const register double tmp2_1 = A_01*w1;
111 const register double tmp7_1 = tmp0_0*w2;
112 const register double tmp3_1 = A_10*w2;
113 const register double tmp6_1 = A_00*w5;
114 const register double tmp1_1 = A_00*w0;
115 const register double tmp4_1 = tmp0_0*w1;
116 const register double tmp9_1 = A_01*w2;
117 const register double tmp5_1 = A_11*w4;
118 const register double tmp8_1 = A_10*w1;
119 const register double tmp0_1 = A_11*w3;
120 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
121 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
122 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
123 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
124 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp4_1 + tmp6_1;
125 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
126 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
127 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp4_1 + tmp5_1;
128 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp5_1 + tmp6_1 + tmp8_1 + tmp9_1;
129 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
130 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1 + tmp5_1 + tmp6_1;
131 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
132 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
133 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp8_1 + tmp9_1;
134 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp5_1 + tmp7_1;
135 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp6_1 + tmp7_1;
136 }
137 }
138 }
139 /**************************************************************/
140 /* process B: */
141 /**************************************************************/
142 if (NULL!=B_p) {
143 add_EM_S=TRUE;
144 for (k=0;k<p.numEqu;k++) {
145 for (m=0;m<p.numComp;m++) {
146 const register double B_0 = B_p[INDEX3(k,0,m, p.numEqu,2)];
147 const register double B_1 = B_p[INDEX3(k,1,m, p.numEqu,2)];
148 const register double tmp2_1 = B_0*w8;
149 const register double tmp0_1 = B_0*w6;
150 const register double tmp3_1 = B_1*w9;
151 const register double tmp1_1 = B_1*w7;
152 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
153 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp2_1;
154 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
155 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
156 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
157 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
158 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
159 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp3_1;
160 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
161 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp3_1;
162 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp2_1;
163 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp3_1;
164 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp3_1;
165 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp2_1;
166 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
167 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp2_1;
168 }
169 }
170 }
171 /**************************************************************/
172 /* process C: */
173 /**************************************************************/
174 if (NULL!=C_p) {
175 add_EM_S=TRUE;
176 for (k=0;k<p.numEqu;k++) {
177 for (m=0;m<p.numComp;m++) {
178 const register double C_0 = C_p[INDEX3(k,m,0, p.numEqu,p.numComp)];
179 const register double C_1 = C_p[INDEX3(k,m,1, p.numEqu,p.numComp)];
180 const register double tmp1_1 = C_1*w7;
181 const register double tmp0_1 = C_0*w8;
182 const register double tmp3_1 = C_0*w6;
183 const register double tmp2_1 = C_1*w9;
184 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
185 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
186 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
187 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp3_1;
188 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp2_1;
189 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp3_1;
190 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
191 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
192 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
193 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp3_1;
194 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp2_1;
195 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp2_1;
196 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
197 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp3_1;
198 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp2_1;
199 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
200 }
201 }
202 }
203 /**************************************************************/
204 /* process D: */
205 /**************************************************************/
206 if (NULL!=D_p) {
207 add_EM_S=TRUE;
208 for (k=0;k<p.numEqu;k++) {
209 for (m=0;m<p.numComp;m++) {
210 const register double D_0 = D_p[INDEX2(k,m, p.numEqu)];
211 const register double tmp0_1 = D_0*w10;
212 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1;
213 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp0_1;
214 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp0_1;
215 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1;
216 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp0_1;
217 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp0_1;
218 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp0_1;
219 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp0_1;
220 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp0_1;
221 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp0_1;
222 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp0_1;
223 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1;
224 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp0_1;
225 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1;
226 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp0_1;
227 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1;
228 }
229 }
230 }
231 /**************************************************************/
232 /* process X: */
233 /**************************************************************/
234 if (NULL!=X_p) {
235 add_EM_F=TRUE;
236 for (k=0;k<p.numEqu;k++) {
237 const register double X_0 = X_p[INDEX2(k,0, p.numEqu)];
238 const register double X_1 = X_p[INDEX2(k,1, p.numEqu)];
239 const register double tmp0_1 = X_0*w11;
240 const register double tmp2_1 = X_0*w13;
241 const register double tmp1_1 = X_1*w12;
242 const register double tmp3_1 = X_1*w14;
243 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1 + tmp1_1;
244 EM_F[INDEX2(k,1,p.numEqu)]+=tmp1_1 + tmp2_1;
245 EM_F[INDEX2(k,2,p.numEqu)]+=tmp0_1 + tmp3_1;
246 EM_F[INDEX2(k,3,p.numEqu)]+=tmp2_1 + tmp3_1;
247 }
248 }
249 /**************************************************************/
250 /* process Y: */
251 /**************************************************************/
252 if (NULL!=Y_p) {
253 add_EM_F=TRUE;
254 for (k=0;k<p.numEqu;k++) {
255 const register double Y_0 = Y_p[k];
256 const register double tmp0_1 = Y_0*w15;
257 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1;
258 EM_F[INDEX2(k,1,p.numEqu)]+=tmp0_1;
259 EM_F[INDEX2(k,2,p.numEqu)]+=tmp0_1;
260 EM_F[INDEX2(k,3,p.numEqu)]+=tmp0_1;
261 }
262 }
263 /* GENERATOR SNIP BOTTOM */
264 /* add to matrix (if add_EM_S) and RHS (if add_EM_F)*/
265 } /* end k0 */
266 } /* end k1 */
267 } /* end k2 */
268 } /* end coloring */
269 } /* end of pointer check */
270 THREAD_MEMFREE(EM_S);
271 THREAD_MEMFREE(EM_F);
272 } /* end parallel region */
273 } /* end of presettings */
274 }

  ViewVC Help
Powered by ViewVC 1.1.26