/[escript]/temp/finley/src/Assemble_jacobeans.c
ViewVC logotype

Diff of /temp/finley/src/Assemble_jacobeans.c

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

revision 776 by gross, Wed Jul 12 00:07:31 2006 UTC revision 777 by gross, Wed Jul 12 08:54:45 2006 UTC
# Line 24  Line 24 
24  #include <omp.h>  #include <omp.h>
25  #endif  #endif
26    
27    
28    
29  /**************************************************************/  /**************************************************************/
30  /*                                                            */  /*                                                            */
31  /*  Jacobean 1D                                               */  /*  Jacobean 1D                                               */
32  /*                                                            */  /*                                                            */
33  void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
34                             dim_t numShape, dim_t numElements,index_t* nodes,                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
35                             double* DSDv, dim_t numTest,double* DTDv,                             double* DSDv, dim_t numTest,double* DTDv,
36                             double* dTdX, double* volume, index_t* element_id) {                             double* dTdX, double* volume, index_t* element_id) {
37       #define DIM 1       #define DIM 1
# Line 42  void Assemble_jacobeans_1D(double* coord Line 44  void Assemble_jacobeans_1D(double* coord
44         double X0[numShape];         double X0[numShape];
45         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
46         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
47             for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];             for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
48             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
49                D=0;                D=0;
50                for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
# Line 51  void Assemble_jacobeans_1D(double* coord Line 53  void Assemble_jacobeans_1D(double* coord
53                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
54                } else {                } else {
55                   invD=1./D;                   invD=1./D;
56                   for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;                   for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
          volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];  
57                }                }
58              volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
59             }             }
60         }         }
61    
# Line 64  void Assemble_jacobeans_1D(double* coord Line 66  void Assemble_jacobeans_1D(double* coord
66  }  }
67  /**************************************************************/  /**************************************************************/
68  /*                                                            */  /*                                                            */
69  /*  Jacobean 2D                                               */  /*  Jacobean 2D with area element                             */
70  /*                                                            */  /*                                                            */
71  void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
72                             dim_t numShape, dim_t numElements,index_t* nodes,                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
73                             double* DSDv, dim_t numTest,double* DTDv,                             double* DSDv, dim_t numTest,double* DTDv,
74                             double* dTdX, double* volume, index_t* element_id) {                             double* dTdX, double* volume, index_t* element_id) {
75       #define DIM 2       #define DIM 2
# Line 82  void Assemble_jacobeans_2D(double* coord Line 84  void Assemble_jacobeans_2D(double* coord
84         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
85         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
86             for (s=0;s<numShape; s++) {             for (s=0;s<numShape; s++) {
87               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
88               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
89             }             }
90             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
91                dXdv00=0;                dXdv00=0;
# Line 108  void Assemble_jacobeans_2D(double* coord Line 110  void Assemble_jacobeans_2D(double* coord
110                   dvdX11= dXdv00*invD;                   dvdX11= dXdv00*invD;
111    
112                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
113                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
114                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
115                   }                   }
          volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];  
116                }                }
117                  volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
118             }             }
119         }         }
120    
# Line 122  void Assemble_jacobeans_2D(double* coord Line 124  void Assemble_jacobeans_2D(double* coord
124  }  }
125  /**************************************************************/  /**************************************************************/
126  /*                                                            */  /*                                                            */
127  /*  Jacobean 3D                                               */  /*  Jacobean 1D manifold in 2D and 1D elements                */
128  /*                                                            */  /*                                                            */
129  void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,double* QuadWeights,
130                             dim_t numShape, dim_t numElements,index_t* nodes,                                     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
131                             double* DSDv, dim_t numTest,double* DTDv,                                     double* DSDv, dim_t numTest,double* DTDv,
132                             double* dTdX, double* volume, index_t* element_id) {                                     double* dTdX, double* volume, index_t* element_id) {
133       #define DIM 3       #define DIM 2
134       #define LOCALDIM 3       #define LOCALDIM 1
135       register int e,q,s;       register int e,q,s;
136       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
137       #pragma omp parallel       #pragma omp parallel
138       {       {
139         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,         register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
140                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;         double X0[numShape], X1[numShape];
        double X0[numShape], X1[numShape], X2[numShape];  
141         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
142         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
143             for (s=0;s<numShape; s++) {             for (s=0;s<numShape; s++) {
144               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
145               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
146               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];             }
147               for (q=0;q<numQuad;q++) {
148                  dXdv00=0;
149                  dXdv10=0;
150                  for (s=0;s<numShape; s++) {
151                     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
152                     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
153                  }
154                  D=dXdv00*dXdv00+dXdv10*dXdv10;
155                  if (D==0.) {
156                      sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
157                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
158                  } else {
159                     invD=1./D;
160                     dvdX00=dXdv00*invD;
161                     dvdX01=dXdv10*invD;
162                     for (s=0;s<numTest; s++) {
163                       dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00;
164                       dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01;
165                     }
166                 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
167                  }
168               }
169           }
170    
171         }
172         #undef DIM
173         #undef LOCALDIM
174    }
175    /**************************************************************/
176    /*                                                            */
177    /*  Jacobean 1D manifold in 2D and 2D elements                */
178    /*                                                            */
179    void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
180                                       dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
181                                       double* DSDv, dim_t numTest,double* DTDv,
182                                       double* dTdX, double* volume, index_t* element_id) {
183         #define DIM 2
184         #define LOCALDIM 2
185         register int e,q,s;
186         char error_msg[LenErrorMsg_MAX];
187         #pragma omp parallel
188         {
189           register double dXdv00,dXdv10,dXdv01,dXdv11,
190                           dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
191           double X0[numShape], X1[numShape];
192           #pragma omp for private(e,q,s) schedule(static)
193           for(e=0;e<numElements;e++){
194               for (s=0;s<numShape; s++) {
195                 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
196                 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
197             }             }
198             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
199                dXdv00=0;                dXdv00=0;
200                dXdv10=0;                dXdv10=0;
               dXdv20=0;  
201                dXdv01=0;                dXdv01=0;
202                dXdv11=0;                dXdv11=0;
               dXdv21=0;  
               dXdv02=0;  
               dXdv12=0;  
               dXdv22=0;  
203                for (s=0;s<numShape; s++) {                for (s=0;s<numShape; s++) {
204                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
205                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
                  dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];  
206                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
207                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
                  dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];  
                  dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];  
                  dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];  
                  dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];  
208                }                }
209                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
210                if (D==0.) {                if (D==0.) {
211                    sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);                    sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
212                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
213                } else {                } else {
214                   invD=1./D;                   invD=1./D;
215                   dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;                   dvdX00= dXdv11*invD;
216                   dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;                   dvdX10=-dXdv10*invD;
217                   dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;                   dvdX01=-dXdv01*invD;
218                   dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;                   dvdX11= dXdv00*invD;
                  dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;  
                  dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;  
                  dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;  
                  dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;  
                  dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;  
219    
220                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
221                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX20;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10;
222                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX21;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11;
                    dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX22;  
223                   }                   }
              volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];  
224                }                }
225              volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
226             }             }
227         }         }
228    
# Line 197  void Assemble_jacobeans_3D(double* coord Line 232  void Assemble_jacobeans_3D(double* coord
232  }  }
233  /**************************************************************/  /**************************************************************/
234  /*                                                            */  /*                                                            */
235  /*  Jacobean 2D manifold in 3D                                */  /*  Jacobean 3D                                               */
236  /*                                                            */  /*                                                            */
237  void Assemble_jacobeans_3D_M2D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
238                                 dim_t numShape, dim_t numElements,index_t* nodes,                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
239                                 double* DSDv, dim_t numTest,double* DTDv,                             double* DSDv, dim_t numTest,double* DTDv,
240                                 double* dTdX, double* volume, index_t* element_id) {                             double* dTdX, double* volume, index_t* element_id) {
241       #define DIM 3       #define DIM 3
242       #define LOCALDIM 2       #define LOCALDIM 3
243       register int e,q,s;       register int e,q,s;
244       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
245       #pragma omp parallel       #pragma omp parallel
246       {       {
247         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
248                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
249         double X0[numShape], X1[numShape], X2[numShape];         double X0[numShape], X1[numShape], X2[numShape];
250         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
251         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
252             for (s=0;s<numShape; s++) {             for (s=0;s<numShape; s++) {
253               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
254               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
255               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
256             }             }
257             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
258                dXdv00=0;                dXdv00=0;
# Line 226  void Assemble_jacobeans_3D_M2D(double* c Line 261  void Assemble_jacobeans_3D_M2D(double* c
261                dXdv01=0;                dXdv01=0;
262                dXdv11=0;                dXdv11=0;
263                dXdv21=0;                dXdv21=0;
264                  dXdv02=0;
265                  dXdv12=0;
266                  dXdv22=0;
267                for (s=0;s<numShape; s++) {                for (s=0;s<numShape; s++) {
268                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
269                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
# Line 233  void Assemble_jacobeans_3D_M2D(double* c Line 271  void Assemble_jacobeans_3D_M2D(double* c
271                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
272                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
273                   dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];                   dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
274                     dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
275                     dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
276                     dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
277                }                }
278                m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
               m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;  
               m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;  
               D=m00*m11-m01*m01;  
279                if (D==0.) {                if (D==0.) {
280                    sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);                    sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
281                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
282                } else {                } else {
283                   invD=1./D;                   invD=1./D;
284                   dvdX00=( m00*dXdv00-m01*dXdv01)*invD;                   dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
285                   dvdX01=( m00*dXdv10-m01*dXdv11)*invD;                   dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
286                   dvdX02=( m00*dXdv20-m01*dXdv21)*invD;                   dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
287                   dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;                   dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
288                   dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;                   dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
289                   dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;                   dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
290                     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
291                     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
292                     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
293    
294                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
295                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX20;
296                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX21;
297                     dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12;                     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX22;
298                   }                   }
299               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];               volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
300                }                }
301             }             }
302         }         }
# Line 265  void Assemble_jacobeans_3D_M2D(double* c Line 307  void Assemble_jacobeans_3D_M2D(double* c
307  }  }
308  /**************************************************************/  /**************************************************************/
309  /*                                                            */  /*                                                            */
310  /*  Jacobean 1D manifold in 3D                                */  /*  Jacobean 2D manifold in 3D with 3D elements               */
311  /*                                                            */  /*                                                            */
312  void Assemble_jacobeans_3D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
313                                 dim_t numShape, dim_t numElements,index_t* nodes,                                     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
314                                 double* DSDv, dim_t numTest,double* DTDv,                                     double* DSDv, dim_t numTest,double* DTDv,
315                                 double* dTdX, double* volume, index_t* element_id) {                                     double* dTdX, double* volume, index_t* element_id) {
316       #define DIM 3       #define DIM 3
317       #define LOCALDIM 1       #define LOCALDIM 3
318       register int e,q,s;       register int e,q,s;
319       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
320       #pragma omp parallel       #pragma omp parallel
321       {       {
322         register double dXdv00,dXdv10,dXdv20,dvdX00,dvdX01,dvdX02,D,invD;         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
323         double X0[numShape], X1[numShape], X2[numShape];                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
324           double X0[numShape], X1[numShape], X2[numShape], m0, m1, m2;
325         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
326         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
327             for (s=0;s<numShape; s++) {             for (s=0;s<numShape; s++) {
328               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
329               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
330               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
331             }             }
332             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
333                dXdv00=0;                dXdv00=0;
334                dXdv10=0;                dXdv10=0;
335                dXdv20=0;                dXdv20=0;
336                  dXdv01=0;
337                  dXdv11=0;
338                  dXdv21=0;
339                  dXdv02=0;
340                  dXdv12=0;
341                  dXdv22=0;
342                for (s=0;s<numShape; s++) {                for (s=0;s<numShape; s++) {
343                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
344                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
345                   dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
346                     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
347                     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
348                     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
349                     dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
350                     dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
351                     dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
352                }                }
353                D=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
354                if (D==0.) {                if (D==0.) {
355                    sprintf(error_msg,"Assemble_jacobeans_3D_M1D: element %d (id %d) has length zero.",e,element_id[e]);                    sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
356                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
357                } else {                } else {
358                   invD=1./D;                   invD=1./D;
359                   dvdX00=dXdv00*invD;                   dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
360                   dvdX01=dXdv10*invD;                   dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
361                   dvdX02=dXdv20*invD;                   dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
362                     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
363                     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
364                     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
365                     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
366                     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
367                     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
368    
369                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
370                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX20;
371                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX21;
372                     dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02;                     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX22;
373                   }                   }
              volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];  
374                }                }
375                  m0=dXdv10*dXdv21-dXdv20*dXdv11;
376                  m1=dXdv20*dXdv01-dXdv00*dXdv21;
377                  m2=dXdv00*dXdv11-dXdv10*dXdv01;
378              volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
379             }             }
380         }         }
381    
# Line 320  void Assemble_jacobeans_3D_M1D(double* c Line 385  void Assemble_jacobeans_3D_M1D(double* c
385  }  }
386  /**************************************************************/  /**************************************************************/
387  /*                                                            */  /*                                                            */
388  /*  Jacobean 1D manifold in 2D                                */  /*  Jacobean 2D manifold in 3D with 2D elements               */
389  /*                                                            */  /*                                                            */
390  void Assemble_jacobeans_2D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
391                                 dim_t numShape, dim_t numElements,index_t* nodes,                                     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
392                                 double* DSDv, dim_t numTest,double* DTDv,                                     double* DSDv, dim_t numTest,double* DTDv,
393                                 double* dTdX, double* volume, index_t* element_id) {                                     double* dTdX, double* volume, index_t* element_id) {
394       #define DIM 2       #define DIM 3
395       #define LOCALDIM 1       #define LOCALDIM 2
396       register int e,q,s;       register int e,q,s;
397       char error_msg[LenErrorMsg_MAX];       char error_msg[LenErrorMsg_MAX];
398       #pragma omp parallel       #pragma omp parallel
399       {       {
400         register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
401         double X0[numShape], X1[numShape];                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
402           double X0[numShape], X1[numShape], X2[numShape];
403         #pragma omp for private(e,q,s) schedule(static)         #pragma omp for private(e,q,s) schedule(static)
404         for(e=0;e<numElements;e++){         for(e=0;e<numElements;e++){
405             for (s=0;s<numShape; s++) {             for (s=0;s<numShape; s++) {
406               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
407               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
408                 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
409             }             }
410             for (q=0;q<numQuad;q++) {             for (q=0;q<numQuad;q++) {
411                dXdv00=0;                dXdv00=0;
412                dXdv10=0;                dXdv10=0;
413                  dXdv20=0;
414                  dXdv01=0;
415                  dXdv11=0;
416                  dXdv21=0;
417                for (s=0;s<numShape; s++) {                for (s=0;s<numShape; s++) {
418                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
419                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
420                     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
421                     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
422                     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
423                     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
424                }                }
425                D=dXdv00*dXdv00+dXdv10*dXdv10;                m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
426                  m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
427                  m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
428                  D=m00*m11-m01*m01;
429                if (D==0.) {                if (D==0.) {
430                    sprintf(error_msg,"Assemble_jacobeans_2D_M1D: element %d (id %d) has length zero.",e,element_id[e]);                    sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
431                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
432                } else {                } else {
433                   invD=1./D;                   invD=1./D;
434                   dvdX00=dXdv00*invD;                   dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
435                   dvdX01=dXdv10*invD;                   dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
436                     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
437                     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
438                     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
439                     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
440                   for (s=0;s<numTest; s++) {                   for (s=0;s<numTest; s++) {
441                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10;
442                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11;
443                       dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12;
444                   }                   }
445               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
446                }                }

Legend:
Removed from v.776  
changed lines
  Added in v.777

  ViewVC Help
Powered by ViewVC 1.1.26