/[escript]/trunk/finley/src/Assemble_PDE_System_2D.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_PDE_System_2D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 21804 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 /****************************************************************************
19
20 Assembles the system of numEqu PDEs into the stiffness matrix S and right
21 hand side F
22
23 -(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
24 and
25 -(X_{k,i})_i + Y_k
26
27 u has p.numComp components in a 2D domain. The shape functions for test and
28 solution must be identical and row_NS == row_NN.
29
30 Shape of the coefficients:
31
32 A = p.numEqu x 2 x p.numComp x 2
33 B = 2 x p.numEqu x p.numComp
34 C = p.numEqu x 2 x p.numComp
35 D = p.numEqu x p.numComp
36 X = p.numEqu x 2
37 Y = p.numEqu
38
39 *****************************************************************************/
40
41 #include "Assemble.h"
42 #include "Util.h"
43
44 #include <escript/index.h>
45
46 namespace finley {
47
48 template<typename Scalar>
49 void Assemble_PDE_System_2D(const AssembleParameters& p,
50 const escript::Data& A, const escript::Data& B,
51 const escript::Data& C, const escript::Data& D,
52 const escript::Data& X, const escript::Data& Y)
53 {
54 const int DIM = 2;
55 bool expandedA = A.actsExpanded();
56 bool expandedB = B.actsExpanded();
57 bool expandedC = C.actsExpanded();
58 bool expandedD = D.actsExpanded();
59 bool expandedX = X.actsExpanded();
60 bool expandedY = Y.actsExpanded();
61 const Scalar zero = static_cast<Scalar>(0);
62 Scalar* F_p = NULL;
63 if (!p.F.isEmpty()) {
64 p.F.requireWrite();
65 F_p = p.F.getSampleDataRW(0, zero);
66 }
67 const std::vector<double>& S(p.row_jac->BasisFunctions->S);
68 const size_t len_EM_S = p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
69 const size_t len_EM_F = p.row_numShapesTotal*p.numEqu;
70
71 #pragma omp parallel
72 {
73 std::vector<Scalar> EM_S(len_EM_S);
74 std::vector<Scalar> EM_F(len_EM_F);
75 IndexVector row_index(p.row_numShapesTotal);
76
77 for (index_t color = p.elements->minColor; color <= p.elements->maxColor; color++) {
78 // loop over all elements
79 #pragma omp for
80 for (index_t e = 0; e < p.elements->numElements; e++) {
81 if (p.elements->Color[e] == color) {
82 for (int isub = 0; isub < p.numSub; isub++) {
83 const double* Vol = &(p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)]);
84 const double* DSDX = &(p.row_jac->DSDX[INDEX5(0,0,0,isub,e, p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)]);
85 std::fill(EM_S.begin(), EM_S.end(), zero);
86 std::fill(EM_F.begin(), EM_F.end(), zero);
87 bool add_EM_F = false;
88 bool add_EM_S = false;
89 ///////////////
90 // process A //
91 ///////////////
92 if (!A.isEmpty()) {
93 const Scalar* A_p = A.getSampleDataRO(e, zero);
94 add_EM_S = true;
95 if (expandedA) {
96 const Scalar* A_q = &A_p[INDEX6(0,0,0,0,0,isub,p.numEqu,DIM,p.numComp,DIM,p.numQuadSub)];
97 for (int s = 0; s < p.row_numShapes; s++) {
98 for (int r = 0; r < p.col_numShapes; r++) {
99 for (int k = 0; k < p.numEqu; k++) {
100 for (int m = 0; m < p.numComp; m++) {
101 Scalar f = zero;
102 for (int q = 0; q < p.numQuadSub; q++) {
103 f += Vol[q]*(
104 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
105 + DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,0,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
106 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
107 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX5(k,1,m,1,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
108 }
109 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
110 }
111 }
112 }
113 }
114 } else { // constant A
115 for (int s = 0; s < p.row_numShapes; s++) {
116 for (int r = 0; r < p.col_numShapes; r++) {
117 Scalar f00 = zero;
118 Scalar f01 = zero;
119 Scalar f10 = zero;
120 Scalar f11 = zero;
121 for (int q = 0; q < p.numQuadSub; q++) {
122 const Scalar f0 = Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
123 const Scalar f1 = Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
124 f00 += f0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
125 f01 += f0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
126 f10 += f1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
127 f11 += f1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
128 }
129 for (int k = 0; k < p.numEqu; k++) {
130 for (int m = 0; m < p.numComp; m++) {
131 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
132 f00 * A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)]
133 + f01 * A_p[INDEX4(k,0,m,1,p.numEqu,DIM,p.numComp)]
134 + f10 * A_p[INDEX4(k,1,m,0,p.numEqu,DIM,p.numComp)]
135 + f11 * A_p[INDEX4(k,1,m,1,p.numEqu,DIM,p.numComp)];
136 }
137 }
138 }
139 }
140 }
141 }
142 ///////////////
143 // process B //
144 ///////////////
145 if (!B.isEmpty()) {
146 const Scalar* B_p = B.getSampleDataRO(e, zero);
147 add_EM_S = true;
148 if (expandedB) {
149 const Scalar* B_q = &B_p[INDEX5(0,0,0,0,isub,p.numEqu,DIM,p.numComp,p.numQuadSub)];
150 for (int s = 0; s < p.row_numShapes; s++) {
151 for (int r = 0; r < p.col_numShapes; r++) {
152 for (int k = 0; k < p.numEqu; k++) {
153 for (int m = 0; m < p.numComp; m++) {
154 Scalar f = zero;
155 for (int q = 0; q < p.numQuadSub; q++) {
156 f += Vol[q]*S[INDEX2(r,q,p.row_numShapes)]*(
157 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]
158 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX4(k,1,m,q,p.numEqu,DIM,p.numComp)]);
159 }
160 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
161 }
162 }
163 }
164 }
165 } else { // constant B
166 for (int s = 0; s < p.row_numShapes; s++) {
167 for (int r = 0; r < p.col_numShapes; r++) {
168 Scalar f0 = zero;
169 Scalar f1 = zero;
170 for (int q = 0; q < p.numQuadSub; q++) {
171 const Scalar f = Vol[q]*S[INDEX2(r,q,p.row_numShapes)];
172 f0 += f * DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
173 f1 += f * DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
174 }
175 for (int k = 0; k < p.numEqu; k++) {
176 for (int m = 0; m < p.numComp; m++) {
177 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
178 f0 * B_p[INDEX3(k,0,m,p.numEqu,DIM)]
179 + f1 * B_p[INDEX3(k,1,m,p.numEqu,DIM)];
180 }
181 }
182 }
183 }
184 }
185 }
186 ///////////////
187 // process C //
188 ///////////////
189 if (!C.isEmpty()) {
190 const Scalar* C_p = C.getSampleDataRO(e, zero);
191 add_EM_S = true;
192 if (expandedC) {
193 const Scalar* C_q = &C_p[INDEX5(0,0,0,0,isub,p.numEqu,p.numComp,DIM,p.numQuadSub)];
194 for (int s = 0; s < p.row_numShapes; s++) {
195 for (int r = 0; r < p.col_numShapes; r++) {
196 for (int k = 0; k < p.numEqu; k++) {
197 for (int m = 0; m < p.numComp; m++) {
198 Scalar f = zero;
199 for (int q = 0; q < p.numQuadSub; q++) {
200 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*(
201 C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
202 + C_q[INDEX4(k,m,1,q,p.numEqu,p.numComp,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
203 }
204 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
205 }
206 }
207 }
208 }
209 } else { // constant C
210 for (int s = 0; s < p.row_numShapes; s++) {
211 for (int r = 0; r < p.col_numShapes; r++) {
212 Scalar f0 = zero;
213 Scalar f1 = zero;
214 for (int q = 0; q < p.numQuadSub; q++) {
215 const Scalar f = Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
216 f0 += f * DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
217 f1 += f * DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
218 }
219 for (int k = 0; k < p.numEqu; k++) {
220 for (int m = 0; m < p.numComp; m++) {
221 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
222 f0 * C_p[INDEX3(k,m,0,p.numEqu,p.numComp)]
223 + f1 * C_p[INDEX3(k,m,1,p.numEqu,p.numComp)];
224 }
225 }
226 }
227 }
228 }
229 }
230 ///////////////
231 // process D //
232 ///////////////
233 if (!D.isEmpty()) {
234 const Scalar* D_p = D.getSampleDataRO(e, zero);
235 add_EM_S = true;
236 if (expandedD) {
237 const Scalar* D_q = &D_p[INDEX4(0,0,0,isub,p.numEqu,p.numComp,p.numQuadSub)];
238 for (int s = 0; s < p.row_numShapes; s++) {
239 for (int r = 0; r < p.col_numShapes; r++) {
240 for (int k = 0; k < p.numEqu; k++) {
241 for (int m = 0; m < p.numComp; m++) {
242 Scalar f = zero;
243 for (int q = 0; q < p.numQuadSub; q++) {
244 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
245 }
246 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
247 }
248 }
249 }
250 }
251 } else { // constant D
252 for (int s = 0; s < p.row_numShapes; s++) {
253 for (int r = 0; r < p.col_numShapes; r++) {
254 Scalar f = zero;
255 for (int q = 0; q < p.numQuadSub; q++) {
256 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
257 }
258 for (int k = 0; k < p.numEqu; k++) {
259 for (int m = 0; m < p.numComp; m++) {
260 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*D_p[INDEX2(k,m,p.numEqu)];
261 }
262 }
263 }
264 }
265 }
266 }
267 ///////////////
268 // process X //
269 ///////////////
270 if (!X.isEmpty()) {
271 const Scalar* X_p = X.getSampleDataRO(e, zero);
272 add_EM_F = true;
273 if (expandedX) {
274 const Scalar* X_q = &X_p[INDEX4(0,0,0,isub,p.numEqu,DIM,p.numQuadSub)];
275 for (int s = 0; s < p.row_numShapes; s++) {
276 for (int k = 0; k < p.numEqu; k++) {
277 Scalar f = zero;
278 for (int q = 0; q < p.numQuadSub; q++) {
279 f += Vol[q]*(
280 DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)]
281 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,1,q,p.numEqu,DIM)]);
282 }
283 EM_F[INDEX2(k,s,p.numEqu)] += f;
284 }
285 }
286 } else { // constant X
287 for (int s = 0; s < p.row_numShapes; s++) {
288 Scalar f0 = zero;
289 Scalar f1 = zero;
290 for (int q = 0; q < p.numQuadSub; q++) {
291 f0 += Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
292 f1 += Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
293 }
294 for (int k = 0; k < p.numEqu; k++) {
295 EM_F[INDEX2(k,s,p.numEqu)] +=
296 f0 * X_p[INDEX2(k,0,p.numEqu)]
297 + f1 * X_p[INDEX2(k,1,p.numEqu)];
298 }
299 }
300 }
301 }
302 ///////////////
303 // process Y //
304 ///////////////
305 if (!Y.isEmpty()) {
306 const Scalar* Y_p = Y.getSampleDataRO(e, zero);
307 add_EM_F = true;
308 if (expandedY) {
309 const Scalar* Y_q = &Y_p[INDEX3(0,0,isub,p.numEqu,p.numQuadSub)];
310 for (int s = 0; s < p.row_numShapes; s++) {
311 for (int k = 0; k < p.numEqu; k++) {
312 Scalar f = zero;
313 for (int q = 0; q < p.numQuadSub; q++)
314 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
315 EM_F[INDEX2(k,s,p.numEqu)]+=f;
316 }
317 }
318 } else { // constant Y
319 for (int s = 0; s < p.row_numShapes; s++) {
320 Scalar f = zero;
321 for (int q = 0; q < p.numQuadSub; q++) {
322 f += Vol[q] * S[INDEX2(s,q,p.row_numShapes)];
323 }
324 for (int k = 0; k < p.numEqu; k++)
325 EM_F[INDEX2(k,s,p.numEqu)] += f * Y_p[k];
326 }
327 }
328 }
329 // add the element matrices onto the matrix and
330 // right hand side
331 for (int q = 0; q < p.row_numShapesTotal; q++) {
332 row_index[q] = p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
333 }
334 if (add_EM_F) {
335 util::addScatter(p.row_numShapesTotal,
336 &row_index[0], p.numEqu, &EM_F[0], F_p,
337 p.row_DOF_UpperBound);
338 }
339 if (add_EM_S) {
340 Assemble_addToSystemMatrix(p.S,
341 p.row_numShapesTotal, &row_index[0],
342 p.numEqu, p.col_numShapesTotal,
343 &row_index[0], p.numComp, &EM_S[0]);
344 }
345 } // end of isub
346 } // end color check
347 } // end element loop
348 } // end color loop
349 } // end parallel region
350 }
351
352 // instantiate our two supported versions
353 template void Assemble_PDE_System_2D<escript::DataTypes::real_t>(
354 const AssembleParameters& p,
355 const escript::Data& A, const escript::Data& B,
356 const escript::Data& C, const escript::Data& D,
357 const escript::Data& X, const escript::Data& Y);
358 template void Assemble_PDE_System_2D<escript::DataTypes::cplx_t>(
359 const AssembleParameters& p,
360 const escript::Data& A, const escript::Data& B,
361 const escript::Data& C, const escript::Data& D,
362 const escript::Data& X, const escript::Data& Y);
363
364 } // namespace finley
365

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_PDE_System_2D.cpp:5567-5588 /branches/lapack2681/finley/src/Assemble_PDE_System2_2D.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_PDE_System2_2D.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_PDE_System2_2D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_System2_2D.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_PDE_System2_2D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_System2_2D.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_PDE_System2_2D.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Assemble_PDE_System_2D.cpp:5898-6118 /release/3.0/finley/src/Assemble_PDE_System2_2D.cpp:2591-2601 /release/4.0/finley/src/Assemble_PDE_System_2D.cpp:5380-5406 /trunk/finley/src/Assemble_PDE_System2_2D.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26