/[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 853 by gross, Wed Sep 20 05:56:36 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 LOCDIM 1  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double D,invD, X0_loc;  
        #pragma omp for private(e,q,s,D,invD,X0_loc) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               D=0;  
               for (s=0;s<numShape; s++) {  
                  X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
                  D+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
               }  
               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,numTest,LOCDIM)]*invD;  
               }  
           volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];  
            }  
        }  
25    
26       }  double* coordinates[DIM*(*)]
27       #undef DIM  dim_t numQuad
28       #undef LOCDIM  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:
 /**************************************************************/  
 /*                                                            */  
 /*  Jacobean 2D with area element                             */  
 /*                                                            */  
 void Assemble_jacobeans_2D(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: 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,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)]=ABS(D)*QuadWeights[q];  
            }  
        }  
39    
40       }  double* dTdX[DIM*numTest*NUMSIDES*numQuad*numElements]
41       #undef DIM  double* volume[numQuad*numElements]
      #undef LOCDIM  
 }  
 /**************************************************************/  
 /*                                                            */  
 /*  Jacobean 1D manifold in 2D and 1D elements                */  
 /*                                                            */  
 void Assemble_jacobeans_2D_M1D_E1D(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 1  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD,  
                        X0_loc, X1_loc;  
        #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               dXdv00=0;  
               dXdv10=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)];  
               }  
               D=dXdv00*dXdv00+dXdv10*dXdv10;  
               if (D==0.) {  
                   sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD=1./D;  
                  dvdX00=dXdv00*invD;  
                  dvdX01=dXdv10*invD;  
                  for (s=0;s<numTest; s++) {  
                    dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00;  
                    dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01;  
                  }  
              volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];  
               }  
            }  
        }  
42    
43       }  */
44       #undef DIM  
45       #undef LOCDIM  #include "ShapeTable.h"
46  }  
47  /**************************************************************/  #define SCALING(_nsub_,_dim_) pow(1./(double)(_nsub_),1./(double)(_dim_))
 /*                                                            */  
 /*  Jacobean 1D manifold in 2D and 1D elements woth contact   */  
 /*                                                            */  
 void Assemble_jacobeans_2D_M1D_E1D_C(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 1  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00_0,dXdv10_0,dvdX00_0,dvdX01_0,D_0,invD_0,  
                        dXdv00_1,dXdv10_1,dvdX00_1,dvdX01_1,D_1,invD_1,  
                        X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1;  
        #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dvdX00_0,dvdX01_0,D_0,invD_0,dXdv00_1,dXdv10_1,dvdX00_1,dvdX01_1,D_1,invD_1,X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               dXdv00_0=0;  
               dXdv10_0=0;  
               dXdv00_1=0;  
               dXdv10_1=0;  
               for (s=0;s<numShape; s++) {  
                  X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
               }  
               D_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0;  
               D_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1;  
               if (D_0 == 0.i || D_1 == 0.) {  
                   sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD_0=1./D_0;  
                  dvdX00_0=dXdv00_0*invD_0;  
                  dvdX01_0=dXdv10_0*invD_0;  
                  dvdX00_1=dXdv00_1*invD_1;  
                  dvdX01_1=dXdv10_1*invD_1;  
                  for (s=0;s<numTest; s++) {  
                    dTdX[INDEX4(        s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0;  
                    dTdX[INDEX4(        s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0;  
                    dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1;  
                    dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1;  
                  }  
              volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];  
               }  
            }  
        }  
48    
      }  
      #undef DIM  
      #undef LOCDIM  
 }  
49  /**************************************************************/  /**************************************************************/
50  /*                                                            */  /*                                                            */
51  /*  Jacobean 1D manifold in 2D and 2D elements                */  /*  Jacobean 2D with area element                             */
52  /*                                                            */  /*                                                            */
53  void Assemble_jacobeans_2D_M1D_E2D(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 LOCDIM 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                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD,      register double dXdv00,dXdv10,dXdv01,dXdv11,
65                         X0_loc, X1_loc;                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
66         #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)
67         for(e=0;e<numElements;e++){      for(e=0;e<numElements;e++)
68             for (q=0;q<numQuad;q++) {      {
69                dXdv00=0;  #define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
70                dXdv10=0;  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
71                dXdv01=0;  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
72                dXdv11=0;  
73                for (s=0;s<numShape; s++) {  #define COMPDXDV1(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
74                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
75                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
76                   dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
77                   dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv00=0;
78                   dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];              dXdv10=0;
79                   dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];              dXdv01=0;
80                }              dXdv11=0;
81                D  =  dXdv00*dXdv11 - dXdv01*dXdv10;          dXdv00=COMPDXDV0(0);
82                if (D==0.) {          dXdv10=COMPDXDV0(1);
83                    sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);          dXdv01=COMPDXDV1(0);
84                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          dXdv11=COMPDXDV1(1);
85                } else {              D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
86                   invD=1./D;          absD[e]=ABS(D);
87                   dvdX00= dXdv11*invD;              if (D==0.)
88                   dvdX10=-dXdv10*invD;          {
89                   dvdX01=-dXdv01*invD;          sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
90                   dvdX11= dXdv00*invD;                  Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
91                }
92                   for (s=0;s<numTest; s++) {          else
93                     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;          {
94                     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;          invD=1./D;
95                   }          dvdX00= dXdv11*invD;
96                }          dvdX10=-dXdv10*invD;
97            volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];          dvdX01=-dXdv01*invD;
98             }          dvdX11= dXdv00*invD;
99         }          if (numQuad==1)
100            {
101       }              for (s=0;s<3; s++)
102       #undef DIM              {
103       #undef LOCDIM  #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
126         #undef LOCDIM
127         #undef DTDXSET
128         #undef COMPDXDV0
129         #undef COMPDXDV1
130  }  }
131  /**************************************************************/  /**************************************************************/
132  /*                                                            */  /*                                                            */
133  /*  Jacobean 1D manifold in 2D and 2D elements with contact   */  /*  Jacobean 1D manifold in 2D and 1D elements                */
134  /*                                                            */  /*                                                            */
135  void Assemble_jacobeans_2D_M1D_E2D_C(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 LOCDIM 2      #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_0,dXdv10_0,dXdv01_0,dXdv11_0,dvdX00_0,dvdX10_0,dvdX01_0,dvdX11_0, D_0,invD_0,      // numQuad is 1 or 2
146                         dXdv00_1,dXdv10_1,dXdv01_1,dXdv11_1,dvdX00_1,dvdX10_1,dvdX01_1,dvdX11_1, D_1,invD_1,      #pragma omp parallel
147                         X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1;      {
148         #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv01_0,dXdv11_0,dvdX00_0,dvdX10_0,dvdX01_0,dvdX11_0, D_0,invD_0,dXdv00_1,dXdv10_1,dXdv01_1,dXdv11_1,dvdX00_1,dvdX10_1,dvdX01_1,dvdX11_1, D_1,invD_1,X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1) schedule(static)      register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
149         for(e=0;e<numElements;e++){      #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
150             for (q=0;q<numQuad;q++) {      for(e=0;e<numElements;e++)
151                dXdv00_0=0;      {
152                dXdv10_0=0;          dXdv00=0;
153                dXdv01_0=0;          dXdv10=0;
154                dXdv11_0=0;          dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
155                dXdv00_1=0;          dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
156                dXdv10_1=0;          D=dXdv00*dXdv00+dXdv10*dXdv10;
157                dXdv01_1=0;          if (D==0.)
158                dXdv11_1=0;          {
159                for (s=0;s<numShape; s++) {          sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
160                   X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s         ,e,numNodes)],DIM)];          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
161                   X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s         ,e,numNodes)],DIM)];          } else {        
162                   X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];          invD=1./D;
163                   X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];          dvdX00=dXdv00*invD;
164                   dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dvdX01=dXdv10*invD;
165                   dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          // The number of quad points is 1 or 2
166                   dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
167                   dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
168                   dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
169                   dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
170                   dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  //      volume[INDEX2(0,e,numQuad)]=sqrt(D)*(*quadweight);
171                   dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          absD[e]=sqrt(D);
172                }          if (numQuad==2)
173                D_0  =  dXdv00_0*dXdv11_0 - dXdv01_0*dXdv10_0;          {
174                D_1  =  dXdv00_1*dXdv11_1 - dXdv01_1*dXdv10_1;              dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
175                if ( (D_0 ==0.) || (D_1 ==0.) ) {              dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
176                    sprintf(error_msg,"Assemble_jacobeans_2D_E2D_C: element %d (id %d) has area zero.",e,element_id[e]);              dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
177                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);              dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
178                } else {  //          volume[INDEX2(1,e,numQuad)]=sqrt(D)*(*quadweight);
179                   invD_0=1./D_0;          }
180                   dvdX00_0= dXdv11_0*invD_0;          }
181                   dvdX10_0=-dXdv10_0*invD_0;      }
182                   dvdX01_0=-dXdv01_0*invD_0;      }   // end parallel
183                   dvdX11_0= dXdv00_0*invD_0;      #undef DIM
184                   invD_1=1./D_1;      #undef LOCDIM
                  dvdX00_1= dXdv11_1*invD_1;  
                  dvdX10_1=-dXdv10_1*invD_1;  
                  dvdX01_1=-dXdv01_1*invD_1;  
                  dvdX11_1= dXdv00_1*invD_1;  
                  for (s=0;s<numTest; s++) {  
                    dTdX[INDEX4(        s,0,q,e,2*numTest,DIM,numQuad)]=  
                                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;  
                    dTdX[INDEX4(        s,1,q,e,2*numTest,DIM,numQuad)]=  
                                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;  
                    dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=  
                                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;  
                    dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=  
                                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;  
                  }  
               }  
           volume[INDEX2(q,e,numQuad)]=(sqrt(dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0)+sqrt(dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1))/2.*QuadWeights[q];  
            }  
        }  
   
      }  
      #undef DIM  
      #undef LOCDIM  
185  }  }
186    
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 LOCDIM 3      int e,q,s;
197       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        #pragma omp parallel
203        {
204        register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
205                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,                         dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
206                         X0_loc,X1_loc,X2_loc;                         X0_loc,X1_loc,X2_loc;
207        #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        #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)      for(e=0;e<numElements;e++)
209         for(e=0;e<numElements;e++){      {
210             for (q=0;q<numQuad;q++) {          dXdv00=0;
211                dXdv00=0;          dXdv10=0;
212                dXdv10=0;          dXdv20=0;
213                dXdv20=0;          dXdv01=0;
214                dXdv01=0;          dXdv11=0;
215                dXdv11=0;          dXdv21=0;
216                dXdv21=0;          dXdv02=0;
217                dXdv02=0;          dXdv12=0;
218                dXdv12=0;          dXdv22=0;
219                dXdv22=0;          for (s=0;s<numShape; s++)
220                for (s=0;s<numShape; s++) {          {
221                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];          X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
222                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];          X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
223                   X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];          X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
224                   dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv00+=X0_loc*DTDV_3D[s][0];
225                   dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv10+=X1_loc*DTDV_3D[s][0];
226                   dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv20+=X2_loc*DTDV_3D[s][0];
227                   dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv01+=X0_loc*DTDV_3D[s][1];
228                   dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv11+=X1_loc*DTDV_3D[s][1];
229                   dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv21+=X2_loc*DTDV_3D[s][1];
230                   dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];          dXdv02+=X0_loc*DTDV_3D[s][2];
231                   dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];          dXdv12+=X1_loc*DTDV_3D[s][2];
232                   dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];          dXdv22+=X2_loc*DTDV_3D[s][2];
233                }          }
234                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);
235                if (D==0.) {          absD[e]=ABS(D);
236                    sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);          if (D==0.)
237                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          {
238                } else {          sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
239                   invD=1./D;          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
240                   dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;          } else {
241                   dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;          invD=1./D;
242                   dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;          dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
243                   dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;          dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
244                   dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;          dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
245                   dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;          dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
246                   dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;          dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
247                   dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;          dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
248                   dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;          dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
249            dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
250                   for (s=0;s<numTest; s++) {          dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
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;          for (q=0;q<numQuad;q++)
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;          {
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;              for (s=0;s<numTest; s++)
254                   }              {
255               volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];              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 LOCDIM       #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 LOCDIM 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, m0, m1, m2,  
                        dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,  
                        X0_loc, X1_loc, X2_loc;  
        #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,X0_loc, X1_loc, X2_loc) schedule(static)  
        for(e=0;e<numElements;e++){  
            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++) {  
                  X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];  
                  X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];  
                  X2_loc=coordinates[INDEX2(2,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)];  
                  dXdv20+=X2_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)];  
                  dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
               }  
               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,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;  
                    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,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;  
                  }  
               }  
               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 LOCDIM  
 }  
 /**************************************************************/  
 /*                                                            */  
 /*  Jacobean 2D manifold in 3D with 3D elements on contact    */  
 /*                                                            */  
 void Assemble_jacobeans_3D_M2D_E3D_C(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 LOCDIM 3  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,dXdv02_0,dXdv12_0,dXdv22_0, m0_0, m1_0, m2_0,  
                        dvdX00_0,dvdX10_0,dvdX20_0,dvdX01_0,dvdX11_0,dvdX21_0,dvdX02_0,dvdX12_0,dvdX22_0, D_0,invD_0,  
                        dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,dXdv02_1,dXdv12_1,dXdv22_1, m0_1, m1_1, m2_1,  
                        dvdX00_1,dvdX10_1,dvdX20_1,dvdX01_1,dvdX11_1,dvdX21_1,dvdX02_1,dvdX12_1,dvdX22_1, D_1,invD_1,  
                        X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1;  
        #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,dXdv02_0,dXdv12_0,dXdv22_0, m0_0, m1_0, m2_0,dvdX00_0,dvdX10_0,dvdX20_0,dvdX01_0,dvdX11_0,dvdX21_0,dvdX02_0,dvdX12_0,dvdX22_0, D_0,invD_0,dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,dXdv02_1,dXdv12_1,dXdv22_1, m0_1, m1_1, m2_1,dvdX00_1,dvdX10_1,dvdX20_1,dvdX01_1,dvdX11_1,dvdX21_1,dvdX02_1,dvdX12_1,dvdX22_1, D_1,invD_1,X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               dXdv00_0=0;  
               dXdv10_0=0;  
               dXdv20_0=0;  
               dXdv01_0=0;  
               dXdv11_0=0;  
               dXdv21_0=0;  
               dXdv02_0=0;  
               dXdv12_0=0;  
               dXdv22_0=0;  
               dXdv00_1=0;  
               dXdv10_1=0;  
               dXdv20_1=0;  
               dXdv01_1=0;  
               dXdv11_1=0;  
               dXdv21_1=0;  
               dXdv02_1=0;  
               dXdv12_1=0;  
               dXdv22_1=0;  
               for (s=0;s<numShape; s++) {  
                  X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X2_loc_0=coordinates[INDEX2(2,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  X2_loc_1=coordinates[INDEX2(2,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv20_0+=X2_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv21_0+=X2_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv02_0+=X0_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv12_0+=X1_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv22_0+=X2_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv20_1+=X2_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv21_1+=X2_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv02_1+=X0_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv12_1+=X1_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
                  dXdv22_1+=X2_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];  
               }  
   
               D_0=dXdv00_0*(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)+dXdv01_0*(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)+dXdv02_0*(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0);  
               D_1=dXdv00_1*(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)+dXdv01_1*(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)+dXdv02_1*(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1);  
               if ( (D_0==0.) || (D_1 == 0.)) {  
                   sprintf(error_msg,"Assemble_jacobeans_3D_C: element %d (id %d) has volume zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD_0=1./D_0;  
                  dvdX00_0=(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)*invD_0;  
                  dvdX10_0=(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)*invD_0;  
                  dvdX20_0=(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0)*invD_0;  
                  dvdX01_0=(dXdv02_0*dXdv21_0-dXdv01_0*dXdv22_0)*invD_0;  
                  dvdX11_0=(dXdv00_0*dXdv22_0-dXdv20_0*dXdv02_0)*invD_0;  
                  dvdX21_0=(dXdv01_0*dXdv20_0-dXdv00_0*dXdv21_0)*invD_0;  
                  dvdX02_0=(dXdv01_0*dXdv12_0-dXdv02_0*dXdv11_0)*invD_0;  
                  dvdX12_0=(dXdv02_0*dXdv10_0-dXdv00_0*dXdv12_0)*invD_0;  
                  dvdX22_0=(dXdv00_0*dXdv11_0-dXdv01_0*dXdv10_0)*invD_0;  
                  invD_1=1./D_1;  
                  dvdX00_1=(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)*invD_1;  
                  dvdX10_1=(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)*invD_1;  
                  dvdX20_1=(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1)*invD_1;  
                  dvdX01_1=(dXdv02_1*dXdv21_1-dXdv01_1*dXdv22_1)*invD_1;  
                  dvdX11_1=(dXdv00_1*dXdv22_1-dXdv20_1*dXdv02_1)*invD_1;  
                  dvdX21_1=(dXdv01_1*dXdv20_1-dXdv00_1*dXdv21_1)*invD_1;  
                  dvdX02_1=(dXdv01_1*dXdv12_1-dXdv02_1*dXdv11_1)*invD_1;  
                  dvdX12_1=(dXdv02_1*dXdv10_1-dXdv00_1*dXdv12_1)*invD_1;  
                  dvdX22_1=(dXdv00_1*dXdv11_1-dXdv01_1*dXdv10_1)*invD_1;  
   
                  for (s=0;s<numTest; s++) {  
                    dTdX[INDEX4(        s,0,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_0;  
                    dTdX[INDEX4(        s,1,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_0;  
                    dTdX[INDEX4(        s,2,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_0;  
                    dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_1;  
                    dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_1;  
                    dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=  
                       DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_1;  
                  }  
               }  
               m0_0=dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0;  
               m1_0=dXdv20_0*dXdv01_0-dXdv00_0*dXdv21_0;  
               m2_0=dXdv00_0*dXdv11_0-dXdv10_0*dXdv01_0;  
               m0_1=dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1;  
               m1_1=dXdv20_1*dXdv01_1-dXdv00_1*dXdv21_1;  
               m2_1=dXdv00_1*dXdv11_1-dXdv10_1*dXdv01_1;  
           volume[INDEX2(q,e,numQuad)]=(sqrt(m0_0*m0_0+m1_0*m1_0+m2_0*m2_0)+sqrt(m0_1*m0_1+m1_1*m1_1+m2_1*m2_1))/2.*QuadWeights[q];  
            }  
        }  
   
      }  
      #undef DIM  
      #undef LOCDIM  
 }  
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 LOCDIM 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        #pragma omp parallel
284        {
285        register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
286                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,                         dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
287                         X0_loc, X1_loc, X2_loc;                         X0_loc, X1_loc, X2_loc;
288         #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)      #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         for(e=0;e<numElements;e++){      for(e=0;e<numElements;e++)
290             for (q=0;q<numQuad;q++) {      {
291                dXdv00=0;          dXdv00=0;
292                dXdv10=0;          dXdv10=0;
293                dXdv20=0;          dXdv20=0;
294                dXdv01=0;          dXdv01=0;
295                dXdv11=0;          dXdv11=0;
296                dXdv21=0;          dXdv21=0;
297                for (s=0;s<numShape; s++) {          for (s=0;s<numShape; s++)
298                   X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];          {
299                   X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];          X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
300                   X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];          X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
301                   dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
302                   dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv00+=X0_loc*DTDV[s][0];
303                   dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];          dXdv10+=X1_loc*DTDV[s][0];
304                   dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv20+=X2_loc*DTDV[s][0];
305                   dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv01+=X0_loc*DTDV[s][1];
306                   dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];          dXdv11+=X1_loc*DTDV[s][1];
307                }          dXdv21+=X2_loc*DTDV[s][1];
308                m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;          }
309                m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;          m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
310                m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;          m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
311                D=m00*m11-m01*m01;          m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
312                if (D==0.) {          D=m00*m11-m01*m01;
313                    sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);          absD[e]=sqrt(D);
314                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);          if (D==0.)
315                } else {          {
316                   invD=1./D;          sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
317                   dvdX00=( m00*dXdv00-m01*dXdv01)*invD;          Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
318                   dvdX01=( m00*dXdv10-m01*dXdv11)*invD;          } else {
319                   dvdX02=( m00*dXdv20-m01*dXdv21)*invD;          invD=1./D;
320                   dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;          dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
321                   dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;          dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
322                   dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;          dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
323                   for (s=0;s<numTest; s++) {          dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
324                     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;          dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
325                     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;          dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
326                     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;          for (q=0;q<numQuad;q++)
327                   }          {
328               volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];              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       }              }
334       #undef DIM  //          volume[INDEX2(q,e,numQuad)]=sqrt(D)*(*quadweight);
335       #undef LOCDIM          }
336  }          }
337  /**************************************************************/      }
338  /*                                                            */      }   // end parallel section
339  /*  Jacobean 2D manifold in 3D with 2D elements  with contact */      #undef DIM
340  /*                                                            */      #undef LOCDIM
 void Assemble_jacobeans_3D_M2D_E2D_C(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 LOCDIM 2  
      register int e,q,s;  
      char error_msg[LenErrorMsg_MAX];  
      #pragma omp parallel  
      {  
        register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,m00_0,m01_0,m11_0,  
                        dvdX00_0,dvdX01_0,dvdX02_0,dvdX10_0,dvdX11_0,dvdX12_0,D_0,invD_0,  
                        dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,m00_1,m01_1,m11_1,  
                        dvdX00_1,dvdX01_1,dvdX02_1,dvdX10_1,dvdX11_1,dvdX12_1,D_1,invD_1,  
                        X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1;  
        #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,m00_0,m01_0,m11_0,dvdX00_0,dvdX01_0,dvdX02_0,dvdX10_0,dvdX11_0,dvdX12_0,D_0,invD_0,dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,m00_1,m01_1,m11_1,dvdX00_1,dvdX01_1,dvdX02_1,dvdX10_1,dvdX11_1,dvdX12_1,D_1,invD_1,X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1) schedule(static)  
        for(e=0;e<numElements;e++){  
            for (q=0;q<numQuad;q++) {  
               dXdv00_0=0;  
               dXdv10_0=0;  
               dXdv20_0=0;  
               dXdv01_0=0;  
               dXdv11_0=0;  
               dXdv21_0=0;  
               dXdv00_1=0;  
               dXdv10_1=0;  
               dXdv20_1=0;  
               dXdv01_1=0;  
               dXdv11_1=0;  
               dXdv21_1=0;  
               for (s=0;s<numShape; s++) {  
                  X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X2_loc_0=coordinates[INDEX2(2,nodes[INDEX2(s         ,e,numNodes)],DIM)];  
                  X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  X2_loc_1=coordinates[INDEX2(2,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];  
                  dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv20_0+=X2_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv21_0+=X2_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv20_1+=X2_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];  
                  dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
                  dXdv21_1+=X2_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];  
               }  
               m00_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0+dXdv20_0*dXdv20_0;  
               m01_0=dXdv00_0*dXdv01_0+dXdv10_0*dXdv11_0+dXdv20_0*dXdv21_0;  
               m11_0=dXdv01_0*dXdv01_0+dXdv11_0*dXdv11_0+dXdv21_0*dXdv21_0;  
               D_0=m00_0*m11_0-m01_0*m01_0;  
               m00_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1+dXdv20_1*dXdv20_1;  
               m01_1=dXdv00_1*dXdv01_1+dXdv10_1*dXdv11_1+dXdv20_1*dXdv21_1;  
               m11_1=dXdv01_1*dXdv01_1+dXdv11_1*dXdv11_1+dXdv21_1*dXdv21_1;  
               D_1=m00_1*m11_1-m01_1*m01_1;  
               if ( (D_0==0.) || (D_1 == 0.) ) {  
                   sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);  
                   Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
               } else {  
                  invD_0=1./D_0;  
                  dvdX00_0=( m00_0*dXdv00_0-m01_0*dXdv01_0)*invD_0;  
                  dvdX01_0=( m00_0*dXdv10_0-m01_0*dXdv11_0)*invD_0;  
                  dvdX02_0=( m00_0*dXdv20_0-m01_0*dXdv21_0)*invD_0;  
                  dvdX10_0=(-m01_0*dXdv00_0+m11_0*dXdv01_0)*invD_0;  
                  dvdX11_0=(-m01_0*dXdv10_0+m11_0*dXdv11_0)*invD_0;  
                  dvdX12_0=(-m01_0*dXdv20_0+m11_0*dXdv21_0)*invD_0;  
                  invD_1=1./D_1;  
                  dvdX00_1=( m00_1*dXdv00_1-m01_1*dXdv01_1)*invD_1;  
                  dvdX01_1=( m00_1*dXdv10_1-m01_1*dXdv11_1)*invD_1;  
                  dvdX02_1=( m00_1*dXdv20_1-m01_1*dXdv21_1)*invD_1;  
                  dvdX10_1=(-m01_1*dXdv00_1+m11_1*dXdv01_1)*invD_1;  
                  dvdX11_1=(-m01_1*dXdv10_1+m11_1*dXdv11_1)*invD_1;  
                  dvdX12_1=(-m01_1*dXdv20_1+m11_1*dXdv21_1)*invD_1;  
                  for (s=0;s<numTest; s++) {  
                    dTdX[INDEX4(        s,0,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;  
                    dTdX[INDEX4(        s,1,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;  
                    dTdX[INDEX4(        s,2,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0;  
                    dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;  
                    dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;  
                    dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=  
                                  DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1;  
                  }  
              volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];  
               }  
            }  
        }  
   
      }  
      #undef DIM  
      #undef LOCDIM  
341  }  }
 /*  
  * $Log:$  
  *  
  */  

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

  ViewVC Help
Powered by ViewVC 1.1.26