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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4442 - (show annotations)
Fri Jun 7 07:04:44 2013 UTC (6 years, 1 month ago) by caltinay
File size: 16129 byte(s)
cosmetic changes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by 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 since 2012 by School of Earth Sciences
13 *
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 finley {
41
42 void Assemble_PDE_Single_2D(const AssembleParameters& p,
43 escript::Data& A, escript::Data& B,
44 escript::Data& C, escript::Data& D,
45 escript::Data& X, escript::Data& Y)
46 {
47 const int DIM = 2;
48 bool expandedA=A.actsExpanded();
49 bool expandedB=B.actsExpanded();
50 bool expandedC=C.actsExpanded();
51 bool expandedD=D.actsExpanded();
52 bool expandedX=X.actsExpanded();
53 bool expandedY=Y.actsExpanded();
54 p.F.requireWrite();
55 double *F_p=p.F.getSampleDataRW(0);
56 const double *S=p.row_jac->BasisFunctions->S;
57 const int len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal;
58 const int len_EM_F=p.row_numShapesTotal;
59
60 #pragma omp parallel
61 {
62 for (int color=p.elements->minColor;color<=p.elements->maxColor;color++) {
63 // loop over all elements:
64 #pragma omp for
65 for (int e=0; e<p.elements->numElements; e++) {
66 if (p.elements->Color[e]==color) {
67 for (int isub=0; isub<p.numSub; isub++) {
68 const double *Vol=&(p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)]);
69 const double *DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,isub,e, p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)]);
70 std::vector<double> EM_S(len_EM_S);
71 std::vector<double> EM_F(len_EM_F);
72 bool add_EM_F=false;
73 bool add_EM_S=false;
74 ///////////////
75 // process A //
76 ///////////////
77 if (!A.isEmpty()) {
78 const double *A_p=A.getSampleDataRO(e);
79 add_EM_S=true;
80 if (expandedA) {
81 const double *A_q=&(A_p[INDEX4(0,0,0,isub,DIM,DIM,p.numQuadSub)]);
82 for (int s=0; s<p.row_numShapes; s++) {
83 for (int r=0; r<p.col_numShapes; r++) {
84 double f=0.;
85 for (int q=0; q<p.numQuadSub; q++) {
86 f+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
87 + DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]
88 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(1,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
89 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(1,1,q,DIM,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
90 }
91 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
92 }
93 }
94 } else { // constant A
95 for (int s=0; s<p.row_numShapes; s++) {
96 for (int r=0; r<p.col_numShapes; r++) {
97 double f00=0.;
98 double f01=0.;
99 double f10=0.;
100 double f11=0.;
101 for (int q=0; q<p.numQuadSub; q++) {
102 const double f0=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
103 const double f1=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
104 f00+=f0*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
105 f01+=f0*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
106 f10+=f1*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
107 f11+=f1*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
108 }
109 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)] +=
110 f00*A_p[INDEX2(0,0,DIM)]
111 + f01*A_p[INDEX2(0,1,DIM)]
112 + f10*A_p[INDEX2(1,0,DIM)]
113 + f11*A_p[INDEX2(1,1,DIM)];
114 }
115 }
116 }
117 }
118 ///////////////
119 // process B //
120 ///////////////
121 if (!B.isEmpty()) {
122 const double *B_p=B.getSampleDataRO(e);
123 add_EM_S=true;
124 if (expandedB) {
125 const double *B_q=&(B_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
126 for (int s=0; s<p.row_numShapes; s++) {
127 for (int r=0; r<p.col_numShapes; r++) {
128 double f=0.;
129 for (int q=0; q<p.numQuadSub; q++) {
130 f+=Vol[q]*S[INDEX2(r,q,p.row_numShapes)]*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(0,q,DIM)]
131 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(1,q,DIM)]);
132 }
133 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
134 }
135 }
136 } else { // constant B
137 for (int s=0; s<p.row_numShapes; s++) {
138 for (int r=0; r<p.col_numShapes; r++) {
139 double f0=0.;
140 double f1=0.;
141 for (int q=0; q<p.numQuadSub; q++) {
142 const double f=Vol[q]*S[INDEX2(r,q,p.row_numShapes)];
143 f0+=f*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
144 f1+=f*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
145 }
146 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f0*B_p[0]+f1*B_p[1];
147 }
148 }
149 }
150 }
151 ///////////////
152 // process C //
153 ///////////////
154 if (!C.isEmpty()) {
155 const double *C_p=C.getSampleDataRO(e);
156 add_EM_S=true;
157 if (expandedC) {
158 const double *C_q=&(C_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
159 for (int s=0; s<p.row_numShapes; s++) {
160 for (int r=0; r<p.col_numShapes; r++) {
161 double f=0.;
162 for (int q=0; q<p.numQuadSub; q++) {
163 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*(C_q[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)]
164 + C_q[INDEX2(1,q,DIM)]*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)]);
165 }
166 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
167 }
168 }
169 } else { // constant C
170 for (int s=0; s<p.row_numShapes; s++) {
171 for (int r=0; r<p.col_numShapes; r++) {
172 double f0=0.;
173 double f1=0.;
174 for (int q=0; q<p.numQuadSub; q++) {
175 const double f=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
176 f0+=f*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
177 f1+=f*DSDX[INDEX3(r,1,q,p.row_numShapesTotal,DIM)];
178 }
179 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f0*C_p[0]+f1*C_p[1];
180 }
181 }
182 }
183 }
184 ///////////////
185 // process D //
186 ///////////////
187 if (!D.isEmpty()) {
188 const double *D_p=D.getSampleDataRO(e);
189 add_EM_S=true;
190 if (expandedD) {
191 const double *D_q=&(D_p[INDEX2(0,isub,p.numQuadSub)]);
192 for (int s=0; s<p.row_numShapes; s++) {
193 for (int r=0; r<p.col_numShapes; r++) {
194 double f=0.;
195 for (int q=0; q<p.numQuadSub; q++) {
196 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[q]*S[INDEX2(r,q,p.row_numShapes)];
197 }
198 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
199 }
200 }
201 } else { // constant D
202 for (int s=0; s<p.row_numShapes; s++) {
203 for (int r=0; r<p.col_numShapes; r++) {
204 double f=0.;
205 for (int q=0; q<p.numQuadSub; q++) {
206 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
207 }
208 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*D_p[0];
209 }
210 }
211 }
212 }
213 ///////////////
214 // process X //
215 ///////////////
216 if (!X.isEmpty()) {
217 const double *X_p=X.getSampleDataRO(e);
218 add_EM_F=true;
219 if (expandedX) {
220 const double *X_q=&(X_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
221 for (int s=0; s<p.row_numShapes; s++) {
222 double f=0.;
223 for (int q=0; q<p.numQuadSub; q++) {
224 f+=Vol[q]*(DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(0,q,DIM)]
225 + DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(1,q,DIM)]);
226 }
227 EM_F[INDEX2(0,s,p.numEqu)]+=f;
228 }
229 } else { // constant X
230 for (int s=0; s<p.row_numShapes; s++) {
231 double f0=0.;
232 double f1=0.;
233 for (int q=0; q<p.numQuadSub; q++) {
234 f0+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
235 f1+=Vol[q]*DSDX[INDEX3(s,1,q,p.row_numShapesTotal,DIM)];
236 }
237 EM_F[INDEX2(0,s,p.numEqu)]+=f0*X_p[0]+f1*X_p[1];
238 }
239 }
240 }
241 ///////////////
242 // process Y //
243 ///////////////
244 if (!Y.isEmpty()) {
245 const double *Y_p=Y.getSampleDataRO(e);
246 add_EM_F=true;
247 if (expandedY) {
248 const double *Y_q=&(Y_p[INDEX2(0,isub,p.numQuadSub)]);
249 for (int s=0; s<p.row_numShapes; s++) {
250 double f=0.;
251 for (int q=0; q<p.numQuadSub; q++) {
252 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[q];
253 }
254 EM_F[INDEX2(0,s,p.numEqu)]+=f;
255 }
256 } else { // constant Y
257 for (int s=0; s<p.row_numShapes; s++) {
258 double f=0.;
259 for (int q=0; q<p.numQuadSub; q++) {
260 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
261 }
262 EM_F[INDEX2(0,s,p.numEqu)]+=f*Y_p[0];
263 }
264 }
265 }
266 // add the element matrices onto the matrix and
267 // right hand side
268 std::vector<int> row_index(p.row_numShapesTotal);
269 for (int q=0; q<p.row_numShapesTotal; q++)
270 row_index[q]=p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
271
272 if (add_EM_F)
273 util::addScatter(p.row_numShapesTotal,
274 &row_index[0], p.numEqu, &EM_F[0], F_p,
275 p.row_DOF_UpperBound);
276 if (add_EM_S)
277 Assemble_addToSystemMatrix(p.S,
278 p.row_numShapesTotal, &row_index[0],
279 p.numEqu, p.col_numShapesTotal,
280 &row_index[0], p.numComp, &EM_S[0]);
281
282 } // end of isub
283 } // end color check
284 } // end element loop
285 } // end color loop
286 } // end parallel region
287 }
288
289 } // namespace finley
290

Properties

Name Value
svn:mergeinfo /branches/lapack2681/finley/src/Assemble_PDE_Single2_2D.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_PDE_Single2_2D.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_PDE_Single2_2D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_Single2_2D.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_PDE_Single2_2D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_Single2_2D.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_PDE_Single2_2D.cpp:3471-3974 /release/3.0/finley/src/Assemble_PDE_Single2_2D.cpp:2591-2601 /trunk/finley/src/Assemble_PDE_Single2_2D.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26