/[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 3169 by jfenwick, Thu Sep 9 03:21:35 2010 UTC revision 3170 by jfenwick, Thu Sep 9 05:40:25 2010 UTC
# Line 85  void Assemble_jacobeans_1D(double* coord Line 85  void Assemble_jacobeans_1D(double* coord
85  /*                                                            */  /*                                                            */
86  /*  Jacobean 2D with area element                             */  /*  Jacobean 2D with area element                             */
87  /*                                                            */  /*                                                            */
88  void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
89                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,                             double* dTdX, double* volume, index_t* element_id)
90                             double* DSDv, dim_t numTest, double* DTDv,  {
91                             double* dTdX, double* volume, index_t* element_id) {      #define DIM 2
92       #define DIM 2      #define LOCDIM 2
93       #define LOCDIM 2      register int e,q,s;
94       register int e,q,s;      char error_msg[LenErrorMsg_MAX];
95       char error_msg[LenErrorMsg_MAX];      const double quadweight=(numQuad==1)?1./2:1./6; /* numQuad is 1 or 3 */
96       double quadweight=(numQuad==1)?1./2:1./6;  /* numQuad is 1 or 3 */      const dim_t numTest=3;      // hoping this is used in constant folding
97      /* ignoring numTest param - in this case it will always be 3 */      static const int DTDV[3][2]={{-1,-1}, {1,0}, {0,1}};    // if constant folding is not applied here will need to write
98       #pragma omp parallel                              // out in full
99       {      #pragma omp parallel
100         register double dXdv00,dXdv10,dXdv01,dXdv11,      {
101                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD,      register double dXdv00,dXdv10,dXdv01,dXdv11,
102                         X0_loc, X1_loc;                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
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++)
105                dXdv00=0;      {
               dXdv10=0;  
               dXdv01=0;  
               dXdv11=0;  
   
   
106  #define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\  #define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
107  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
108  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
# Line 116  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu Line 111  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu
111  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
112  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
113    
114          dXdv00=COMPDXDV0(0);          dXdv00=0;
115          dXdv10=COMPDXDV0(1);              dXdv10=0;
116          dXdv01=COMPDXDV1(0);              dXdv01=0;
117          dXdv11=COMPDXDV1(1);              dXdv11=0;
118                for (s=0;s<3; s++) {          dXdv00=COMPDXDV0(0);
119                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];          dXdv10=COMPDXDV0(1);
120                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];          dXdv01=COMPDXDV1(0);
121  //                 dXdv00+=X0_loc*DSDv[INDEX3(s,0,0,numShape,LOCDIM)];          dXdv11=COMPDXDV1(1);
122  //                 dXdv10+=X1_loc*DSDv[INDEX3(s,0,0,numShape,LOCDIM)];              D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
123  //                 dXdv01+=X0_loc*DSDv[INDEX3(s,1,0,numShape,LOCDIM)];              if (D==0.)
124  //                 dXdv11+=X1_loc*DSDv[INDEX3(s,1,0,numShape,LOCDIM)];          {
125                }          sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
126                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;                  Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
127                if (D==0.) {              }
128                    sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);          else
129                    Dudley_setError(ZERO_DIVISION_ERROR,error_msg);          {
130                } else {          invD=1./D;
131                   invD=1./D;          dvdX00= dXdv11*invD;
132                   dvdX00= dXdv11*invD;          dvdX10=-dXdv10*invD;
133                   dvdX10=-dXdv10*invD;          dvdX01=-dXdv01*invD;
134                   dvdX01=-dXdv01*invD;          dvdX11= dXdv00*invD;
135                   dvdX11= dXdv00*invD;          if (numQuad==1)
136           if (numQuad==1)          {
137           {              for (s=0;s<3; s++)
             for (s=0;s<3; s++) {  
             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;  
             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;  
       
             }  
             volume[INDEX2(0,e,numQuad)]=ABS(D)*quadweight;  
          }  
          else   // numQuad==3  
          {  
             for (q=0;q<3;++q)   // relying on unroll loops here  
138              {              {
139              for (s=0;s<3; s++) {  #define DTDXSET(P,Q) dTdX[INDEX4(s,P,Q,e,numTest,DIM,numQuad)] = DTDV[s][0]*dvdX0##P+DTDV[s][1]*dvdX1##P
140              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;  
141              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;              DTDXSET(0,0);
142                    DTDXSET(1,0);
143                }
144                volume[INDEX2(0,e,1)]=ABS(D)*quadweight;
145            }
146            else    // numQuad==3
147            {
148                for (q=0;q<numTest;++q) // relying on unroll loops to optimise this
149                {
150                for (s=0;s<3; s++)
151                {
152                    DTDXSET(0,q);
153                    DTDXSET(1,q);
154              }              }
155              volume[INDEX2(q,e,numQuad)]=ABS(D)*quadweight;                volume[INDEX2(q,e,3)]=ABS(D)*quadweight;    
156              }              }
157          }          }
158             }             }
# Line 164  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu Line 160  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu
160       }       }
161       #undef DIM       #undef DIM
162       #undef LOCDIM       #undef LOCDIM
163         #undef DTDXSET
164       #undef COMPDXDV0       #undef COMPDXDV0
165       #undef COMPDXDV1       #undef COMPDXDV1
166  }  }
# Line 216  void Assemble_jacobeans_2D_M1D_E1D(doubl Line 212  void Assemble_jacobeans_2D_M1D_E1D(doubl
212       #undef DIM       #undef DIM
213       #undef LOCDIM       #undef LOCDIM
214  }  }
215    
216  /**************************************************************/  /**************************************************************/
217  /*                                                            */  /*                                                            */
218  /*  Jacobean 1D manifold in 2D and 1D elements woth contact   */  /*  Jacobean 1D manifold in 2D and 1D elements woth contact   */

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

  ViewVC Help
Powered by ViewVC 1.1.26