/[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 3171 by jfenwick, Fri Sep 10 00:31:11 2010 UTC revision 3172 by jfenwick, Fri Sep 10 01:38:04 2010 UTC
# Line 129  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu Line 129  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu
129  /*                                                            */  /*                                                            */
130  /*  Jacobean 1D manifold in 2D and 1D elements                */  /*  Jacobean 1D manifold in 2D and 1D elements                */
131  /*                                                            */  /*                                                            */
132  void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,/*double* bogus_QuadWeights,  void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,
133                                     dim_t bogus_numShape,*/ dim_t numElements, dim_t numNodes, index_t* nodes,                                     dim_t numElements, dim_t numNodes, index_t* nodes,
                                    /*double* bogus_DSDv, dim_t bogus_numTest, double* DTDv,*/  
134                                     double* dTdX, double* volume, index_t* element_id) {                                     double* dTdX, double* volume, index_t* element_id) {
135       #define DIM 2       #define DIM 2
136       #define LOCDIM 1       #define LOCDIM 1
# Line 186  void Assemble_jacobeans_2D_M1D_E1D(doubl Line 185  void Assemble_jacobeans_2D_M1D_E1D(doubl
185    
186  /**************************************************************/  /**************************************************************/
187  /*                                                            */  /*                                                            */
 /*  Jacobean 1D manifold in 2D and 2D elements                */  
 /*                                                            */  
 void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,  
                                    dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,  
                                    double* DSDv, dim_t numTest,double* DTDv,  
                                    double* dTdX, double* volume, index_t* element_id) {  
      #define DIM 2  
      #define LOCDIM 2  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00,dXdv10,dXdv01,dXdv11,  
                        dvdX00,dvdX10,dvdX01,dvdX11, D,invD,  
                        X0_loc, X1_loc;  
        #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               dXdv00=0;  
               dXdv10=0;  
               dXdv01=0;  
               dXdv11=0;  
               for (s=0;s<numShape; s++) {  
                  X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
                  X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  
                  dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
               }  
               D  =  dXdv00*dXdv11 - dXdv01*dXdv10;  
               if (D==0.) {  
                   sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);  
                   Dudley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD=1./D;  
                  dvdX00= dXdv11*invD;  
                  dvdX10=-dXdv10*invD;  
                  dvdX01=-dXdv01*invD;  
                  dvdX11= dXdv00*invD;  
   
                  for (s=0;s<numTest; s++) {  
                    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;  
                    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;  
                  }  
               }  
           volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];  
            }  
        }  
   
      }  
      #undef DIM  
      #undef LOCDIM  
 }  
   
 /**************************************************************/  
 /*                                                            */  
188  /*  Jacobean 3D                                               */  /*  Jacobean 3D                                               */
189  /*                                                            */  /*                                                            */
190  void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
                            dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,  
                            double* DSDv, dim_t numTest, double* DTDv,  
191                             double* dTdX, double* volume, index_t* element_id) {                             double* dTdX, double* volume, index_t* element_id) {
192       #define DIM 3       #define DIM 3
193       #define LOCDIM 3       #define LOCDIM 3
194       int e,q,s;       int e,q,s;
195       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
196         // numQuad is 1 or 4
197         const dim_t numShape=4, numTest=4;
198         const double quadweight=(numQuad==1)?1./6:1./24;
199         const double DTDV[4][3]={{-1, -1, -1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
200       #pragma omp parallel       #pragma omp parallel
201       {       {
202         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
# Line 261  void Assemble_jacobeans_3D(double* coord Line 205  void Assemble_jacobeans_3D(double* coord
205    
206        #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22,D,invD,X0_loc,X1_loc,X2_loc) schedule(static)        #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22,D,invD,X0_loc,X1_loc,X2_loc) schedule(static)
207         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
            for (q=0;q<numQuad;q++) {  
208                dXdv00=0;                dXdv00=0;
209                dXdv10=0;                dXdv10=0;
210                dXdv20=0;                dXdv20=0;
# Line 275  void Assemble_jacobeans_3D(double* coord Line 218  void Assemble_jacobeans_3D(double* coord
218                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
219                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
220                   X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];                   X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
221                   dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];                   dXdv00+=X0_loc*DTDV[s][0];
222                   dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];                   dXdv10+=X1_loc*DTDV[s][0];
223                   dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];                   dXdv20+=X2_loc*DTDV[s][0];
224                   dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];                   dXdv01+=X0_loc*DTDV[s][1];
225                   dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];                   dXdv11+=X1_loc*DTDV[s][1];
226                   dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];                   dXdv21+=X2_loc*DTDV[s][1];
227                   dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];                   dXdv02+=X0_loc*DTDV[s][2];
228                   dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];                   dXdv12+=X1_loc*DTDV[s][2];
229                   dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];                   dXdv22+=X2_loc*DTDV[s][2];
230                }                }
231                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
232                if (D==0.) {                if (D==0.) {
# Line 301  void Assemble_jacobeans_3D(double* coord Line 244  void Assemble_jacobeans_3D(double* coord
244                   dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;                   dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
245                   dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;                   dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
246    
247               for (q=0;q<numQuad;q++) {
248    
249    
250                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
251                     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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10+DTDV[s][2]*dvdX20;
252                     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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11+DTDV[s][2]*dvdX21;
253                     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;                     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12+DTDV[s][2]*dvdX22;
254                   }                   }
255               volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];               volume[INDEX2(q,e,numQuad)]=ABS(D)*quadweight;
256                }                }
257             }             }
258         }         }

Legend:
Removed from v.3171  
changed lines
  Added in v.3172

  ViewVC Help
Powered by ViewVC 1.1.26