Annotation of /branches/trilinos_from_5897/dudley/src/Assemble_PDE_System2_2D.cpp

Revision 3187 - (hide annotations)
Thu Sep 16 02:57:17 2010 UTC (8 years, 10 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Assemble_PDE_System2_2D.c
File MIME type: text/plain
File size: 14971 byte(s)
```Enforcing indent -linux -i4 -bl -bli0 -l120

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

 ViewVC Help Powered by ViewVC 1.1.26