/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_PDE_System_3D.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_System_3D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 24500 byte(s)
Much needed sync with trunk...

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
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 dudley {
44
45 void Assemble_PDE_System_3D(const Assemble_Parameters& p, const Dudley_ElementFile* elements,
46 escript::ASM_ptr mat, escript::Data& F,
47 const escript::Data& A, const escript::Data& B,
48 const escript::Data& C, const escript::Data& D,
49 const escript::Data& X, const escript::Data& Y)
50 {
51 const int DIM = 3;
52 bool expandedA = A.actsExpanded();
53 bool expandedB = B.actsExpanded();
54 bool expandedC = C.actsExpanded();
55 bool expandedD = D.actsExpanded();
56 bool expandedX = X.actsExpanded();
57 bool expandedY = Y.actsExpanded();
58 double* F_p = NULL;
59 if (!F.isEmpty()) {
60 F.requireWrite();
61 F_p = F.getSampleDataRW(0);
62 }
63 const double* S = p.shapeFns;
64 const size_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp;
65 const size_t len_EM_F = p.numShapes * p.numEqu;
66
67 #pragma omp parallel
68 {
69 std::vector<double> EM_S(len_EM_S);
70 std::vector<double> EM_F(len_EM_F);
71 std::vector<index_t> row_index(p.numShapes);
72
73 for (int color = elements->minColor; color <= elements->maxColor; color++) {
74 // loop over all elements
75 #pragma omp for
76 for (index_t e = 0; e < elements->numElements; e++) {
77 if (elements->Color[e] == color) {
78 const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
79 const double* DSDX = &p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)];
80 std::fill(EM_S.begin(), EM_S.end(), 0);
81 std::fill(EM_F.begin(), EM_F.end(), 0);
82 bool add_EM_F = false;
83 bool add_EM_S = false;
84
85 ///////////////
86 // process A //
87 ///////////////
88 if (!A.isEmpty()) {
89 const double* A_p = A.getSampleDataRO(e);
90 add_EM_S = true;
91 if (expandedA) {
92 const double* A_q = &A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)];
93 for (int s = 0; s < p.numShapes; s++) {
94 for (int r = 0; r < p.numShapes; r++) {
95 for (int k = 0; k < p.numEqu; k++) {
96 for (int m = 0; m < p.numComp; m++) {
97 double f = 0;
98 for (int q = 0; q < p.numQuad; q++) {
99 f +=
100 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
101 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
102 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
103 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
104 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
105 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
106 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
107 A_q[INDEX5(k, 0, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
108 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
109 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
110 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
111 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
112 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
113 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
114 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
115 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
116 A_q[INDEX5(k, 1, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
117 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] +
118 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
119 A_q[INDEX5(k, 2, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
120 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
121 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
122 A_q[INDEX5(k, 2, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
123 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
124 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
125 A_q[INDEX5(k, 2, m, 2, q, p.numEqu, DIM, p.numComp, DIM)]
126 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
127
128 }
129 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
130 }
131 }
132 }
133 }
134 }
135 else
136 {
137 for (int s = 0; s < p.numShapes; s++) {
138 for (int r = 0; r < p.numShapes; r++) {
139 double f00 = 0;
140 double f01 = 0;
141 double f02 = 0;
142 double f10 = 0;
143 double f11 = 0;
144 double f12 = 0;
145 double f20 = 0;
146 double f21 = 0;
147 double f22 = 0;
148 for (int q = 0; q < p.numQuad; q++) {
149 const double f0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
150 f00 += f0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
151 f01 += f0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
152 f02 += f0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
153
154 const double f1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
155 f10 += f1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
156 f11 += f1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
157 f12 += f1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
158
159 const double f2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
160 f20 += f2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
161 f21 += f2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
162 f22 += f2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
163 }
164 for (int k = 0; k < p.numEqu; k++) {
165 for (int m = 0; m < p.numComp; m++) {
166 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
167 f00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] +
168 f01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] +
169 f02 * A_p[INDEX4(k, 0, m, 2, p.numEqu, DIM, p.numComp)] +
170 f10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)] +
171 f11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)] +
172 f12 * A_p[INDEX4(k, 1, m, 2, p.numEqu, DIM, p.numComp)] +
173 f20 * A_p[INDEX4(k, 2, m, 0, p.numEqu, DIM, p.numComp)] +
174 f21 * A_p[INDEX4(k, 2, m, 1, p.numEqu, DIM, p.numComp)] +
175 f22 * A_p[INDEX4(k, 2, m, 2, p.numEqu, DIM, p.numComp)];
176 }
177 }
178 }
179 }
180 }
181 }
182 ///////////////
183 // process B //
184 ///////////////
185 if (!B.isEmpty()) {
186 const double* B_p = B.getSampleDataRO(e);
187 add_EM_S = true;
188 if (expandedB) {
189 const double* B_q = &B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)];
190 for (int s = 0; s < p.numShapes; s++) {
191 for (int r = 0; r < p.numShapes; r++) {
192 for (int k = 0; k < p.numEqu; k++) {
193 for (int m = 0; m < p.numComp; m++) {
194 double f = 0;
195 for (int q = 0; q < p.numQuad; q++) {
196 f +=
197 vol * S[INDEX2(r, q, p.numShapes)] *
198 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
199 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
200 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
201 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)] +
202 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
203 B_q[INDEX4(k, 2, m, q, p.numEqu, DIM, p.numComp)]);
204 }
205 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
206 }
207 }
208 }
209 }
210 } else {
211 for (int s = 0; s < p.numShapes; s++) {
212 for (int r = 0; r < p.numShapes; r++) {
213 double f0 = 0;
214 double f1 = 0;
215 double f2 = 0;
216 for (int q = 0; q < p.numQuad; q++) {
217 const double f = vol * S[INDEX2(r, q, p.numShapes)];
218 f0 += f * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
219 f1 += f * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
220 f2 += f * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
221 }
222 for (int k = 0; k < p.numEqu; k++) {
223 for (int m = 0; m < p.numComp; m++) {
224 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
225 f0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
226 f1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)] +
227 f2 * B_p[INDEX3(k, 2, m, p.numEqu, DIM)];
228 }
229 }
230 }
231 }
232 }
233 }
234 ///////////////
235 // process C //
236 ///////////////
237 if (!C.isEmpty()) {
238 const double* C_p = C.getSampleDataRO(e);
239 add_EM_S = true;
240 if (expandedC) {
241 const double* C_q = &C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)];
242 for (int s = 0; s < p.numShapes; s++) {
243 for (int r = 0; r < p.numShapes; r++) {
244 for (int k = 0; k < p.numEqu; k++) {
245 for (int m = 0; m < p.numComp; m++) {
246 double f = 0;
247 for (int q = 0; q < p.numQuad; q++) {
248 f +=
249 vol * S[INDEX2(s, q, p.numShapes)] *
250 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
251 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
252 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
253 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
254 C_q[INDEX4(k, m, 2, q, p.numEqu, p.numComp, DIM)] *
255 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]);
256 }
257 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
258 }
259 }
260 }
261 }
262 } else {
263 for (int s = 0; s < p.numShapes; s++) {
264 for (int r = 0; r < p.numShapes; r++) {
265 double f0 = 0;
266 double f1 = 0;
267 double f2 = 0;
268 for (int q = 0; q < p.numQuad; q++) {
269 const double f = vol * S[INDEX2(s, q, p.numShapes)];
270 f0 += f * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
271 f1 += f * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
272 f2 += f * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)];
273 }
274 for (int k = 0; k < p.numEqu; k++) {
275 for (int m = 0; m < p.numComp; m++)
276 {
277 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
278 f0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
279 f1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)] +
280 f2 * C_p[INDEX3(k, m, 2, p.numEqu, p.numComp)];
281 }
282 }
283 }
284 }
285 }
286 }
287 ///////////////
288 // process D //
289 ///////////////
290 if (!D.isEmpty()) {
291 const double* D_p = D.getSampleDataRO(e);
292 add_EM_S = true;
293 if (expandedD) {
294 const double* D_q = &D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)];
295 for (int s = 0; s < p.numShapes; s++) {
296 for (int r = 0; r < p.numShapes; r++) {
297 for (int k = 0; k < p.numEqu; k++) {
298 for (int m = 0; m < p.numComp; m++) {
299 double f = 0;
300 for (int q = 0; q < p.numQuad; q++) {
301 f +=
302 vol * S[INDEX2(s, q, p.numShapes)] *
303 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
304 S[INDEX2(r, q, p.numShapes)];
305 }
306 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
307 }
308 }
309 }
310 }
311 } else {
312 for (int s = 0; s < p.numShapes; s++) {
313 for (int r = 0; r < p.numShapes; r++) {
314 double f = 0;
315 for (int q = 0; q < p.numQuad; q++)
316 f += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
317 for (int k = 0; k < p.numEqu; k++) {
318 for (int m = 0; m < p.numComp; m++) {
319 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
320 f * D_p[INDEX2(k, m, p.numEqu)];
321 }
322 }
323 }
324 }
325 }
326 }
327 ///////////////
328 // process X //
329 ///////////////
330 if (!X.isEmpty()) {
331 const double* X_p = X.getSampleDataRO(e);
332 add_EM_F = true;
333 if (expandedX) {
334 const double* X_q = &X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)];
335 for (int s = 0; s < p.numShapes; s++) {
336 for (int k = 0; k < p.numEqu; k++) {
337 double f = 0;
338 for (int q = 0; q < p.numQuad; q++) {
339 f +=
340 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
341 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
342 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
343 X_q[INDEX3(k, 1, q, p.numEqu, DIM)] +
344 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] *
345 X_q[INDEX3(k, 2, q, p.numEqu, DIM)]);
346 }
347 EM_F[INDEX2(k, s, p.numEqu)] += f;
348 }
349 }
350 } else {
351 for (int s = 0; s < p.numShapes; s++) {
352 double f0 = 0;
353 double f1 = 0;
354 double f2 = 0;
355 for (int q = 0; q < p.numQuad; q++) {
356 f0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
357 f1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
358 f2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)];
359 }
360 for (int k = 0; k < p.numEqu; k++) {
361 EM_F[INDEX2(k, s, p.numEqu)] += f0 * X_p[INDEX2(k, 0, p.numEqu)]
362 + f1 * X_p[INDEX2(k, 1, p.numEqu)] + f2 * X_p[INDEX2(k, 2, p.numEqu)];
363 }
364 }
365 }
366 }
367 ///////////////
368 // process Y //
369 ///////////////
370 if (!Y.isEmpty()) {
371 const double* Y_p = Y.getSampleDataRO(e);
372 add_EM_F = true;
373 if (expandedY) {
374 const double* Y_q = &Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)];
375 for (int s = 0; s < p.numShapes; s++) {
376 for (int k = 0; k < p.numEqu; k++) {
377 double f = 0.;
378 for (int q = 0; q < p.numQuad; q++)
379 f += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
380 EM_F[INDEX2(k, s, p.numEqu)] += f;
381 }
382 }
383 } else {
384 for (int s = 0; s < p.numShapes; s++) {
385 double f = 0;
386 for (int q = 0; q < p.numQuad; q++)
387 f += vol * S[INDEX2(s, q, p.numShapes)];
388 for (int k = 0; k < p.numEqu; k++)
389 EM_F[INDEX2(k, s, p.numEqu)] += f * Y_p[k];
390 }
391 }
392 }
393
394 // add the element matrices onto the matrix and right
395 // hand side
396 for (int q = 0; q < p.numShapes; q++)
397 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
398
399 if (add_EM_F)
400 Dudley_Util_AddScatter(p.numShapes, &row_index[0],
401 p.numEqu, &EM_F[0], F_p, p.row_DOF_UpperBound);
402 if (add_EM_S)
403 Assemble_addToSystemMatrix(mat, p.numShapes,
404 &row_index[0], p.numEqu, p.numShapes,
405 &row_index[0], p.numComp, &EM_S[0]);
406 } // end color check
407 } // end element loop
408 } // end color loop
409 } // end parallel region
410 }
411
412 } // namespace dudley
413

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/Assemble_PDE_System2_3D.cpp:5567-5588 /branches/4.0fordebian/dudley/src/Assemble_PDE_System_3D.cpp:5567-5588 /branches/complex/dudley/src/Assemble_PDE_System2_3D.cpp:5866-5937 /branches/complex/dudley/src/Assemble_PDE_System_3D.cpp:5866-5937 /branches/diaplayground/dudley/src/Assemble_PDE_System_3D.cpp:4940-5147 /branches/lapack2681/finley/src/Assemble_PDE_System2_3D.cpp:2682-2741 /branches/lapack2681/finley/src/Assemble_PDE_System_3D.cpp:2682-2741 /branches/pasowrap/dudley/src/Assemble_PDE_System2_3D.cpp:3661-3674 /branches/pasowrap/dudley/src/Assemble_PDE_System_3D.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Assemble_PDE_System2_3D.cpp:3871-3891 /branches/py3_attempt2/dudley/src/Assemble_PDE_System_3D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_System2_3D.cpp:2610-2624 /branches/restext/finley/src/Assemble_PDE_System_3D.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Assemble_PDE_System2_3D.cpp:3669-3791 /branches/ripleygmg_from_3668/dudley/src/Assemble_PDE_System_3D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_System2_3D.cpp:2569-2590 /branches/stage3.0/finley/src/Assemble_PDE_System_3D.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Assemble_PDE_System2_3D.cpp:3471-3974 /branches/symbolic_from_3470/dudley/src/Assemble_PDE_System_3D.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_PDE_System2_3D.cpp:3517-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_PDE_System_3D.cpp:3517-3974 /release/3.0/finley/src/Assemble_PDE_System2_3D.cpp:2591-2601 /release/3.0/finley/src/Assemble_PDE_System_3D.cpp:2591-2601 /release/4.0/dudley/src/Assemble_PDE_System2_3D.cpp:5380-5406 /release/4.0/dudley/src/Assemble_PDE_System_3D.cpp:5380-5406 /trunk/dudley/src/Assemble_PDE_System2_3D.cpp:4257-4344 /trunk/dudley/src/Assemble_PDE_System_3D.cpp:5898-5962,5982-6007 /trunk/ripley/test/python/dudley/src/Assemble_PDE_System2_3D.cpp:3480-3515 /trunk/ripley/test/python/dudley/src/Assemble_PDE_System_3D.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26