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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 20107 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 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 2D domain. The shape functions for test and
27 solution must be identical and and row_NS == row_NN.
28
29 Shape of the coefficients:
30
31 A = p.numEqu x 2 x p.numComp x 2
32 B = 2 x p.numEqu x p.numComp
33 C = p.numEqu x 2 x p.numComp
34 D = p.numEqu x p.numComp
35 X = p.numEqu x 2
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_2D(const AssembleParameters& p, const 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 = 2;
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, 1, q, p.numShapes, DIM)] *
107 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)]
108 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
109 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
110 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)]
111 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
112 }
113 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
114 }
115 }
116 }
117 }
118 } else {
119 for (int s = 0; s < p.numShapes; s++) {
120 for (int r = 0; r < p.numShapes; r++) {
121 double f00 = 0;
122 double f01 = 0;
123 double f10 = 0;
124 double f11 = 0;
125 for (int q = 0; q < p.numQuad; q++) {
126 const double f0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
127 const double f1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
128 f00 += f0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
129 f01 += f0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
130 f10 += f1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
131 f11 += f1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
132 }
133 for (int k = 0; k < p.numEqu; k++) {
134 for (int m = 0; m < p.numComp; m++)
135 {
136 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
137 f00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)]
138 + f01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)]
139 + f10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)]
140 + f11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)];
141 }
142 }
143 }
144 }
145 }
146 }
147 ///////////////
148 // process B //
149 ///////////////
150 if (!B.isEmpty()) {
151 const double* B_p = B.getSampleDataRO(e);
152 add_EM_S = true;
153 if (expandedB) {
154 const double* B_q = &B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)];
155 for (int s = 0; s < p.numShapes; s++) {
156 for (int r = 0; r < p.numShapes; r++) {
157 for (int k = 0; k < p.numEqu; k++) {
158 for (int m = 0; m < p.numComp; m++) {
159 double f = 0;
160 for (int q = 0; q < p.numQuad; q++) {
161 f += vol * S[INDEX2(r, q, p.numShapes)] *
162 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
163 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] +
164 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
165 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]);
166 }
167 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
168 }
169 }
170 }
171 }
172 } else {
173 for (int s = 0; s < p.numShapes; s++) {
174 for (int r = 0; r < p.numShapes; r++) {
175 double f0 = 0;
176 double f1 = 0;
177 for (int q = 0; q < p.numQuad; q++) {
178 const double f = vol * S[INDEX2(r, q, p.numShapes)];
179 f0 += f * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
180 f1 += f * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
181 }
182 for (int k = 0; k < p.numEqu; k++) {
183 for (int m = 0; m < p.numComp; m++) {
184 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
185 f0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] +
186 f1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)];
187 }
188 }
189 }
190 }
191 }
192 }
193 ///////////////
194 // process C //
195 ///////////////
196 if (!C.isEmpty()) {
197 const double* C_p = C.getSampleDataRO(e);
198 add_EM_S = true;
199 if (expandedC) {
200 const double* C_q = &C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)];
201 for (int s = 0; s < p.numShapes; s++) {
202 for (int r = 0; r < p.numShapes; r++) {
203 for (int k = 0; k < p.numEqu; k++) {
204 for (int m = 0; m < p.numComp; m++) {
205 double f = 0;
206 for (int q = 0; q < p.numQuad; q++) {
207 f += vol * S[INDEX2(s, q, p.numShapes)] *
208 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] *
209 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
210 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] *
211 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
212 }
213 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
214 }
215 }
216 }
217 }
218 } else {
219 for (int s = 0; s < p.numShapes; s++) {
220 for (int r = 0; r < p.numShapes; r++) {
221 double f0 = 0;
222 double f1 = 0;
223 for (int q = 0; q < p.numQuad; q++) {
224 const double f = vol * S[INDEX2(s, q, p.numShapes)];
225 f0 += f * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
226 f1 += f * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
227 }
228 for (int k = 0; k < p.numEqu; k++) {
229 for (int m = 0; m < p.numComp; m++) {
230 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
231 f0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] +
232 f1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)];
233 }
234 }
235 }
236 }
237 }
238 }
239 ///////////////
240 // process D //
241 ///////////////
242 if (!D.isEmpty()) {
243 const double* D_p = D.getSampleDataRO(e);
244 add_EM_S = true;
245 if (expandedD) {
246 const double* D_q = &D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)];
247 for (int s = 0; s < p.numShapes; s++) {
248 for (int r = 0; r < p.numShapes; r++) {
249 for (int k = 0; k < p.numEqu; k++) {
250 for (int m = 0; m < p.numComp; m++) {
251 double f = 0;
252 for (int q = 0; q < p.numQuad; q++) {
253 f +=
254 vol * S[INDEX2(s, q, p.numShapes)] *
255 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] *
256 S[INDEX2(r, q, p.numShapes)];
257 }
258 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
259 }
260 }
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 for (int k = 0; k < p.numEqu; k++) {
270 for (int m = 0; m < p.numComp; m++) {
271 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] +=
272 f * D_p[INDEX2(k, m, p.numEqu)];
273 }
274 }
275 }
276 }
277 }
278 }
279 ///////////////
280 // process X //
281 ///////////////
282 if (!X.isEmpty()) {
283 const double* X_p = X.getSampleDataRO(e);
284 add_EM_F = true;
285 if (expandedX) {
286 const double* X_q = &X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)];
287 for (int s = 0; s < p.numShapes; s++) {
288 for (int k = 0; k < p.numEqu; k++) {
289 double f = 0;
290 for (int q = 0; q < p.numQuad; q++) {
291 f +=
292 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
293 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] +
294 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
295 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]);
296 }
297 EM_F[INDEX2(k, s, p.numEqu)] += f;
298 }
299 }
300 } else {
301 for (int s = 0; s < p.numShapes; s++) {
302 double f0 = 0;
303 double f1 = 0;
304 for (int q = 0; q < p.numQuad; q++) {
305 f0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
306 f1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
307 }
308 for (int k = 0; k < p.numEqu; k++)
309 EM_F[INDEX2(k, s, p.numEqu)] +=
310 f0 * X_p[INDEX2(k, 0, p.numEqu)] + f1 * X_p[INDEX2(k, 1, p.numEqu)];
311 }
312 }
313 }
314 ///////////////
315 // process Y //
316 ///////////////
317 if (!Y.isEmpty()) {
318 const double* Y_p = Y.getSampleDataRO(e);
319 add_EM_F = true;
320 if (expandedY) {
321 const double* Y_q = &Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)];
322 for (int s = 0; s < p.numShapes; s++) {
323 for (int k = 0; k < p.numEqu; k++) {
324 double f = 0;
325 for (int q = 0; q < p.numQuad; q++)
326 f += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)];
327 EM_F[INDEX2(k, s, p.numEqu)] += f;
328 }
329 }
330 } else {
331 for (int s = 0; s < p.numShapes; s++) {
332 double f = 0;
333 for (int q = 0; q < p.numQuad; q++)
334 f += vol * S[INDEX2(s, q, p.numShapes)];
335 for (int k = 0; k < p.numEqu; k++)
336 EM_F[INDEX2(k, s, p.numEqu)] += f * Y_p[k];
337 }
338 }
339 }
340 // add the element matrices onto the matrix and right
341 // hand side
342 for (int q = 0; q < p.numShapes; q++)
343 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
344
345 if (add_EM_F)
346 util::addScatter(p.numShapes, &row_index[0],
347 p.numEqu, &EM_F[0], F_p, p.row_DOF_UpperBound);
348 if (add_EM_S)
349 Assemble_addToSystemMatrix(mat, p.numShapes,
350 &row_index[0], p.numEqu, p.numShapes,
351 &row_index[0], p.numComp, &EM_S[0]);
352
353 } // end color check
354 } // end element loop
355 } // end color loop
356 } // end parallel region
357 }
358
359 } // namespace dudley
360

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26