/[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 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 10 months ago) by caltinay
File size: 19341 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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

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