/[escript]/branches/domexper/dudley/src/Assemble_jacobeans.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3086 by jfenwick, Thu Aug 5 05:07:58 2010 UTC revision 3169 by jfenwick, Thu Sep 9 03:21:35 2010 UTC
# Line 93  void Assemble_jacobeans_2D(double* coord Line 93  void Assemble_jacobeans_2D(double* coord
93       #define LOCDIM 2       #define LOCDIM 2
94       register int e,q,s;       register int e,q,s;
95       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
96         double quadweight=(numQuad==1)?1./2:1./6;  /* numQuad is 1 or 3 */
97        /* ignoring numTest param - in this case it will always be 3 */
98       #pragma omp parallel       #pragma omp parallel
99       {       {
100         register double dXdv00,dXdv10,dXdv01,dXdv11,         register double dXdv00,dXdv10,dXdv01,dXdv11,
# Line 100  void Assemble_jacobeans_2D(double* coord Line 102  void Assemble_jacobeans_2D(double* coord
102                         X0_loc, X1_loc;                         X0_loc, X1_loc;
103         #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)         #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
104         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
            for (q=0;q<numQuad;q++) {  
105                dXdv00=0;                dXdv00=0;
106                dXdv10=0;                dXdv10=0;
107                dXdv01=0;                dXdv01=0;
108                dXdv11=0;                dXdv11=0;
109                for (s=0;s<numShape; s++) {  
110    
111    #define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
112    coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
113    coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
114    
115    #define COMPDXDV1(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
116    coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
117    coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
118    
119            dXdv00=COMPDXDV0(0);
120            dXdv10=COMPDXDV0(1);
121            dXdv01=COMPDXDV1(0);
122            dXdv11=COMPDXDV1(1);
123                  for (s=0;s<3; s++) {
124                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
125                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
126                   dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  //                 dXdv00+=X0_loc*DSDv[INDEX3(s,0,0,numShape,LOCDIM)];
127                   dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  //                 dXdv10+=X1_loc*DSDv[INDEX3(s,0,0,numShape,LOCDIM)];
128                   dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  //                 dXdv01+=X0_loc*DSDv[INDEX3(s,1,0,numShape,LOCDIM)];
129                   dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  //                 dXdv11+=X1_loc*DSDv[INDEX3(s,1,0,numShape,LOCDIM)];
130                }                }
131                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
132                if (D==0.) {                if (D==0.) {
# Line 123  void Assemble_jacobeans_2D(double* coord Line 138  void Assemble_jacobeans_2D(double* coord
138                   dvdX10=-dXdv10*invD;                   dvdX10=-dXdv10*invD;
139                   dvdX01=-dXdv01*invD;                   dvdX01=-dXdv01*invD;
140                   dvdX11= dXdv00*invD;                   dvdX11= dXdv00*invD;
141             if (numQuad==1)
142                   for (s=0;s<numTest; s++) {           {
143                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;              for (s=0;s<3; s++) {
144                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;              dTdX[INDEX4(s,0,0,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,0,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,0,numTest,LOCDIM)]*dvdX10;
145                   }              dTdX[INDEX4(s,1,0,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,0,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,0,numTest,LOCDIM)]*dvdX11;
146        
147                }              }
148                volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];              volume[INDEX2(0,e,numQuad)]=ABS(D)*quadweight;
149             }
150             else   // numQuad==3
151             {
152                for (q=0;q<3;++q)   // relying on unroll loops here
153                {
154                for (s=0;s<3; s++) {
155                dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;
156                dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;
157        
158                }
159                volume[INDEX2(q,e,numQuad)]=ABS(D)*quadweight;  
160                }
161            }
162             }             }
163         }      }
   
164       }       }
165       #undef DIM       #undef DIM
166       #undef LOCDIM       #undef LOCDIM
167    
168         #undef COMPDXDV0
169         #undef COMPDXDV1
170  }  }
171  /**************************************************************/  /**************************************************************/
172  /*                                                            */  /*                                                            */

Legend:
Removed from v.3086  
changed lines
  Added in v.3169

  ViewVC Help
Powered by ViewVC 1.1.26