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

Contents of /trunk/finley/src/Assemble_PDE_System_1D.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: 17172 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
21 Assembles the system of numEqu PDEs into the stiffness matrix S and right
22 hand side F
23
24 -(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
25 and
26 -(X_{k,i})_i + Y_k
27
28 u has p.numComp components in a 1D domain. The shape functions for test and
29 solution must be identical and row_NS == row_NN.
30
31 Shape of the coefficients:
32
33 A = p.numEqu x 1 x p.numComp x 1
34 B = 1 x numEqu x p.numComp
35 C = p.numEqu x 1 x p.numComp
36 D = p.numEqu x p.numComp
37 X = p.numEqu x 1
38 Y = p.numEqu
39
40 *****************************************************************************/
41
42 #include "Assemble.h"
43 #include "Util.h"
44
45 #include <escript/index.h>
46
47 namespace finley {
48
49 void Assemble_PDE_System_1D(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 = 1;
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 double *F_p = NULL;
62 if(!p.F.isEmpty()) {
63 p.F.requireWrite();
64 F_p = p.F.getSampleDataRW(0);
65 }
66 const std::vector<double>& S(p.row_jac->BasisFunctions->S);
67 const size_t len_EM_S = p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
68 const size_t len_EM_F = p.row_numShapesTotal*p.numEqu;
69
70 #pragma omp parallel
71 {
72 for (index_t color = p.elements->minColor; color <= p.elements->maxColor; color++) {
73 // loop over all elements
74 #pragma omp for
75 for (index_t e = 0; e < p.elements->numElements; e++) {
76 if (p.elements->Color[e]==color) {
77 for (int isub=0; isub<p.numSub; isub++) {
78 const double *Vol=&(p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)]);
79 const double *DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,isub,e,p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)]);
80 std::vector<double> EM_S(len_EM_S);
81 std::vector<double> EM_F(len_EM_F);
82 bool add_EM_F=false;
83 bool add_EM_S=false;
84 ///////////////
85 // process A //
86 ///////////////
87 if (!A.isEmpty()) {
88 const double *A_p=A.getSampleDataRO(e);
89 add_EM_S=true;
90 if (expandedA) {
91 const double *A_q=&(A_p[INDEX6(0,0,0,0,0,isub,p.numEqu,DIM,p.numComp,DIM,p.numQuadSub)]);
92 for (int s=0; s<p.row_numShapes; s++) {
93 for (int r=0; r<p.col_numShapes; r++) {
94 for (int k=0; k<p.numEqu; k++) {
95 for (int m=0; m<p.numComp; m++) {
96 double f=0.;
97 for (int q=0; q<p.numQuadSub; q++) {
98 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*
99 A_q[INDEX5(k,0,m,0,q,p.numEqu,DIM,p.numComp,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
100 }
101 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
102 }
103 }
104 }
105 }
106 } else { // constant A
107 for (int s=0; s<p.row_numShapes; s++) {
108 for (int r=0; r<p.col_numShapes; r++) {
109 double f=0.;
110 for (int q=0; q<p.numQuadSub; q++)
111 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
112 for (int k=0; k<p.numEqu; k++) {
113 for (int m=0; m<p.numComp; m++) {
114 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
115 f*A_p[INDEX4(k,0,m,0,p.numEqu,DIM,p.numComp)];
116 }
117 }
118 }
119 }
120 }
121 }
122 ///////////////
123 // process B //
124 ///////////////
125 if (!B.isEmpty()) {
126 const double *B_p=B.getSampleDataRO(e);
127 add_EM_S=true;
128 if (expandedB) {
129 const double *B_q=&(B_p[INDEX5(0,0,0,0,isub,p.numEqu,DIM,p.numComp,p.numQuadSub)]);
130 for (int s=0; s<p.row_numShapes; s++) {
131 for (int r=0; r<p.col_numShapes; r++) {
132 for (int k=0; k<p.numEqu; k++) {
133 for (int m=0; m<p.numComp; m++) {
134 double f=0.;
135 for (int q=0; q<p.numQuadSub; q++) {
136 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)] *
137 B_q[INDEX4(k,0,m,q,p.numEqu,DIM,p.numComp)]*S[INDEX2(r,q,p.row_numShapes)];
138 }
139 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
140 }
141 }
142 }
143 }
144 } else { // constant B
145 for (int s=0; s<p.row_numShapes; s++) {
146 for (int r=0; r<p.col_numShapes; r++) {
147 double f=0.;
148 for (int q=0; q<p.numQuadSub; q++)
149 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*S[INDEX2(r,q,p.row_numShapes)];
150 for (int k=0; k<p.numEqu; k++) {
151 for (int m=0; m<p.numComp; m++) {
152 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
153 f*B_p[INDEX3(k,0,m,p.numEqu,DIM)];
154 }
155 }
156 }
157 }
158 }
159 }
160 ///////////////
161 // process C //
162 ///////////////
163 if (!C.isEmpty()) {
164 const double *C_p=C.getSampleDataRO(e);
165 add_EM_S=true;
166 if (expandedC) {
167 const double *C_q=&(C_p[INDEX5(0,0,0,0,isub,p.numEqu,p.numComp,DIM,p.numQuadSub)]);
168 for (int s=0; s<p.row_numShapes; s++) {
169 for (int r=0; r<p.col_numShapes; r++) {
170 for (int k=0; k<p.numEqu; k++) {
171 for (int m=0; m<p.numComp; m++) {
172 double f=0.;
173 for (int q=0; q<p.numQuadSub; q++) {
174 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)] *
175 C_q[INDEX4(k,m,0,q,p.numEqu,p.numComp,DIM)] *
176 DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
177 }
178 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
179 }
180 }
181 }
182 }
183 } else { // constant C
184 for (int s=0; s<p.row_numShapes; s++) {
185 for (int r=0; r<p.col_numShapes; r++) {
186 double f=0.;
187 for (int q=0; q<p.numQuadSub; q++)
188 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
189 for (int k=0; k<p.numEqu; k++) {
190 for (int m=0; m<p.numComp; m++) {
191 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
192 f*C_p[INDEX3(k,m,0,p.numEqu,p.numComp)];
193 }
194 }
195 }
196 }
197 }
198 }
199 ///////////////
200 // process D //
201 ///////////////
202 if (!D.isEmpty()) {
203 const double *D_p=D.getSampleDataRO(e);
204 add_EM_S=true;
205 if (expandedD) {
206 const double *D_q=&(D_p[INDEX4(0,0,0,isub,p.numEqu,p.numComp,p.numQuadSub)]);
207 for (int s=0; s<p.row_numShapes; s++) {
208 for (int r=0; r<p.col_numShapes; r++) {
209 for (int k=0; k<p.numEqu; k++) {
210 for (int m=0; m<p.numComp; m++) {
211 double f=0.;
212 for (int q=0; q<p.numQuadSub; q++) {
213 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)];
214 }
215 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
216 }
217 }
218 }
219 }
220 } else { // constant D
221 for (int s=0; s<p.row_numShapes; s++) {
222 for (int r=0; r<p.col_numShapes; r++) {
223 double f=0.;
224 for (int q=0; q<p.numQuadSub; q++)
225 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
226 for (int k=0; k<p.numEqu; k++) {
227 for (int m=0; m<p.numComp; m++) {
228 EM_S[INDEX4(k,m,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*D_p[INDEX2(k,m,p.numEqu)];
229 }
230 }
231 }
232 }
233 }
234 }
235 ///////////////
236 // process X //
237 ///////////////
238 if (!X.isEmpty()) {
239 const double *X_p=X.getSampleDataRO(e);
240 add_EM_F=true;
241 if (expandedX) {
242 const double *X_q=&(X_p[INDEX4(0,0,0,isub,p.numEqu,DIM,p.numQuadSub)]);
243 for (int s=0; s<p.row_numShapes; s++) {
244 for (int k=0; k<p.numEqu; k++) {
245 double f=0.;
246 for (int q=0; q<p.numQuadSub; q++)
247 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX3(k,0,q,p.numEqu,DIM)];
248 EM_F[INDEX2(k,s,p.numEqu)]+=f;
249 }
250 }
251 } else { // constant X
252 for (int s=0; s<p.row_numShapes; s++) {
253 double f=0.;
254 for (int q=0; q<p.numQuadSub; q++)
255 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
256 for (int k=0; k<p.numEqu; k++)
257 EM_F[INDEX2(k,s,p.numEqu)]+=f*X_p[INDEX2(k,0,p.numEqu)];
258 }
259 }
260 }
261 ///////////////
262 // process Y //
263 ///////////////
264 if (!Y.isEmpty()) {
265 const double *Y_p=Y.getSampleDataRO(e);
266 add_EM_F=true;
267 if (expandedY) {
268 const double *Y_q=&(Y_p[INDEX3(0,0,isub,p.numEqu,p.numQuadSub)]);
269 for (int s=0; s<p.row_numShapes; s++) {
270 for (int k=0; k<p.numEqu; k++) {
271 double f=0.;
272 for (int q=0; q<p.numQuadSub; q++)
273 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[INDEX2(k,q,p.numEqu)];
274 EM_F[INDEX2(k,s,p.numEqu)]+=f;
275 }
276 }
277 } else { // constant Y
278 for (int s=0; s<p.row_numShapes; s++) {
279 double f=0.;
280 for (int q=0; q<p.numQuadSub; q++)
281 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
282 for (int k=0; k<p.numEqu; k++)
283 EM_F[INDEX2(k,s,p.numEqu)]+=f*Y_p[k];
284 }
285 }
286 }
287 // add the element matrices onto the matrix and
288 // right hand side
289 std::vector<index_t> row_index(p.row_numShapesTotal);
290 for (int q=0; q<p.row_numShapesTotal; q++)
291 row_index[q]=p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
292
293 if (add_EM_F)
294 util::addScatter(p.row_numShapesTotal,
295 &row_index[0], p.numEqu, &EM_F[0], F_p,
296 p.row_DOF_UpperBound);
297 if (add_EM_S)
298 Assemble_addToSystemMatrix(p.S,
299 p.row_numShapesTotal, &row_index[0],
300 p.numEqu, p.col_numShapesTotal,
301 &row_index[0], p.numComp, &EM_S[0]);
302 } // end of isub
303 } // end color check
304 } // end element loop
305 } // end color loop
306 } // end parallel region
307 }
308
309 } // namespace finley
310

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26