/[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 6090 - (show annotations)
Wed Mar 23 06:35:54 2016 UTC (3 years ago) by caltinay
File size: 24271 byte(s)
Simplified dudley PDE routines.

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

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