/[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

trunk/finley/src/Assemble_jacobeans.c revision 777 by gross, Wed Jul 12 08:54:45 2006 UTC branches/domexper/dudley/src/Assemble_jacobeans.c revision 3189 by jfenwick, Thu Sep 16 05:20:34 2010 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
   
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    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$ */  
   
 /**************************************************************/  
14    
15  #include "Assemble.h"  #include "Assemble.h"
16  #include "Util.h"  #include "Util.h"
# Line 24  Line 18 
18  #include <omp.h>  #include <omp.h>
19  #endif  #endif
20    
21    // Unless the loops in here get complicated again, this file should be compiled with loop unrolling
22    
23    
24  /**************************************************************/  /* input:
 /*                                                            */  
 /*  Jacobean 1D                                               */  
 /*                                                            */  
 void Assemble_jacobeans_1D(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 1  
      #define LOCALDIM 1  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double D,invD;  
        double X0[numShape];  
        #pragma omp for private(e,q,s) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
            for (q=0;q<numQuad;q++) {  
               D=0;  
               for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];  
               if (D==0.) {  
                   sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD=1./D;  
                  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];  
            }  
        }  
25    
26       }  double* coordinates[DIM*(*)]
27       #undef DIM  dim_t numQuad
28       #undef LOCALDIM  double* QuadWeights[numQuad]
29    dim_t numShape
30    dim_t numElements
31    dim_t numNodes
32    index_t* nodes[numNodes*numElements]  where NUMSIDES*numShape<=numNodes
33    double* DSDv[numShape*DIM*numQuad]
34    dim_t numTest
35    double* DTDv[LOCDIM*numTest*numQuad]
36    index_t* element_id[numElements]
37    
38    output:
39    
40    double* dTdX[DIM*numTest*NUMSIDES*numQuad*numElements]
41    double* volume[numQuad*numElements]
42    
43    */
44    
45    #include "ShapeTable.h"
46    
47    #define SCALING(_nsub_,_dim_) pow(1./(double)(_nsub_),1./(double)(_dim_))
48    
 }  
49  /**************************************************************/  /**************************************************************/
50  /*                                                            */  /*                                                            */
51  /*  Jacobean 2D with area element                             */  /*  Jacobean 2D with area element                             */
52  /*                                                            */  /*                                                            */
53  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,
54                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,                             double* dTdX, double* absD, double* quadweight, index_t* element_id)
55                             double* DSDv, dim_t numTest,double* DTDv,  {
56                             double* dTdX, double* volume, index_t* element_id) {      #define DIM 2
57       #define DIM 2      #define LOCDIM 2
58       #define LOCALDIM 2      register int e,q,s;
59       register int e,q,s;      char error_msg[LenErrorMsg_MAX];
60       char error_msg[LenErrorMsg_MAX];      *quadweight=(numQuad==1)?1./2:1./6; /* numQuad is 1 or 3 */
61       #pragma omp parallel      const dim_t numTest=3;      // hoping this is used in constant folding
62       {      #pragma omp parallel
63         register double dXdv00,dXdv10,dXdv01,dXdv11,      {
64        register double dXdv00,dXdv10,dXdv01,dXdv11,
65                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD;                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
66         double X0[numShape], X1[numShape];      #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
67         #pragma omp for private(e,q,s) schedule(static)      for(e=0;e<numElements;e++)
68         for(e=0;e<numElements;e++){      {
69             for (s=0;s<numShape; s++) {  #define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
70               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
71               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
72             }  
73             for (q=0;q<numQuad;q++) {  #define COMPDXDV1(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
74                dXdv00=0;  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
75                dXdv10=0;  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
76                dXdv01=0;  
77                dXdv11=0;          dXdv00=0;
78                for (s=0;s<numShape; s++) {              dXdv10=0;
79                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];              dXdv01=0;
80                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];              dXdv11=0;
81                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv00=COMPDXDV0(0);
82                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv10=COMPDXDV0(1);
83                }          dXdv01=COMPDXDV1(0);
84                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;          dXdv11=COMPDXDV1(1);
85                if (D==0.) {              D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
86                    sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);          absD[e]=ABS(D);
87                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);              if (D==0.)
88                } else {          {
89                   invD=1./D;          sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
90                   dvdX00= dXdv11*invD;                  Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
91                   dvdX10=-dXdv10*invD;              }
92                   dvdX01=-dXdv01*invD;          else
93                   dvdX11= dXdv00*invD;          {
94            invD=1./D;
95                   for (s=0;s<numTest; s++) {          dvdX00= dXdv11*invD;
96                     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;          dvdX10=-dXdv10*invD;
97                     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;          dvdX01=-dXdv01*invD;
98                   }          dvdX11= dXdv00*invD;
99                }          if (numQuad==1)
100                volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];          {
101             }              for (s=0;s<3; s++)
102         }              {
103    #define DTDXSET(P,Q) dTdX[INDEX4(s,P,Q,e,numTest,DIM,numQuad)] = DTDV_2D[s][0]*dvdX0##P+DTDV_2D[s][1]*dvdX1##P
104       }  
105                DTDXSET(0,0);
106                DTDXSET(1,0);
107                }
108    //          volume[INDEX2(0,e,1)]=ABS(D)*(*quadweight);
109            }
110            else    // numQuad==3
111            {
112                for (q=0;q<numTest;++q) // relying on unroll loops to optimise this
113                {
114                for (s=0;s<3; s++)
115                {
116                    DTDXSET(0,q);
117                    DTDXSET(1,q);
118                }
119    //          volume[INDEX2(q,e,3)]=ABS(D)*(*quadweight);
120                }
121            }
122            }
123        }
124         }  // end parallel
125       #undef DIM       #undef DIM
126       #undef LOCALDIM       #undef LOCDIM
127         #undef DTDXSET
128         #undef COMPDXDV0
129         #undef COMPDXDV1
130  }  }
131  /**************************************************************/  /**************************************************************/
132  /*                                                            */  /*                                                            */
133  /*  Jacobean 1D manifold in 2D and 1D elements                */  /*  Jacobean 1D manifold in 2D and 1D elements                */
134  /*                                                            */  /*                                                            */
135  void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,
136                                     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,                                     dim_t numElements, dim_t numNodes, index_t* nodes,
137                                     double* DSDv, dim_t numTest,double* DTDv,                                     double* dTdX, double* absD, double* quadweight, index_t* element_id)
138                                     double* dTdX, double* volume, index_t* element_id) {  {
139       #define DIM 2      #define DIM 2
140       #define LOCALDIM 1      #define LOCDIM 1
141       register int e,q,s;      register int e;
142       char error_msg[LenErrorMsg_MAX];      char error_msg[LenErrorMsg_MAX];
143       #pragma omp parallel      const dim_t numTest=2;
144       {      *quadweight=(numQuad==1)?1.0:0.5;
145         register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;      // numQuad is 1 or 2
146         double X0[numShape], X1[numShape];      #pragma omp parallel
147         #pragma omp for private(e,q,s) schedule(static)      {
148         for(e=0;e<numElements;e++){      register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
149             for (s=0;s<numShape; s++) {      #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
150               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];      for(e=0;e<numElements;e++)
151               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];      {
152             }          dXdv00=0;
153             for (q=0;q<numQuad;q++) {          dXdv10=0;
154                dXdv00=0;          dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
155                dXdv10=0;          dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
156                for (s=0;s<numShape; s++) {          D=dXdv00*dXdv00+dXdv10*dXdv10;
157                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          if (D==0.)
158                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          {
159                }          sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
160                D=dXdv00*dXdv00+dXdv10*dXdv10;          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
161                if (D==0.) {          } else {        
162                    sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);          invD=1./D;
163                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          dvdX00=dXdv00*invD;
164                } else {          dvdX01=dXdv10*invD;
165                   invD=1./D;          // The number of quad points is 1 or 2
166                   dvdX00=dXdv00*invD;          dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
167                   dvdX01=dXdv10*invD;          dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
168                   for (s=0;s<numTest; s++) {          dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
169                     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00;          dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
170                     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01;  //      volume[INDEX2(0,e,numQuad)]=sqrt(D)*(*quadweight);
171                   }          absD[e]=sqrt(D);
172               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];          if (numQuad==2)
173                }          {
174             }              dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
175         }              dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
176                dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
177       }              dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
178       #undef DIM  //          volume[INDEX2(1,e,numQuad)]=sqrt(D)*(*quadweight);
179       #undef LOCALDIM          }
180            }
181        }
182        }   // end parallel
183        #undef DIM
184        #undef LOCDIM
185  }  }
 /**************************************************************/  
 /*                                                            */  
 /*  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 LOCALDIM 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;  
        double X0[numShape], X1[numShape];  
        #pragma omp for private(e,q,s) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (s=0;s<numShape; s++) {  
              X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
              X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  
            }  
            for (q=0;q<numQuad;q++) {  
               dXdv00=0;  
               dXdv10=0;  
               dXdv01=0;  
               dXdv11=0;  
               for (s=0;s<numShape; s++) {  
                  dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];  
                  dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];  
                  dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];  
                  dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];  
               }  
               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]);  
                   Finley_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,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10;  
                    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;  
                  }  
               }  
           volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];  
            }  
        }  
186    
      }  
      #undef DIM  
      #undef LOCALDIM  
 }  
187  /**************************************************************/  /**************************************************************/
188  /*                                                            */  /*                                                            */
189  /*  Jacobean 3D                                               */  /*  Jacobean 3D                                               */
190  /*                                                            */  /*                                                            */
191  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,
192                             dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,                             double* dTdX, double* absD, double* quadweight, index_t* element_id)
193                             double* DSDv, dim_t numTest,double* DTDv,  {
194                             double* dTdX, double* volume, index_t* element_id) {      #define DIM 3
195       #define DIM 3      #define LOCDIM 3
196       #define LOCALDIM 3      int e,q,s;
197       register int e,q,s;      char error_msg[LenErrorMsg_MAX];
198       char error_msg[LenErrorMsg_MAX];      // numQuad is 1 or 4
199       #pragma omp parallel      const dim_t numShape=4, numTest=4;
200       {      *quadweight=(numQuad==1)?1./6:1./24;
201         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,  
202                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;      #pragma omp parallel
203         double X0[numShape], X1[numShape], X2[numShape];      {
204         #pragma omp for private(e,q,s) schedule(static)      register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
205         for(e=0;e<numElements;e++){                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
206             for (s=0;s<numShape; s++) {                         X0_loc,X1_loc,X2_loc;
207               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];      #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)
208               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];      for(e=0;e<numElements;e++)
209               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];      {
210             }          dXdv00=0;
211             for (q=0;q<numQuad;q++) {          dXdv10=0;
212                dXdv00=0;          dXdv20=0;
213                dXdv10=0;          dXdv01=0;
214                dXdv20=0;          dXdv11=0;
215                dXdv01=0;          dXdv21=0;
216                dXdv11=0;          dXdv02=0;
217                dXdv21=0;          dXdv12=0;
218                dXdv02=0;          dXdv22=0;
219                dXdv12=0;          for (s=0;s<numShape; s++)
220                dXdv22=0;          {
221                for (s=0;s<numShape; s++) {          X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
222                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
223                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
224                   dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          dXdv00+=X0_loc*DTDV_3D[s][0];
225                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv10+=X1_loc*DTDV_3D[s][0];
226                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv20+=X2_loc*DTDV_3D[s][0];
227                   dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv01+=X0_loc*DTDV_3D[s][1];
228                   dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];          dXdv11+=X1_loc*DTDV_3D[s][1];
229                   dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];          dXdv21+=X2_loc*DTDV_3D[s][1];
230                   dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];          dXdv02+=X0_loc*DTDV_3D[s][2];
231                }          dXdv12+=X1_loc*DTDV_3D[s][2];
232                D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);          dXdv22+=X2_loc*DTDV_3D[s][2];
233                if (D==0.) {          }
234                    sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);          D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
235                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          absD[e]=ABS(D);
236                } else {          if (D==0.)
237                   invD=1./D;          {
238                   dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;          sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
239                   dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
240                   dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;          } else {
241                   dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;          invD=1./D;
242                   dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;          dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
243                   dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;          dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
244                   dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;          dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
245                   dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;          dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
246                   dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;          dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
247            dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
248                   for (s=0;s<numTest; s++) {          dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
249                     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;          dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
250                     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;          dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
251                     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;          for (q=0;q<numQuad;q++)
252                   }          {
253               volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];              for (s=0;s<numTest; s++)
254                }              {
255             }              dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX00+DTDV_3D[s][1]*dvdX10+DTDV_3D[s][2]*dvdX20;
256         }              dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX01+DTDV_3D[s][1]*dvdX11+DTDV_3D[s][2]*dvdX21;
257                dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX02+DTDV_3D[s][1]*dvdX12+DTDV_3D[s][2]*dvdX22;
258       }              }
259    //          volume[INDEX2(q,e,numQuad)]=ABS(D)*(*quadweight);
260                    }
261            }
262        }
263         }  // end parallel
264       #undef DIM       #undef DIM
265       #undef LOCALDIM       #undef LOCDIM
266  }  }
 /**************************************************************/  
 /*                                                            */  
 /*  Jacobean 2D manifold in 3D with 3D elements               */  
 /*                                                            */  
 void Assemble_jacobeans_3D_M2D_E3D(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 3  
      #define LOCALDIM 3  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,  
                        dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;  
        double X0[numShape], X1[numShape], X2[numShape], m0, m1, m2;  
        #pragma omp for private(e,q,s) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (s=0;s<numShape; s++) {  
              X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
              X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  
              X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];  
            }  
            for (q=0;q<numQuad;q++) {  
               dXdv00=0;  
               dXdv10=0;  
               dXdv20=0;  
               dXdv01=0;  
               dXdv11=0;  
               dXdv21=0;  
               dXdv02=0;  
               dXdv12=0;  
               dXdv22=0;  
               for (s=0;s<numShape; s++) {  
                  dXdv00+=X0[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)];  
                  dXdv01+=X0[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)];  
               }  
               D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);  
               if (D==0.) {  
                   sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD=1./D;  
                  dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;  
                  dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;  
                  dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;  
                  dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*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;  
   
                  for (s=0;s<numTest; s++) {  
                    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;  
                    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;  
                    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;  
                  }  
               }  
               m0=dXdv10*dXdv21-dXdv20*dXdv11;  
               m1=dXdv20*dXdv01-dXdv00*dXdv21;  
               m2=dXdv00*dXdv11-dXdv10*dXdv01;  
           volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];  
            }  
        }  
267    
      }  
      #undef DIM  
      #undef LOCALDIM  
 }  
268  /**************************************************************/  /**************************************************************/
269  /*                                                            */  /*                                                            */
270  /*  Jacobean 2D manifold in 3D with 2D elements               */  /*  Jacobean 2D manifold in 3D with 2D elements               */
271  /*                                                            */  /*                                                            */
272  void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,  void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
273                                     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,                                     double* dTdX, double* absD, double* quadweight, index_t* element_id)
274                                     double* DSDv, dim_t numTest,double* DTDv,  {
275                                     double* dTdX, double* volume, index_t* element_id) {      #define DIM 3
276       #define DIM 3      #define LOCDIM 2
277       #define LOCALDIM 2      register int e,q,s;
278       register int e,q,s;      char error_msg[LenErrorMsg_MAX];
279       char error_msg[LenErrorMsg_MAX];      // numQuad is 1 or 3
280       #pragma omp parallel      *quadweight=(numQuad==1)?1./2:1./6;
281       {      const double DTDV[3][2]={{-1.,-1.},{1.,0.},{0.,1.}};
282         register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,      const dim_t numShape=3, numTest=3;
283                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;      #pragma omp parallel
284         double X0[numShape], X1[numShape], X2[numShape];      {
285         #pragma omp for private(e,q,s) schedule(static)      register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
286         for(e=0;e<numElements;e++){                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
287             for (s=0;s<numShape; s++) {                         X0_loc, X1_loc, X2_loc;
288               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];      #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD, X0_loc, X1_loc, X2_loc) schedule(static)
289               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];      for(e=0;e<numElements;e++)
290               X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];      {
291             }          dXdv00=0;
292             for (q=0;q<numQuad;q++) {          dXdv10=0;
293                dXdv00=0;          dXdv20=0;
294                dXdv10=0;          dXdv01=0;
295                dXdv20=0;          dXdv11=0;
296                dXdv01=0;          dXdv21=0;
297                dXdv11=0;          for (s=0;s<numShape; s++)
298                dXdv21=0;          {
299                for (s=0;s<numShape; s++) {          X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
300                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
301                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
302                   dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];          dXdv00+=X0_loc*DTDV[s][0];
303                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv10+=X1_loc*DTDV[s][0];
304                   dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv20+=X2_loc*DTDV[s][0];
305                   dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];          dXdv01+=X0_loc*DTDV[s][1];
306                }          dXdv11+=X1_loc*DTDV[s][1];
307                m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;          dXdv21+=X2_loc*DTDV[s][1];
308                m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;          }
309                m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;          m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
310                D=m00*m11-m01*m01;          m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
311                if (D==0.) {          m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
312                    sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);          D=m00*m11-m01*m01;
313                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          absD[e]=sqrt(D);
314                } else {          if (D==0.)
315                   invD=1./D;          {
316                   dvdX00=( m00*dXdv00-m01*dXdv01)*invD;          sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
317                   dvdX01=( m00*dXdv10-m01*dXdv11)*invD;          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
318                   dvdX02=( m00*dXdv20-m01*dXdv21)*invD;          } else {
319                   dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;          invD=1./D;
320                   dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;          dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
321                   dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;          dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
322                   for (s=0;s<numTest; s++) {          dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
323                     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;          dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
324                     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;          dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
325                     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;          dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
326                   }          for (q=0;q<numQuad;q++)
327               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];          {
328                }              for (s=0;s<numTest; s++)
329             }              {
330         }              dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10;
331                dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11;
332       }              dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12;
333       #undef DIM              }
334       #undef LOCALDIM  //          volume[INDEX2(q,e,numQuad)]=sqrt(D)*(*quadweight);
335            }
336            }
337        }
338        }   // end parallel section
339        #undef DIM
340        #undef LOCDIM
341  }  }
 /*  
  * $Log$  
  *  
  */  

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

  ViewVC Help
Powered by ViewVC 1.1.26