# Contents of /trunk/finley/src/Assemble_PDE_System2_2D.c

Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years, 2 months ago) by jfenwick
File MIME type: text/plain
File size: 19691 byte(s)
```Merging dudley and scons updates from branches

```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 by University of Queensland 5 * 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 14 15 /**************************************************************/ 16 17 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */ 18 /* the shape functions for test and solution must be identical */ 19 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 38 39 #include "Assemble.h" 40 #include "Util.h" 41 #ifdef _OPENMP 42 #include 43 #endif 44 45 46 /**************************************************************/ 47 48 void Finley_Assemble_PDE_System2_2D(Finley_Assemble_Parameters p, 49 Finley_ElementFile* elements, 50 Paso_SystemMatrix* Mat, escriptDataC* F, 51 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) { 52 53 #define DIM 2 54 index_t color; 55 dim_t e, isub; 56 __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; 57 double *EM_S, *EM_F, *Vol, *DSDX; 58 index_t *row_index; 59 register dim_t q, s,r,k,m; 60 register double rtmp, rtmp0, rtmp1, rtmp00, rtmp10, rtmp01, rtmp11; 61 bool_t add_EM_F, add_EM_S; 62 63 bool_t extendedA=isExpanded(A); 64 bool_t extendedB=isExpanded(B); 65 bool_t extendedC=isExpanded(C); 66 bool_t extendedD=isExpanded(D); 67 bool_t extendedX=isExpanded(X); 68 bool_t extendedY=isExpanded(Y); 69 double *F_p=(requireWrite(F), getSampleDataRW(F,0)); /* use comma, to get around the mixed code and declarations thing */ 70 double *S=p.row_jac->BasisFunctions->S; 71 dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp; 72 dim_t len_EM_F=p.row_numShapesTotal*p.numEqu; 73 74 #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, isub) 75 { 76 77 EM_S=THREAD_MEMALLOC(len_EM_S,double); 78 EM_F=THREAD_MEMALLOC(len_EM_F,double); 79 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t); 80 81 if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) { 82 83 for (color=elements->minColor;color<=elements->maxColor;color++) { 84 /* open loop over all elements: */ 85 #pragma omp for private(e) schedule(static) 86 for(e=0;enumElements;e++){ 87 if (elements->Color[e]==color) { 88 89 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 96 for (isub=0; isubvolume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]); 99 DSDX=&(p.row_jac->DSDX[INDEX5(0,0,0,isub,e,p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)]); 100 for (q=0;qNodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]]; 336 337 if (add_EM_F) Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound); 338 if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_numShapesTotal,row_index,p.numEqu,p.col_numShapesTotal,row_index,p.numComp,EM_S); 339 340 } /* end of isub */ 341 } /* end color check */ 342 } /* end element loop */ 343 } /* end color loop */ 344 345 THREAD_MEMFREE(EM_S); 346 THREAD_MEMFREE(EM_F); 347 THREAD_MEMFREE(row_index); 348 349 } /* end of pointer check */ 350 } /* end parallel region */ 351 } 352 /* 353 * \$Log\$ 354 */

 ViewVC Help Powered by ViewVC 1.1.26