/[escript]/trunk/finley/src/Assemble_PDE_System2_3D.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_PDE_System2_3D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4442 - (show annotations)
Fri Jun 7 07:04:44 2013 UTC (6 years, 1 month ago) by caltinay
File size: 23475 byte(s)
cosmetic changes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /****************************************************************************
18
19 Assembles the system of numEqu PDEs into the stiffness matrix S and right
20 hand side F
21
22 -(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
23 and
24 -(X_{k,i})_i + Y_k
25
26 u has p.numComp components in a 3D domain. The shape functions for test and
27 solution must be identical and row_NS == row_NN.
28
29 Shape of the coefficients:
30
31 A = p.numEqu x 3 x p.numComp x 3
32 B = 3 x p.numEqu x p.numComp
33 C = p.numEqu x 3 x p.numComp
34 D = p.numEqu x p.numComp
35 X = p.numEqu x 3
36 Y = p.numEqu
37
38 *****************************************************************************/
39
40 #include "Assemble.h"
41 #include "Util.h"
42
43 namespace finley {
44
45 void Assemble_PDE_System_3D(const AssembleParameters& p,
46 escript::Data& A, escript::Data& B,
47 escript::Data& C, escript::Data& D,
48 escript::Data& X, escript::Data& Y)
49 {
50 const int DIM = 3;
51 bool expandedA=A.actsExpanded();
52 bool expandedB=B.actsExpanded();
53 bool expandedC=C.actsExpanded();
54 bool expandedD=D.actsExpanded();
55 bool expandedX=X.actsExpanded();
56 bool expandedY=Y.actsExpanded();
57 p.F.requireWrite();
58 double *F_p=p.F.getSampleDataRW(0);
59 const double *S=p.row_jac->BasisFunctions->S;
60 const size_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
61 const size_t len_EM_F=p.row_numShapesTotal*p.numEqu;
62
63 #pragma omp parallel
64 {
65 for (int color=p.elements->minColor; color<=p.elements->maxColor; color++) {
66 // loop over all elements:
67 #pragma omp for
68 for (int e=0; e<p.elements->numElements; e++) {
69 if (p.elements->Color[e]==color) {
70 for (int isub=0; isub<p.numSub; isub++) {
71 const double *Vol=&(p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)]);
72 const double *DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,isub,e,p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)]);
73 std::vector<double> EM_S(len_EM_S);
74 std::vector<double> EM_F(len_EM_F);
75 bool add_EM_F=false;
76 bool add_EM_S=false;
77 ///////////////
78 // process A //
79 ///////////////
80 if (!A.isEmpty()) {
81 const double *A_p=A.getSampleDataRO(e);
82 add_EM_S=true;
83 if (expandedA) {
84 const double *A_q=&(A_p[INDEX6(0,0,0,0,0,isub,p.numEqu,DIM,p.numComp,DIM,p.numQuadSub)]);
85 for (int s=0; s<p.row_numShapes; s++) {
86 for (int r=0; r<p.col_numShapes; r++) {
87 for (int k=0; k<p.numEqu; k++) {
88 for (int m=0; m<p.numComp; m++) {
89 double f=0.;
90 for (int q=0; q<p.numQuadSub; q++) {
91 f+=Vol[q]*(
92 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
93 +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
94 +DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]
95 +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
96 +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
97 +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]
98 +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
99 +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
100 +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,2,m,2,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);
101 }
102 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
103 }
104 }
105 }
106 }
107 } else { // constant A
108 for (int s=0; s<p.row_numShapes; s++) {
109 for (int r=0; r<p.col_numShapes; r++) {
110 double f00=0;
111 double f01=0;
112 double f02=0;
113 double f10=0;
114 double f11=0;
115 double f12=0;
116 double f20=0;
117 double f21=0;
118 double f22=0;
119 for (int q=0; q<p.numQuadSub; q++) {
120 const double f0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
121 f00+=f0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
122 f01+=f0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
123 f02+=f0*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
124
125 const double f1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
126 f10+=f1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
127 f11+=f1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
128 f12+=f1*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
129
130 const double f2=Vol[q]*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
131 f20+=f2*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
132 f21+=f2*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
133 f22+=f2*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
134 }
135 for (int k=0; k<p.numEqu; k++) {
136 for (int m=0; m<p.numComp; m++) {
137 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
138 f00*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
139 +f01*A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
140 +f02*A_p[INDEX4(k,0,m,2,p.numEqu,DIM,p.numComp)]
141 +f10*A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
142 +f11*A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)]
143 +f12*A_p[INDEX4(k,1,m,2,p.numEqu,DIM,p.numComp)]
144 +f20*A_p[INDEX4(k,2,m,0,p.numEqu,DIM,p.numComp)]
145 +f21*A_p[INDEX4(k,2,m,1,p.numEqu,DIM,p.numComp)]
146 +f22*A_p[INDEX4(k,2,m,2,p.numEqu,DIM,p.numComp)];
147 }
148 }
149 }
150 }
151 }
152 }
153 ///////////////
154 // process B //
155 ///////////////
156 if (!B.isEmpty()) {
157 const double *B_p=B.getSampleDataRO(e);
158 add_EM_S=true;
159 if (expandedB) {
160 const double *B_q=&(B_p[INDEX5(0,0,0,0,isub,p.numEqu,DIM,p.numComp,p.numQuadSub)]);
161 for (int s=0; s<p.row_numShapes; s++) {
162 for (int r=0; r<p.col_numShapes; r++) {
163 for (int k=0; k<p.numEqu; k++) {
164 for (int m=0; m<p.numComp; m++) {
165 double f=0.;
166 for (int q=0; q<p.numQuadSub; q++) {
167 f+=Vol[q]*S[INDEX2(r,q,p.row_numShapes)]*(
168 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
169 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]
170 + DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,2,m,q,p.numEqu,DIM,p.numComp)]);
171 }
172 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
173 }
174 }
175 }
176 }
177 } else { // constant B
178 for (int s=0; s<p.row_numShapes; s++) {
179 for (int r=0; r<p.col_numShapes; r++) {
180 double f0=0;
181 double f1=0;
182 double f2=0;
183 for (int q=0; q<p.numQuadSub; q++) {
184 const double f=Vol[q]*S[INDEX2(r,q,p.row_numShapes)];
185 f0+=f*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
186 f1+=f*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
187 f2+=f*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
188 }
189 for (int k=0; k<p.numEqu; k++) {
190 for (int m=0; m<p.numComp; m++) {
191 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
192 f0*B_p[INDEX3(k,0,m,p.numEqu,DIM)]
193 +f1*B_p[INDEX3(k,1,m,p.numEqu,DIM)]
194 +f2*B_p[INDEX3(k,2,m,p.numEqu,DIM)];
195 }
196 }
197 }
198 }
199 }
200 }
201 ///////////////
202 // process C //
203 ///////////////
204 if (!C.isEmpty()) {
205 const double *C_p=C.getSampleDataRO(e);
206 add_EM_S=true;
207 if (expandedC) {
208 const double *C_q=&(C_p[INDEX5(0,0,0,0,isub,p.numEqu,p.numComp,DIM,p.numQuadSub)]);
209 for (int s=0; s<p.row_numShapes; s++) {
210 for (int r=0; r<p.col_numShapes; r++) {
211 for (int k=0; k<p.numEqu; k++) {
212 for (int m=0; m<p.numComp; m++) {
213 double f=0.;
214 for (int q=0; q<p.numQuadSub; q++) {
215 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*(
216 C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
217 +C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
218 +C_q[INDEX4(k,m,2,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)]);
219 }
220 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
221 }
222 }
223 }
224 }
225 } else { // constant C
226 for (int s=0; s<p.row_numShapes; s++) {
227 for (int r=0; r<p.col_numShapes; r++) {
228 double f0=0.;
229 double f1=0.;
230 double f2=0.;
231 for (int q=0; q<p.numQuadSub; q++) {
232 const double f=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
233 f0+=f*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
234 f1+=f*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
235 f2+=f*DSDX[INDEX3(r,2,q,p.row_numShapesTotal,DIM)];
236 }
237 for (int k=0; k<p.numEqu; k++) {
238 for (int m=0; m<p.numComp; m++) {
239 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
240 f0*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]
241 +f1*C_p[INDEX3(k,m,1,p.numEqu,p.numComp)]
242 +f2*C_p[INDEX3(k,m,2,p.numEqu,p.numComp)];
243 }
244 }
245 }
246 }
247 }
248 }
249 ///////////////
250 // process D //
251 ///////////////
252 if (!D.isEmpty()) {
253 const double *D_p=D.getSampleDataRO(e);
254 add_EM_S=true;
255 if (expandedD) {
256 const double *D_q=&(D_p[INDEX4(0,0,0,isub,p.numEqu,p.numComp,p.numQuadSub)]);
257 for (int s=0; s<p.row_numShapes; s++) {
258 for (int r=0; r<p.col_numShapes; r++) {
259 for (int k=0; k<p.numEqu; k++) {
260 for (int m=0; m<p.numComp; m++) {
261 double f=0.;
262 for (int q=0; q<p.numQuadSub; q++) {
263 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
264 }
265 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
266 }
267 }
268 }
269 }
270 } else { // constant D
271 for (int s=0; s<p.row_numShapes; s++) {
272 for (int r=0; r<p.col_numShapes; r++) {
273 double f=0.;
274 for (int q=0; q<p.numQuadSub; q++)
275 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
276 for (int k=0; k<p.numEqu; k++) {
277 for (int m=0; m<p.numComp; m++) {
278 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*D_p[INDEX2(k,m,p.numEqu)];
279 }
280 }
281 }
282 }
283 }
284 }
285 ///////////////
286 // process X //
287 ///////////////
288 if (!X.isEmpty()) {
289 const double *X_p=X.getSampleDataRO(e);
290 add_EM_F=true;
291 if (expandedX) {
292 const double *X_q=&(X_p[INDEX4(0,0,0,isub,p.numEqu,DIM,p.numQuadSub)]);
293 for (int s=0; s<p.row_numShapes; s++) {
294 for (int k=0; k<p.numEqu; k++) {
295 double f=0.;
296 for (int q=0; q<p.numQuadSub; q++) {
297 f+=Vol[q]*(
298 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]
299 +DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]
300 +DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,2,q,p.numEqu,DIM)]);
301 }
302 EM_F[INDEX2(k,s,p.numEqu)]+=f;
303 }
304 }
305 } else { // constant X
306 for (int s=0; s<p.row_numShapes; s++) {
307 double f0=0;
308 double f1=0;
309 double f2=0;
310 for (int q=0; q<p.numQuadSub; q++) {
311 f0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
312 f1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
313 f2+=Vol[q]*DSDX[INDEX3(s,2,q,p.row_numShapesTotal,DIM)];
314 }
315 for (int k=0; k<p.numEqu; k++) {
316 EM_F[INDEX2(k,s,p.numEqu)] +=
317 f0*X_p[INDEX2(k,0,p.numEqu)]
318 +f1*X_p[INDEX2(k,1,p.numEqu)]
319 +f2*X_p[INDEX2(k,2,p.numEqu)];
320 }
321 }
322 }
323 }
324 ///////////////
325 // process Y //
326 ///////////////
327 if (!Y.isEmpty()) {
328 const double *Y_p=Y.getSampleDataRO(e);
329 add_EM_F=true;
330 if (expandedY) {
331 const double *Y_q=&(Y_p[INDEX3(0,0,isub,p.numEqu,p.numQuadSub)]);
332 for (int s=0; s<p.row_numShapes; s++) {
333 for (int k=0; k<p.numEqu; k++) {
334 double f=0.;
335 for (int q=0; q<p.numQuadSub; q++)
336 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
337 EM_F[INDEX2(k,s,p.numEqu)]+=f;
338 }
339 }
340 } else { // constant Y
341 for (int s=0; s<p.row_numShapes; s++) {
342 double f=0.;
343 for (int q=0; q<p.numQuadSub; q++)
344 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
345 for (int k=0; k<p.numEqu; k++)
346 EM_F[INDEX2(k,s,p.numEqu)]+=f*Y_p[k];
347 }
348 }
349 }
350 // add the element matrices onto the matrix and
351 // right hand side
352 std::vector<int> row_index(p.row_numShapesTotal);
353 for (int q=0; q<p.row_numShapesTotal; q++)
354 row_index[q]=p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
355
356 if (add_EM_F)
357 util::addScatter(p.row_numShapesTotal,
358 &row_index[0], p.numEqu, &EM_F[0], F_p,
359 p.row_DOF_UpperBound);
360 if (add_EM_S)
361 Assemble_addToSystemMatrix(p.S,
362 p.row_numShapesTotal, &row_index[0],
363 p.numEqu, p.col_numShapesTotal,
364 &row_index[0], p.numComp, &EM_S[0]);
365 } // end of isub
366 } // end color check
367 } // end element loop
368 } // end color loop
369 } // end parallel region
370 }
371
372 } // namespace finley
373

Properties

Name Value
svn:mergeinfo /branches/lapack2681/finley/src/Assemble_PDE_System2_3D.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_PDE_System2_3D.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_PDE_System2_3D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_System2_3D.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_PDE_System2_3D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_System2_3D.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_PDE_System2_3D.cpp:3471-3974 /release/3.0/finley/src/Assemble_PDE_System2_3D.cpp:2591-2601 /trunk/finley/src/Assemble_PDE_System2_3D.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26