# Contents of /branches/doubleplusgood/dudley/src/Assemble_PDE_System2_2D.cpp

Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 14770 byte(s)
```like sand though the hourglass
```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2013 by 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 since 2012 by School of Earth Sciences 13 * 14 *****************************************************************************/ 15 16 /************************************************************************************/ 17 18 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */ 19 /* the shape functions for test and solution must be identical */ 20 21 /* -(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 and -(X_{k,i})_i + Y_k */ 22 23 /* u has p.numComp components 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 = p.numEqu x 2 x p.numComp x 2 */ 29 /* B = 2 x p.numEqu x p.numComp */ 30 /* C = p.numEqu x 2 x p.numComp */ 31 /* D = p.numEqu x p.numComp */ 32 /* X = p.numEqu x 2 */ 33 /* Y = p.numEqu */ 34 35 /************************************************************************************/ 36 37 #include "Assemble.h" 38 #include "Util.h" 39 #ifdef _OPENMP 40 #include 41 #endif 42 43 /************************************************************************************/ 44 45 void Dudley_Assemble_PDE_System2_2D(Dudley_Assemble_Parameters p, Dudley_ElementFile * elements, 46 Paso_SystemMatrix * Mat, escriptDataC * F, 47 escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D, 48 escriptDataC * X, escriptDataC * Y) 49 { 50 51 #define DIM 2 52 index_t color; 53 dim_t e; 54 __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q; 55 double *EM_S, *EM_F, *DSDX; 56 index_t *row_index; 57 register dim_t q, s, r, k, m; 58 register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11; 59 bool_t add_EM_F, add_EM_S; 60 61 bool_t extendedA = isExpanded(A); 62 bool_t extendedB = isExpanded(B); 63 bool_t extendedC = isExpanded(C); 64 bool_t extendedD = isExpanded(D); 65 bool_t extendedX = isExpanded(X); 66 bool_t extendedY = isExpanded(Y); 67 double *F_p = (requireWrite(F), getSampleDataRW(F, 0)); /* use comma, to get around the mixed code and declarations thing */ 68 const double *S = p.shapeFns; 69 dim_t len_EM_S = p.numShapes * p.numShapes * p.numEqu * p.numComp; 70 dim_t len_EM_F = p.numShapes * p.numEqu; 71 72 #pragma omp parallel private(color,EM_S, EM_F, DSDX, A_p, B_p, C_p, D_p, X_p, Y_p, A_q, B_q, C_q, D_q, X_q, Y_q,row_index,q, s,r,k,m,rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11,add_EM_F, add_EM_S) 73 { 74 75 EM_S = new double[len_EM_S]; 76 EM_F = new double[len_EM_F]; 77 row_index = new index_t[p.numShapes]; 78 79 if (!Dudley_checkPtr(EM_S) && !Dudley_checkPtr(EM_F) && !Dudley_checkPtr(row_index)) 80 { 81 82 for (color = elements->minColor; color <= elements->maxColor; color++) 83 { 84 /* open loop over all elements: */ 85 #pragma omp for private(e) schedule(static) 86 for (e = 0; e < elements->numElements; e++) 87 { 88 if (elements->Color[e] == color) 89 { 90 double vol = p.row_jac->absD[e] * p.row_jac->quadweight; 91 92 A_p = getSampleDataRO(A, e); 93 B_p = getSampleDataRO(B, e); 94 C_p = getSampleDataRO(C, e); 95 D_p = getSampleDataRO(D, e); 96 X_p = getSampleDataRO(X, e); 97 Y_p = getSampleDataRO(Y, e); 98 DSDX = &(p.row_jac->DSDX[INDEX5(0, 0, 0, 0, e, p.numShapes, DIM, p.numQuad, 1)]); 99 for (q = 0; q < len_EM_S; ++q) 100 EM_S[q] = 0; 101 for (q = 0; q < len_EM_F; ++q) 102 EM_F[q] = 0; 103 add_EM_F = FALSE; 104 add_EM_S = FALSE; 105 106 /************************************************************************************/ 107 /* process A: */ 108 /************************************************************************************/ 109 A_p = getSampleDataRO(A, e); 110 if (NULL != A_p) 111 { 112 add_EM_S = TRUE; 113 if (extendedA) 114 { 115 A_q = &(A_p[INDEX6(0, 0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, DIM, p.numQuad)]); 116 for (s = 0; s < p.numShapes; s++) 117 { 118 for (r = 0; r < p.numShapes; r++) 119 { 120 for (k = 0; k < p.numEqu; k++) 121 { 122 for (m = 0; m < p.numComp; m++) 123 { 124 rtmp = 0; 125 for (q = 0; q < p.numQuad; q++) 126 { 127 rtmp += 128 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 129 A_q[INDEX5(k, 0, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] 130 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 131 DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 132 A_q[INDEX5(k, 0, m, 1, q, p.numEqu, DIM, p.numComp, DIM)] 133 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)] + 134 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 135 A_q[INDEX5(k, 1, m, 0, q, p.numEqu, DIM, p.numComp, DIM)] 136 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 137 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 138 A_q[INDEX5(k, 1, m, 1, q, p.numEqu, DIM, p.numComp, DIM)] 139 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]); 140 } 141 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp; 142 } 143 } 144 } 145 } 146 } 147 else 148 { 149 for (s = 0; s < p.numShapes; s++) 150 { 151 for (r = 0; r < p.numShapes; r++) 152 { 153 rtmp00 = 0; 154 rtmp01 = 0; 155 rtmp10 = 0; 156 rtmp11 = 0; 157 for (q = 0; q < p.numQuad; q++) 158 { 159 rtmp0 = vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 160 rtmp1 = vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 161 rtmp00 += rtmp0 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 162 rtmp01 += rtmp0 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 163 rtmp10 += rtmp1 * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 164 rtmp11 += rtmp1 * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 165 } 166 for (k = 0; k < p.numEqu; k++) 167 { 168 for (m = 0; m < p.numComp; m++) 169 { 170 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += 171 rtmp00 * A_p[INDEX4(k, 0, m, 0, p.numEqu, DIM, p.numComp)] 172 + rtmp01 * A_p[INDEX4(k, 0, m, 1, p.numEqu, DIM, p.numComp)] 173 + rtmp10 * A_p[INDEX4(k, 1, m, 0, p.numEqu, DIM, p.numComp)] 174 + rtmp11 * A_p[INDEX4(k, 1, m, 1, p.numEqu, DIM, p.numComp)]; 175 } 176 } 177 } 178 } 179 } 180 } 181 /************************************************************************************/ 182 /* process B: */ 183 /************************************************************************************/ 184 B_p = getSampleDataRO(B, e); 185 if (NULL != B_p) 186 { 187 add_EM_S = TRUE; 188 if (extendedB) 189 { 190 B_q = &(B_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, DIM, p.numComp, p.numQuad)]); 191 for (s = 0; s < p.numShapes; s++) 192 { 193 for (r = 0; r < p.numShapes; r++) 194 { 195 for (k = 0; k < p.numEqu; k++) 196 { 197 for (m = 0; m < p.numComp; m++) 198 { 199 rtmp = 0; 200 for (q = 0; q < p.numQuad; q++) 201 { 202 rtmp += vol * S[INDEX2(r, q, p.numShapes)] * 203 (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 204 B_q[INDEX4(k, 0, m, q, p.numEqu, DIM, p.numComp)] + 205 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 206 B_q[INDEX4(k, 1, m, q, p.numEqu, DIM, p.numComp)]); 207 } 208 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp; 209 } 210 } 211 } 212 } 213 } 214 else 215 { 216 for (s = 0; s < p.numShapes; s++) 217 { 218 for (r = 0; r < p.numShapes; r++) 219 { 220 rtmp0 = 0; 221 rtmp1 = 0; 222 for (q = 0; q < p.numQuad; q++) 223 { 224 rtmp = vol * S[INDEX2(r, q, p.numShapes)]; 225 rtmp0 += rtmp * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 226 rtmp1 += rtmp * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 227 } 228 for (k = 0; k < p.numEqu; k++) 229 { 230 for (m = 0; m < p.numComp; m++) 231 { 232 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += 233 rtmp0 * B_p[INDEX3(k, 0, m, p.numEqu, DIM)] + 234 rtmp1 * B_p[INDEX3(k, 1, m, p.numEqu, DIM)]; 235 } 236 } 237 } 238 } 239 } 240 } 241 /************************************************************************************/ 242 /* process C: */ 243 /************************************************************************************/ 244 C_p = getSampleDataRO(C, e); 245 if (NULL != C_p) 246 { 247 add_EM_S = TRUE; 248 if (extendedC) 249 { 250 C_q = &(C_p[INDEX5(0, 0, 0, 0, 0, p.numEqu, p.numComp, DIM, p.numQuad)]); 251 for (s = 0; s < p.numShapes; s++) 252 { 253 for (r = 0; r < p.numShapes; r++) 254 { 255 for (k = 0; k < p.numEqu; k++) 256 { 257 for (m = 0; m < p.numComp; m++) 258 { 259 rtmp = 0; 260 for (q = 0; q < p.numQuad; q++) 261 { 262 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * 263 (C_q[INDEX4(k, m, 0, q, p.numEqu, p.numComp, DIM)] * 264 DSDX[INDEX3(r, 0, q, p.numShapes, DIM)] + 265 C_q[INDEX4(k, m, 1, q, p.numEqu, p.numComp, DIM)] * 266 DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]); 267 } 268 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp; 269 } 270 } 271 } 272 } 273 } 274 else 275 { 276 for (s = 0; s < p.numShapes; s++) 277 { 278 for (r = 0; r < p.numShapes; r++) 279 { 280 rtmp0 = 0; 281 rtmp1 = 0; 282 for (q = 0; q < p.numQuad; q++) 283 { 284 rtmp = vol * S[INDEX2(s, q, p.numShapes)]; 285 rtmp0 += rtmp * DSDX[INDEX3(r, 0, q, p.numShapes, DIM)]; 286 rtmp1 += rtmp * DSDX[INDEX3(r, 1, q, p.numShapes, DIM)]; 287 } 288 for (k = 0; k < p.numEqu; k++) 289 { 290 for (m = 0; m < p.numComp; m++) 291 { 292 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += 293 rtmp0 * C_p[INDEX3(k, m, 0, p.numEqu, p.numComp)] + 294 rtmp1 * C_p[INDEX3(k, m, 1, p.numEqu, p.numComp)]; 295 } 296 } 297 } 298 } 299 } 300 } 301 /*********************************************************************************** */ 302 /* process D */ 303 /************************************************************************************/ 304 D_p = getSampleDataRO(D, e); 305 if (NULL != D_p) 306 { 307 add_EM_S = TRUE; 308 if (extendedD) 309 { 310 D_q = &(D_p[INDEX4(0, 0, 0, 0, p.numEqu, p.numComp, p.numQuad)]); 311 for (s = 0; s < p.numShapes; s++) 312 { 313 for (r = 0; r < p.numShapes; r++) 314 { 315 for (k = 0; k < p.numEqu; k++) 316 { 317 for (m = 0; m < p.numComp; m++) 318 { 319 rtmp = 0; 320 for (q = 0; q < p.numQuad; q++) 321 { 322 rtmp += 323 vol * S[INDEX2(s, q, p.numShapes)] * 324 D_q[INDEX3(k, m, q, p.numEqu, p.numComp)] * 325 S[INDEX2(r, q, p.numShapes)]; 326 } 327 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += rtmp; 328 } 329 } 330 } 331 } 332 } 333 else 334 { 335 for (s = 0; s < p.numShapes; s++) 336 { 337 for (r = 0; r < p.numShapes; r++) 338 { 339 rtmp = 0; 340 for (q = 0; q < p.numQuad; q++) 341 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(r, q, p.numShapes)]; 342 for (k = 0; k < p.numEqu; k++) 343 { 344 for (m = 0; m < p.numComp; m++) 345 { 346 EM_S[INDEX4(k, m, s, r, p.numEqu, p.numComp, p.numShapes)] += 347 rtmp * D_p[INDEX2(k, m, p.numEqu)]; 348 } 349 } 350 } 351 } 352 } 353 } 354 /************************************************************************************/ 355 /* process X: */ 356 /************************************************************************************/ 357 X_p = getSampleDataRO(X, e); 358 if (NULL != X_p) 359 { 360 add_EM_F = TRUE; 361 if (extendedX) 362 { 363 X_q = &(X_p[INDEX4(0, 0, 0, 0, p.numEqu, DIM, p.numQuad)]); 364 for (s = 0; s < p.numShapes; s++) 365 { 366 for (k = 0; k < p.numEqu; k++) 367 { 368 rtmp = 0; 369 for (q = 0; q < p.numQuad; q++) 370 { 371 rtmp += 372 vol * (DSDX[INDEX3(s, 0, q, p.numShapes, DIM)] * 373 X_q[INDEX3(k, 0, q, p.numEqu, DIM)] + 374 DSDX[INDEX3(s, 1, q, p.numShapes, DIM)] * 375 X_q[INDEX3(k, 1, q, p.numEqu, DIM)]); 376 } 377 EM_F[INDEX2(k, s, p.numEqu)] += rtmp; 378 } 379 } 380 } 381 else 382 { 383 for (s = 0; s < p.numShapes; s++) 384 { 385 rtmp0 = 0; 386 rtmp1 = 0; 387 for (q = 0; q < p.numQuad; q++) 388 { 389 rtmp0 += vol * DSDX[INDEX3(s, 0, q, p.numShapes, DIM)]; 390 rtmp1 += vol * DSDX[INDEX3(s, 1, q, p.numShapes, DIM)]; 391 } 392 for (k = 0; k < p.numEqu; k++) 393 EM_F[INDEX2(k, s, p.numEqu)] += 394 rtmp0 * X_p[INDEX2(k, 0, p.numEqu)] + rtmp1 * X_p[INDEX2(k, 1, p.numEqu)]; 395 } 396 } 397 } 398 /************************************************************************************/ 399 /* process Y: */ 400 /************************************************************************************/ 401 Y_p = getSampleDataRO(Y, e); 402 if (NULL != Y_p) 403 { 404 add_EM_F = TRUE; 405 if (extendedY) 406 { 407 Y_q = &(Y_p[INDEX3(0, 0, 0, p.numEqu, p.numQuad)]); 408 for (s = 0; s < p.numShapes; s++) 409 { 410 for (k = 0; k < p.numEqu; k++) 411 { 412 rtmp = 0; 413 for (q = 0; q < p.numQuad; q++) 414 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * Y_q[INDEX2(k, q, p.numEqu)]; 415 EM_F[INDEX2(k, s, p.numEqu)] += rtmp; 416 } 417 } 418 } 419 else 420 { 421 for (s = 0; s < p.numShapes; s++) 422 { 423 rtmp = 0; 424 for (q = 0; q < p.numQuad; q++) 425 rtmp += vol * S[INDEX2(s, q, p.numShapes)]; 426 for (k = 0; k < p.numEqu; k++) 427 EM_F[INDEX2(k, s, p.numEqu)] += rtmp * Y_p[k]; 428 } 429 } 430 } 431 /*********************************************************************************************************************/ 432 /* add the element matrices onto the matrix and right hand side */ 433 /*********************************************************************************************************************/ 434 for (q = 0; q < p.numShapes; q++) 435 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]]; 436 437 if (add_EM_F) 438 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_F, F_p, p.row_DOF_UpperBound); 439 if (add_EM_S) 440 Dudley_Assemble_addToSystemMatrix(Mat, p.numShapes, row_index, p.numEqu, 441 p.numShapes, row_index, p.numComp, EM_S); 442 443 } /* end color check */ 444 } /* end element loop */ 445 } /* end color loop */ 446 447 delete[] EM_S; 448 delete[] EM_F; 449 delete[] row_index; 450 451 } /* end of pointer check */ 452 } /* end parallel region */ 453 } 454 455 /* 456 * \$Log\$ 457 */

 ViewVC Help Powered by ViewVC 1.1.26