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

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

  ViewVC Help
Powered by ViewVC 1.1.26