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

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_Single_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: 15985 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 and 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 2D 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 = 2 x 2
29 B = 2
30 C = 2
31 D = scalar
32 X = 2
33 Y = scalar
34
35 *****************************************************************************/
36
37 #include "Assemble.h"
38 #include "Util.h"
39
40 namespace dudley {
41
42 void Assemble_PDE_Single_2D(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 = 2;
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 // process A //
83 /////////////////
84 if (!A.isEmpty()) {
85 const double* A_p = A.getSampleDataRO(e);
86 add_EM_S = true;
87 if (expandedA) {
88 const double* A_q = &A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)];
89 for (int s = 0; s < p.numShapes; s++) {
90 for (int r = 0; r < p.numShapes; r++) {
91 double f = 0.;
92 for (int q = 0; q < p.numQuad; q++) {
93 f += vol *
94 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
95 A_q[INDEX3(0, 0, q, DIM, DIM)] *
96 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
97 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
98 A_q[INDEX3(0, 1, q, DIM, DIM)] *
99 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] +
100 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
101 A_q[INDEX3(1, 0, q, DIM, DIM)] *
102 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] +
103 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] *
104 A_q[INDEX3(1, 1, q, DIM, DIM)] *
105 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
106 }
107 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
108 }
109 }
110 } else {
111 for (int s = 0; s < p.numShapes; s++) {
112 for (int r = 0; r < p.numShapes; r++) {
113 double f00 = 0;
114 double f01 = 0;
115 double f10 = 0;
116 double f11 = 0;
117 for (int q = 0; q < p.numQuad; q++) {
118 const double f0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
119 const double f1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
120 f00 += f0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
121 f01 += f0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
122 f10 += f1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
123 f11 += f1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
124 }
125 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
126 f00 * A_p[INDEX2(0, 0, DIM)] + f01 * A_p[INDEX2(0, 1, DIM)] +
127 f10 * A_p[INDEX2(1, 0, DIM)] + f11 * A_p[INDEX2(1, 1, DIM)];
128 }
129 }
130 }
131 }
132 ///////////////
133 // process B //
134 ///////////////
135 if (!B.isEmpty()) {
136 const double* B_p = B.getSampleDataRO(e);
137 add_EM_S = true;
138 if (expandedB) {
139 const double* B_q = &B_p[INDEX3(0, 0, 0, DIM, p.numQuad)];
140 for (int s = 0; s < p.numShapes; s++) {
141 for (int r = 0; r < p.numShapes; r++) {
142 double f = 0.;
143 for (int q = 0; q < p.numQuad; q++) {
144 f +=
145 vol * S[INDEX2(r, q, p.numShapes)] *
146 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
147 B_q[INDEX2(0, q, DIM)] +
148 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * B_q[INDEX2(1, q, DIM)]);
149 }
150 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
151 }
152 }
153 } else {
154 for (int s = 0; s < p.numShapes; s++) {
155 for (int r = 0; r < p.numShapes; r++) {
156 double f0 = 0;
157 double f1 = 0;
158 for (int q = 0; q < p.numQuad; q++) {
159 const double f = vol * S[INDEX2(r, q, p.numShapes)];
160 f0 += f * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
161 f1 += f * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
162 }
163 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
164 f0 * B_p[0] + f1 * B_p[1];
165 }
166 }
167 }
168 }
169 ///////////////
170 // process C //
171 ///////////////
172 if (!C.isEmpty())
173 {
174 const double* C_p = C.getSampleDataRO(e);
175 add_EM_S = true;
176 if (expandedC) {
177 const double* C_q = &C_p[INDEX3(0, 0, 0, DIM, p.numQuad)];
178 for (int s = 0; s < p.numShapes; s++) {
179 for (int r = 0; r < p.numShapes; r++) {
180 double f = 0;
181 for (int q = 0; q < p.numQuad; q++) {
182 f += vol * S[INDEX2(s, q, p.numShapes)]*
183 (C_q[INDEX2(0, q, DIM)] *
184 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]
185 + C_q[INDEX2(1, q, DIM)] *
186 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]);
187 }
188 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
189 }
190 }
191 } else {
192 for (int s = 0; s < p.numShapes; s++) {
193 for (int r = 0; r < p.numShapes; r++) {
194 double f0 = 0;
195 double f1 = 0;
196 for (int q = 0; q < p.numQuad; q++) {
197 const double f = vol * S[INDEX2(s, q, p.numShapes)];
198 f0 += f * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)];
199 f1 += f * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)];
200 }
201 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] +=
202 f0 * C_p[0] + f1 * C_p[1];
203 }
204 }
205 }
206 }
207 ///////////////
208 // process D //
209 ///////////////
210 if (!D.isEmpty())
211 {
212 const double* D_p = D.getSampleDataRO(e);
213 add_EM_S = true;
214 if (expandedD) {
215 const double* D_q = &D_p[INDEX2(0, 0, p.numQuad)];
216 for (int s = 0; s < p.numShapes; s++) {
217 for (int r = 0; r < p.numShapes; r++) {
218 double f = 0;
219 for (int q = 0; q < p.numQuad; q++)
220 f +=
221 vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] *
222 S[INDEX2(r, q, p.numShapes)];
223 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += f;
224 }
225 }
226 } else {
227 for (int s = 0; s < p.numShapes; s++) {
228 for (int r = 0; r < p.numShapes; r++) {
229 double f = 0;
230 for (int q = 0; q < p.numQuad; q++)
231 f += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)];
232 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numComp, p.numShapes)] += f * D_p[0];
233 }
234 }
235 }
236 }
237 ///////////////
238 // process X //
239 ///////////////
240 if (!X.isEmpty()) {
241 const double* X_p = X.getSampleDataRO(e);
242 add_EM_F = true;
243 if (expandedX) {
244 const double* X_q = &X_p[INDEX3(0, 0, 0, DIM, p.numQuad)];
245 for (int s = 0; s < p.numShapes; s++) {
246 double f = 0.;
247 for (int q = 0; q < p.numQuad; q++) {
248 f += vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] *
249 X_q[INDEX2(0, q, DIM)] +
250 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * X_q[INDEX2(1, q, DIM)]);
251 }
252 EM_F[INDEX2(0, s, p.numEqu)] += f;
253 }
254 } else {
255 for (int s = 0; s < p.numShapes; s++) {
256 double f0 = 0.;
257 double f1 = 0.;
258 for (int q = 0; q < p.numQuad; q++) {
259 f0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)];
260 f1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)];
261 }
262 EM_F[INDEX2(0, s, p.numEqu)] += f0*X_p[0] + f1*X_p[1];
263 }
264 }
265 }
266 ///////////////
267 // process Y //
268 ///////////////
269 if (!Y.isEmpty()) {
270 const double* Y_p = Y.getSampleDataRO(e);
271 add_EM_F = true;
272 if (expandedY) {
273 const double* Y_q = &Y_p[INDEX2(0, 0, p.numQuad)];
274 for (int s = 0; s < p.numShapes; s++) {
275 double f = 0;
276 for (int q = 0; q < p.numQuad; q++)
277 f += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q];
278 EM_F[INDEX2(0, s, p.numEqu)] += f;
279 }
280 } else {
281 for (int s = 0; s < p.numShapes; s++) {
282 double f = 0;
283 for (int q = 0; q < p.numQuad; q++)
284 f += vol * S[INDEX2(s, q, p.numShapes)];
285 EM_F[INDEX2(0, s, p.numEqu)] += f * Y_p[0];
286 }
287 }
288 }
289 // add the element matrices onto the matrix and
290 // right hand side
291 for (int q = 0; q < p.numShapes; q++)
292 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
293 if (add_EM_F)
294 util::addScatter(p.numShapes, &row_index[0],
295 p.numEqu, &EM_F[0], F_p, p.row_DOF_UpperBound);
296 if (add_EM_S)
297 Assemble_addToSystemMatrix(mat, p.numShapes,
298 &row_index[0], p.numEqu, p.numShapes,
299 &row_index[0], p.numComp, &EM_S[0]);
300 } // end color check
301 } // end element loop
302 } // end color loop
303 } // end parallel region
304 }
305
306 } // namespace dudley
307

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26