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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26