# Contents of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_Single_3D.cpp

Revision 6090 - (show annotations)
Wed Mar 23 06:35:54 2016 UTC (3 years, 1 month ago) by caltinay
File size: 19204 byte(s)
```Simplified dudley PDE routines.

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2016 by The University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 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 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 3D 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 = 3 x 3 29 B = 3 30 C = 3 31 D = scalar 32 X = 3 33 Y = scalar 34 35 *****************************************************************************/ 36 37 #include "Assemble.h" 38 #include "Util.h" 39 40 namespace dudley { 41 42 void Assemble_PDE_Single_3D(const AssembleParameters& p, 43 const escript::Data& A, const escript::Data& B, 44 const escript::Data& C, const escript::Data& D, 45 const escript::Data& X, const escript::Data& Y) 46 { 47 const int DIM = 3; 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 double* F_p = NULL; 55 if (!p.F.isEmpty()) { 56 p.F.requireWrite(); 57 F_p = p.F.getSampleDataRW(0); 58 } 59 const double* S = p.shapeFns; 60 const int len_EM_S = p.numShapes * p.numShapes; 61 const int len_EM_F = p.numShapes; 62 63 #pragma omp parallel 64 { 65 std::vector EM_S(len_EM_S); 66 std::vector EM_F(len_EM_F); 67 std::vector row_index(len_EM_F); 68 69 for (index_t color = p.elements->minColor; color <= p.elements->maxColor; color++) { 70 // loop over all elements 71 #pragma omp for 72 for (index_t e = 0; e < p.elements->numElements; e++) { 73 if (p.elements->Color[e] == color) { 74 const double vol = p.jac->absD[e] * p.jac->quadweight; 75 const double* DSDX = &p.jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]; 76 std::fill(EM_S.begin(), EM_S.end(), 0); 77 std::fill(EM_F.begin(), EM_F.end(), 0); 78 bool add_EM_F = false; 79 bool add_EM_S = false; 80 81 /////////////// 82 // process A // 83 /////////////// 84 if (!A.isEmpty()) 85 { 86 const double* A_p = A.getSampleDataRO(e); 87 add_EM_S = true; 88 if (expandedA) { 89 const double* A_q = &A_p[INDEX4(0, 0, 0, 0, DIM, DIM, p.numQuad)]; 90 for (int s = 0; s < p.numShapes; s++) { 91 for (int r = 0; r < p.numShapes; r++) { 92 double f = 0; 93 for (int q = 0; q < p.numQuad; q++) { 94 f += 95 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 96 A_q[INDEX3(0, 0, q, DIM, DIM)] * 97 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 98 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 99 A_q[INDEX3(0, 1, q, DIM, DIM)] * 100 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] + 101 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 102 A_q[INDEX3(0, 2, q, DIM, DIM)] * 103 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] + 104 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 105 A_q[INDEX3(1, 0, q, DIM, DIM)] * 106 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 107 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 108 A_q[INDEX3(1, 1, q, DIM, DIM)] * 109 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] + 110 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 111 A_q[INDEX3(1, 2, q, DIM, DIM)] * 112 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)] + 113 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * 114 A_q[INDEX3(2, 0, q, DIM, DIM)] * 115 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 116 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * 117 A_q[INDEX3(2, 1, q, DIM, DIM)] * 118 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] + 119 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * 120 A_q[INDEX3(2, 2, q, DIM, DIM)] * 121 DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]); 122 } 123 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += f; 124 } 125 } 126 } else { 127 for (int s = 0; s < p.numShapes; s++) { 128 for (int r = 0; r < p.numShapes; r++) { 129 double f00 = 0; 130 double f01 = 0; 131 double f02 = 0; 132 double f10 = 0; 133 double f11 = 0; 134 double f12 = 0; 135 double f20 = 0; 136 double f21 = 0; 137 double f22 = 0; 138 for (int q = 0; q < p.numQuad; q++) { 139 140 const double f0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 141 f00 += f0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 142 f01 += f0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 143 f02 += f0 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]; 144 145 const double f1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 146 f10 += f1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 147 f11 += f1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 148 f12 += f1 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]; 149 150 const double f2 = vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)]; 151 f20 += f2 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 152 f21 += f2 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 153 f22 += f2 * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]; 154 } 155 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += 156 f00 * A_p[INDEX2(0, 0, DIM)] + f01 * A_p[INDEX2(0, 1, DIM)] + 157 f02 * A_p[INDEX2(0, 2, DIM)] + f10 * A_p[INDEX2(1, 0, DIM)] + 158 f11 * A_p[INDEX2(1, 1, DIM)] + f12 * A_p[INDEX2(1, 2, DIM)] + 159 f20 * A_p[INDEX2(2, 0, DIM)] + f21 * A_p[INDEX2(2, 1, DIM)] + 160 f22 * A_p[INDEX2(2, 2, DIM)]; 161 } 162 } 163 } 164 } 165 /////////////// 166 // process B // 167 /////////////// 168 if (!B.isEmpty()) { 169 const double* B_p = B.getSampleDataRO(e); 170 add_EM_S = true; 171 if (expandedB) { 172 const double* B_q = &B_p[INDEX3(0, 0, 0, DIM, p.numQuad)]; 173 for (int s = 0; s < p.numShapes; s++) { 174 for (int r = 0; r < p.numShapes; r++) { 175 double f = 0; 176 for (int q = 0; q < p.numQuad; q++) { 177 f += vol * S[INDEX2(r, q, p.numShapes)] * 178 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 179 B_q[INDEX2(0, q, DIM)] + 180 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 181 B_q[INDEX2(1, q, DIM)] + 182 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * B_q[INDEX2(2, q, DIM)]); 183 } 184 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += f; 185 } 186 } 187 } else { 188 for (int s = 0; s < p.numShapes; s++) { 189 for (int r = 0; r < p.numShapes; r++) { 190 double f0 = 0; 191 double f1 = 0; 192 double f2 = 0; 193 for (int q = 0; q < p.numQuad; q++) { 194 const double f = vol * S[INDEX2(r, q, p.numShapes)]; 195 f0 += f * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 196 f1 += f * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 197 f2 += f * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)]; 198 } 199 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += 200 f0 * B_p[0] + f1 * B_p[1] + f2 * B_p[2]; 201 } 202 } 203 } 204 } 205 /////////////// 206 // process C // 207 /////////////// 208 if (!C.isEmpty()) { 209 const double* C_p = C.getSampleDataRO(e); 210 add_EM_S = true; 211 if (expandedC) { 212 const double* C_q = &C_p[INDEX3(0, 0, 0, DIM, p.numQuad)]; 213 for (int s = 0; s < p.numShapes; s++) { 214 for (int r = 0; r < p.numShapes; r++) { 215 double f = 0; 216 for (int q = 0; q < p.numQuad; q++) { 217 f += vol * S[INDEX2(s, q, p.numShapes)] * 218 (C_q[INDEX2(0, q, DIM)] * 219 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 220 C_q[INDEX2(1, q, DIM)] * 221 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] + 222 C_q[INDEX2(2, q, DIM)] * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]); 223 } 224 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += f; 225 } 226 } 227 } else { 228 for (int s = 0; s < p.numShapes; s++) { 229 for (int r = 0; r < p.numShapes; r++) { 230 double f0 = 0; 231 double f1 = 0; 232 double f2 = 0; 233 for (int q = 0; q < p.numQuad; q++) { 234 const double f = vol * S[INDEX2(s, q, p.numShapes)]; 235 f0 += f * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 236 f1 += f * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 237 f2 += f * DSDX[INDEX3(r, 2, q, p.numShapes, DIM)]; 238 } 239 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += 240 f0 * C_p[0] + f1 * C_p[1] + f2 * C_p[2]; 241 } 242 } 243 } 244 } 245 /////////////// 246 // process D // 247 /////////////// 248 if (!D.isEmpty()) { 249 const double* D_p = D.getSampleDataRO(e); 250 add_EM_S = true; 251 if (expandedD) { 252 const double* D_q = &D_p[INDEX2(0, 0, p.numQuad)]; 253 for (int s = 0; s < p.numShapes; s++) { 254 for (int r = 0; r < p.numShapes; r++) { 255 double f = 0; 256 for (int q = 0; q < p.numQuad; q++) 257 f += 258 vol * S[INDEX2(s, q, p.numShapes)] * D_q[q] * 259 S[INDEX2(r, q, p.numShapes)]; 260 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += f; 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 EM_S[INDEX4(0, 0, s, r, p.numEqu, p.numEqu, p.numShapes)] += f * D_p[0]; 270 } 271 } 272 } 273 } 274 /////////////// 275 // process X // 276 /////////////// 277 if (!X.isEmpty()) { 278 const double* X_p = X.getSampleDataRO(e); 279 add_EM_F = true; 280 if (expandedX) { 281 const double* X_q = &X_p[INDEX3(0, 0, 0, DIM, p.numQuad)]; 282 for (int s = 0; s < p.numShapes; s++) { 283 double f = 0; 284 for (int q = 0; q < p.numQuad; q++) { 285 f += 286 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 287 X_q[INDEX2(0, q, DIM)] + 288 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 289 X_q[INDEX2(1, q, DIM)] + 290 DSDX[INDEX3(s, 2, q, p.numShapes, DIM)] * X_q[INDEX2(2, q, DIM)]); 291 } 292 EM_F[INDEX2(0, s, p.numEqu)] += f; 293 } 294 } else { 295 for (int s = 0; s < p.numShapes; s++) { 296 double f0 = 0; 297 double f1 = 0; 298 double f2 = 0; 299 for (int q = 0; q < p.numQuad; q++) { 300 f0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 301 f1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 302 f2 += vol * DSDX[INDEX3(s, 2, q, p.numShapes, DIM)]; 303 } 304 EM_F[INDEX2(0, s, p.numEqu)] += f0 * X_p[0] + f1 * X_p[1] + f2 * X_p[2]; 305 } 306 } 307 } 308 /////////////// 309 // process Y // 310 /////////////// 311 if (!Y.isEmpty()) { 312 const double* Y_p = Y.getSampleDataRO(e); 313 add_EM_F = true; 314 if (expandedY) { 315 const double* Y_q = &Y_p[INDEX2(0, 0, p.numQuad)]; 316 for (int s = 0; s < p.numShapes; s++) { 317 double f = 0; 318 for (int q = 0; q < p.numQuad; q++) 319 f += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[q]; 320 EM_F[INDEX2(0, s, p.numEqu)] += f; 321 } 322 } else { 323 for (int s = 0; s < p.numShapes; s++) { 324 double f = 0; 325 for (int q = 0; q < p.numQuad; q++) 326 f += vol * S[INDEX2(s, q, p.numShapes)]; 327 EM_F[INDEX2(0, s, p.numEqu)] += f * Y_p[0]; 328 } 329 } 330 } 331 // add the element matrices onto the matrix and right 332 // hand side 333 for (int q = 0; q < p.numShapes; q++) 334 row_index[q] = p.DOF[p.elements->Nodes[INDEX2(q, e, p.NN)]]; 335 336 if (add_EM_F) 337 util::addScatter(p.numShapes, &row_index[0], p.numEqu, 338 &EM_F[0], F_p, p.DOF_UpperBound); 339 if (add_EM_S) 340 Assemble_addToSystemMatrix(p.S, row_index, p.numEqu, 341 EM_S); 342 343 } // end color check 344 } // end element loop 345 } // end color loop 346 } // end parallel region 347 } 348 349 } // namespace dudley 350

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26