# Diff of /branches/trilinos_from_5897/dudley/src/Assemble_jacobians.cpp

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 18  Line 18
18  #include "ShapeTable.h"  #include "ShapeTable.h"
19  #include "Util.h"  #include "Util.h"
20
21  /* Unless the loops in here get complicated again, this file should be compiled with loop unrolling */  // Unless the loops in here get complicated again, this file should be
22    // compiled with loop unrolling
23
24  /* input:  /* input:
25
# Line 32  index_t* nodes[numNodes*numElements]  wh Line 33  index_t* nodes[numNodes*numElements]  wh
34  dim_t numTest  dim_t numTest
36  index_t* element_id[numElements]  index_t* elementId[numElements]
37
38  output:  output:
39
47  namespace dudley {  namespace dudley {
48
49  /****************************************************************************/  /****************************************************************************/
50  /*                                                            */  //
51  /*  Jacobian 2D with area element                             */  //  Jacobian 2D with area element
52  /*                                                            */  //
53  void Assemble_jacobians_2D(double* coordinates, dim_t numQuad,  void Assemble_jacobians_2D(const double* coordinates, int numQuad,
54                           dim_t numElements, dim_t numNodes, index_t* nodes,                         dim_t numElements, int numNodes, const index_t* nodes,
55                           double* dTdX, double* absD, double* quadweight,                         double* dTdX, double* absD, double* quadweight,
56                           index_t* element_id)                         const index_t* elementId)
57  {  {
58      const int DIM = 2;      const int DIM = 2;
59      int e, q;      const int numTest = 3; // hoping this is used in constant folding
60      const dim_t numTest = 3;    /* hoping this is used in constant folding */      *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6; // numQuad is 1 or 3
61      *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6;     /* numQuad is 1 or 3 */  #pragma omp parallel for
62  #pragma omp parallel      for (index_t e = 0; e < numElements; e++) {
63      {  #define COMPDXDV0(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
double dXdv00, dXdv10, dXdv01, dXdv11, dvdX00, dvdX10, dvdX01, dvdX11, D, invD;
#pragma omp for private(e,q, dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD) schedule(static)
for (e = 0; e < numElements; e++)
{
#define COMPDXDV0(P)  coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
64  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
65  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
66
# Line 72  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu Line 68  coordinates[INDEX2(P,nodes[INDEX2(2,e,nu
68  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\  coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
69  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)  coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
70
71              dXdv00 = 0;          double dXdv00 = COMPDXDV0(0);
72              dXdv10 = 0;          double dXdv10 = COMPDXDV0(1);
73              dXdv01 = 0;          double dXdv01 = COMPDXDV1(0);
74              dXdv11 = 0;          double dXdv11 = COMPDXDV1(1);
75              dXdv00 = COMPDXDV0(0);          const double D = dXdv00 * dXdv11 - dXdv01 * dXdv10;
76              dXdv10 = COMPDXDV0(1);          absD[e] = std::abs(D);
77              dXdv01 = COMPDXDV1(0);          if (D == 0.) {
78              dXdv11 = COMPDXDV1(1);              std::stringstream ss;
79              D = dXdv00 * dXdv11 - dXdv01 * dXdv10;              ss << "Assemble_jacobians_2D: element " << e
80              absD[e] = std::abs(D);                  << " (id " << elementId[e] << ") has area zero.";
81              if (D == 0.) {              throw DudleyException(ss.str());
82                  std::stringstream ss;          } else {
83                  ss << "Assemble_jacobians_2D: element " << e              const double invD = 1. / D;
84                      << " (id " << element_id[e] << ") has area zero.";              const double dvdX00 = dXdv11 * invD;
85                  const std::string msg(ss.str());              const double dvdX10 = -dXdv10 * invD;
86                  throw DudleyException(msg);              const double dvdX01 = -dXdv01 * invD;
87              } else {              const double dvdX11 = dXdv00 * invD;
88                  invD = 1. / D;              if (numQuad == 1) {
89                  dvdX00 = dXdv11 * invD;                  dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] =
90                  dvdX10 = -dXdv10 * invD;                      DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;
91                  dvdX01 = -dXdv01 * invD;                  dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] =
92                  dvdX11 = dXdv00 * invD;                      DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;
93                  if (numQuad == 1) {                  dTdX[INDEX4(2, 0, 0, e, numTest, DIM, numQuad)] =
94                      dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] = DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;                      DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;
95                      dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] = DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;
96                      dTdX[INDEX4(2, 0, 0, e, numTest, DIM, numQuad)] = DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;                  dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] =
97                        DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;
98                      dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] = DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;                  dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] =
99                      dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] = DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;                      DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;
100                      dTdX[INDEX4(2, 1, 0, e, numTest, DIM, numQuad)] = DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;                  dTdX[INDEX4(2, 1, 0, e, numTest, DIM, numQuad)] =
101                        DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;
102                  } else { /* numQuad==3 */
103                      // relying on unroll loops to optimise this              } else { // numQuad == 3
104                      for (q = 0; q < numTest; ++q) {                  // relying on unroll loops to optimise this
105                          dTdX[INDEX4(0, 0, q, e, numTest, DIM, numQuad)] =                  for (int q = 0; q < numTest; ++q) {
106                              DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;                      dTdX[INDEX4(0, 0, q, e, numTest, DIM, numQuad)] =
107                          dTdX[INDEX4(1, 0, q, e, numTest, DIM, numQuad)] =                          DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;
108                              DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;                      dTdX[INDEX4(1, 0, q, e, numTest, DIM, numQuad)] =
109                          dTdX[INDEX4(2, 0, q, e, numTest, DIM, numQuad)] =                          DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;
110                              DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;                      dTdX[INDEX4(2, 0, q, e, numTest, DIM, numQuad)] =
111                            DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;
112                          dTdX[INDEX4(0, 1, q, e, numTest, DIM, numQuad)] =
113                              DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;                      dTdX[INDEX4(0, 1, q, e, numTest, DIM, numQuad)] =
114                          dTdX[INDEX4(1, 1, q, e, numTest, DIM, numQuad)] =                          DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;
115                              DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;                      dTdX[INDEX4(1, 1, q, e, numTest, DIM, numQuad)] =
116                          dTdX[INDEX4(2, 1, q, e, numTest, DIM, numQuad)] =                          DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;
117                              DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;                      dTdX[INDEX4(2, 1, q, e, numTest, DIM, numQuad)] =
118                            DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;
119
}
120                  }                  }
121              }              }
122          }          }
123      } // end parallel      } // end parallel for
#undef DTDXSET
124  #undef COMPDXDV0  #undef COMPDXDV0
125  #undef COMPDXDV1  #undef COMPDXDV1
126  }  }
127
128  /****************************************************************************/  //
129  /*                                                            */  // Jacobian 1D manifold in 2D and 1D elements
130  /*  Jacobean 1D manifold in 2D and 1D elements                */  //
131  /*                                                            */  void Assemble_jacobians_2D_M1D_E1D(const double* coordinates, int numQuad,
132  void Assemble_jacobians_2D_M1D_E1D(double *coordinates, dim_t numQuad,                                     dim_t numElements, int numNodes,
133                                     dim_t numElements, dim_t numNodes, index_t * nodes,                                     const index_t* nodes, double* dTdX,
134                                     double *dTdX, double *absD, double *quadweight, index_t * element_id)                                     double* absD, double* quadweight,
135                                       const index_t* elementId)
136  {  {
137      const int DIM = 2;      const int DIM = 2;
138      int e;      const int numTest = 2;
139      const dim_t numTest = 2;      *quadweight = (numQuad == 1) ? 1.0 : 0.5; // numQuad is 1 or 2
140      *quadweight = (numQuad == 1) ? 1.0 : 0.5;  #pragma omp parallel for
141      /* numQuad is 1 or 2 */      for (index_t e = 0; e < numElements; e++) {
142  #pragma omp parallel          double dXdv00 =
143      {              coordinates[INDEX2(0, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +
144          double dXdv00, dXdv10, dvdX00, dvdX01, D, invD;              coordinates[INDEX2(0, nodes[INDEX2(1, e, numNodes)], DIM)];
145  #pragma omp for private(e,dXdv00,dXdv10,dvdX00,dvdX01,D,invD) schedule(static)          double dXdv10 =
146          for (e = 0; e < numElements; e++) {              coordinates[INDEX2(1, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +
147              dXdv00 = 0;              coordinates[INDEX2(1, nodes[INDEX2(1, e, numNodes)], DIM)];
148              dXdv10 = 0;          const double D = dXdv00 * dXdv00 + dXdv10 * dXdv10;
149              dXdv00 +=          if (D == 0.) {
150                  coordinates[INDEX2(0, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +              std::stringstream ss;
151                  coordinates[INDEX2(0, nodes[INDEX2(1, e, numNodes)], DIM)];              ss << "Assemble_jacobians_2D_M1D_E1D: element " << e
152              dXdv00 +=                  << " (id " << elementId[e] << ") has length zero.";
153                  coordinates[INDEX2(1, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +              throw DudleyException(ss.str());
154                  coordinates[INDEX2(1, nodes[INDEX2(1, e, numNodes)], DIM)];          } else {
155              D = dXdv00 * dXdv00 + dXdv10 * dXdv10;              const double invD = 1. / D;
156              if (D == 0.) {              const double dvdX00 = dXdv00 * invD;
157                  std::stringstream ss;              const double dvdX01 = dXdv10 * invD;
158                  ss << "Assemble_jacobians_2D_M1D_E1D: element " << e              // The number of quad points is 1 or 2
159                      << " (id " << element_id[e] << ") has length zero.";              dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;
160                  const std::string msg(ss.str());              dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;
161                  throw DudleyException(msg);              dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;
162              } else {              dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;
163                  invD = 1. / D;              absD[e] = sqrt(D);
164                  dvdX00 = dXdv00 * invD;              if (numQuad == 2) {
165                  dvdX01 = dXdv10 * invD;                  dTdX[INDEX4(0, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
166                  /* The number of quad points is 1 or 2 */                  dTdX[INDEX4(0, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
167                  dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;                  dTdX[INDEX4(1, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
168                  dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;                  dTdX[INDEX4(1, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;
dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;
absD[e] = sqrt(D);
dTdX[INDEX4(0, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
dTdX[INDEX4(0, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
dTdX[INDEX4(1, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
dTdX[INDEX4(1, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
}
169              }              }
170          }          }
171      } // end parallel      } // end parallel for
172  }  }
173
174  /*                                                            */  //
175  /*  Jacobian 3D                                               */  // Jacobian 3D
176  /*                                                            */  //
177  void Assemble_jacobians_3D(double *coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t * nodes,  void Assemble_jacobians_3D(const double* coordinates, int numQuad,
178                             double *dTdX, double *absD, double *quadweight, index_t * element_id)                             dim_t numElements, int numNodes,
179                               const index_t* nodes, double* dTdX, double* absD,
180                               double* quadweight, const index_t* elementId)
181  {  {
182      const int DIM = 3;      const int DIM = 3;
183      int e, q, s;      const int numShape = 4;
184      /* numQuad is 1 or 4 */      const int numTest = 4;
185      const dim_t numShape = 4, numTest = 4;      *quadweight = (numQuad == 1) ? 1. / 6 : 1. / 24; // numQuad is 1 or 4
186      *quadweight = (numQuad == 1) ? 1. / 6 : 1. / 24;
187    #pragma omp parallel for
188  #pragma omp parallel      for (index_t e = 0; e < numElements; e++) {
189      {          double dXdv00 = 0;
190          double dXdv00, dXdv10, dXdv20, dXdv01, dXdv11, dXdv21, dXdv02, dXdv12, dXdv22,          double dXdv10 = 0;
191              dvdX00, dvdX10, dvdX20, dvdX01, dvdX11, dvdX21, dvdX02, dvdX12, dvdX22, D, invD, X0_loc, X1_loc, X2_loc;          double dXdv20 = 0;
192  #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)          double dXdv01 = 0;
193          for (e = 0; e < numElements; e++) {          double dXdv11 = 0;
194              dXdv00 = 0;          double dXdv21 = 0;
195              dXdv10 = 0;          double dXdv02 = 0;
196              dXdv20 = 0;          double dXdv12 = 0;
197              dXdv01 = 0;          double dXdv22 = 0;
198              dXdv11 = 0;          for (int s = 0; s < numShape; s++) {
199              dXdv21 = 0;              const double X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];
200              dXdv02 = 0;              const double X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];
201              dXdv12 = 0;              const double X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];
202              dXdv22 = 0;              dXdv00 += X0_loc * DTDV_3D[s][0];
203              for (s = 0; s < numShape; s++) {              dXdv10 += X1_loc * DTDV_3D[s][0];
204                  X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv20 += X2_loc * DTDV_3D[s][0];
205                  X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv01 += X0_loc * DTDV_3D[s][1];
206                  X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv11 += X1_loc * DTDV_3D[s][1];
207                  dXdv00 += X0_loc * DTDV_3D[s][0];              dXdv21 += X2_loc * DTDV_3D[s][1];
208                  dXdv10 += X1_loc * DTDV_3D[s][0];              dXdv02 += X0_loc * DTDV_3D[s][2];
209                  dXdv20 += X2_loc * DTDV_3D[s][0];              dXdv12 += X1_loc * DTDV_3D[s][2];
210                  dXdv01 += X0_loc * DTDV_3D[s][1];              dXdv22 += X2_loc * DTDV_3D[s][2];
211                  dXdv11 += X1_loc * DTDV_3D[s][1];          }
212                  dXdv21 += X2_loc * DTDV_3D[s][1];          const double D = dXdv00 * (dXdv11 * dXdv22 - dXdv12 * dXdv21)
213                  dXdv02 += X0_loc * DTDV_3D[s][2];                         + dXdv01 * (dXdv20 * dXdv12 - dXdv10 * dXdv22)
214                  dXdv12 += X1_loc * DTDV_3D[s][2];                         + dXdv02 * (dXdv10 * dXdv21 - dXdv20 * dXdv11);
215                  dXdv22 += X2_loc * DTDV_3D[s][2];          absD[e] = std::abs(D);
216              }          if (D == 0.) {
217              D = dXdv00 * (dXdv11 * dXdv22 - dXdv12 * dXdv21) + dXdv01 * (dXdv20 * dXdv12 - dXdv10 * dXdv22) +              std::stringstream ss;
218                  dXdv02 * (dXdv10 * dXdv21 - dXdv20 * dXdv11);              ss << "Assemble_jacobians_3D: element " << e
219              absD[e] = std::abs(D);                  << " (id " << elementId[e] << ") has volume zero.";
220              if (D == 0.) {              throw DudleyException(ss.str());
221                  std::stringstream ss;          } else {
222                  ss << "Assemble_jacobians_3D: element " << e              const double invD = 1. / D;
223                      << " (id " << element_id[e] << ") has volume zero.";              const double dvdX00 = (dXdv11 * dXdv22 - dXdv12 * dXdv21) * invD;
224                  const std::string msg(ss.str());              const double dvdX10 = (dXdv20 * dXdv12 - dXdv10 * dXdv22) * invD;
225                  throw DudleyException(msg);              const double dvdX20 = (dXdv10 * dXdv21 - dXdv20 * dXdv11) * invD;
226              } else {              const double dvdX01 = (dXdv02 * dXdv21 - dXdv01 * dXdv22) * invD;
227                  invD = 1. / D;              const double dvdX11 = (dXdv00 * dXdv22 - dXdv20 * dXdv02) * invD;
228                  dvdX00 = (dXdv11 * dXdv22 - dXdv12 * dXdv21) * invD;              const double dvdX21 = (dXdv01 * dXdv20 - dXdv00 * dXdv21) * invD;
229                  dvdX10 = (dXdv20 * dXdv12 - dXdv10 * dXdv22) * invD;              const double dvdX02 = (dXdv01 * dXdv12 - dXdv02 * dXdv11) * invD;
230                  dvdX20 = (dXdv10 * dXdv21 - dXdv20 * dXdv11) * invD;              const double dvdX12 = (dXdv02 * dXdv10 - dXdv00 * dXdv12) * invD;
231                  dvdX01 = (dXdv02 * dXdv21 - dXdv01 * dXdv22) * invD;              const double dvdX22 = (dXdv00 * dXdv11 - dXdv01 * dXdv10) * invD;
232                  dvdX11 = (dXdv00 * dXdv22 - dXdv20 * dXdv02) * invD;              for (int q = 0; q < numQuad; q++) {
233                  dvdX21 = (dXdv01 * dXdv20 - dXdv00 * dXdv21) * invD;                  for (int s = 0; s < numTest; s++) {
234                  dvdX02 = (dXdv01 * dXdv12 - dXdv02 * dXdv11) * invD;                      dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] =
235                  dvdX12 = (dXdv02 * dXdv10 - dXdv00 * dXdv12) * invD;                          DTDV_3D[s][0] * dvdX00 + DTDV_3D[s][1] * dvdX10
236                  dvdX22 = (dXdv00 * dXdv11 - dXdv01 * dXdv10) * invD;                          + DTDV_3D[s][2] * dvdX20;
237                  for (q = 0; q < numQuad; q++) {                      dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] =
238                      for (s = 0; s < numTest; s++) {                          DTDV_3D[s][0] * dvdX01 + DTDV_3D[s][1] * dvdX11
239                          dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] =                          + DTDV_3D[s][2] * dvdX21;
240                              DTDV_3D[s][0] * dvdX00 + DTDV_3D[s][1] * dvdX10 + DTDV_3D[s][2] * dvdX20;                      dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] =
241                          dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] =                          DTDV_3D[s][0] * dvdX02 + DTDV_3D[s][1] * dvdX12
242                              DTDV_3D[s][0] * dvdX01 + DTDV_3D[s][1] * dvdX11 + DTDV_3D[s][2] * dvdX21;                          + DTDV_3D[s][2] * dvdX22;
dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] =
DTDV_3D[s][0] * dvdX02 + DTDV_3D[s][1] * dvdX12 + DTDV_3D[s][2] * dvdX22;
}
243                  }                  }
244              }              }
245          }          }
246      } // end parallel      } // end parallel for
247  }  }
248
249  /*                                                            */  //
250  /*  Jacobian 2D manifold in 3D with 2D elements               */  // Jacobian 2D manifold in 3D with 2D elements
251  /*                                                            */  //
252  void Assemble_jacobians_3D_M2D_E2D(double *coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes,  void Assemble_jacobians_3D_M2D_E2D(const double* coordinates, int numQuad,
253                                     index_t * nodes, double *dTdX, double *absD, double *quadweight,                                     dim_t numElements, int numNodes,
254                                     index_t * element_id)                                     const index_t* nodes, double* dTdX,
256                                       const index_t* elementId)
257  {  {
258      const int DIM = 3;      const int DIM = 3;
int e, q, s;
259      const double DTDV[3][2] = { {-1., -1.}, {1., 0.}, {0., 1.} };      const double DTDV[3][2] = { {-1., -1.}, {1., 0.}, {0., 1.} };
260      const dim_t numShape = 3, numTest = 3;      const int numShape = 3;
261      /* numQuad is 1 or 3 */      const int numTest = 3;
262      *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6;      *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6; // numQuad is 1 or 3
263  #pragma omp parallel  #pragma omp parallel for
264      {      for (index_t e = 0; e < numElements; e++) {
265          double dXdv00, dXdv10, dXdv20, dXdv01, dXdv11, dXdv21, m00, m01, m11,          double dXdv00 = 0;
266              dvdX00, dvdX01, dvdX02, dvdX10, dvdX11, dvdX12, D, invD, X0_loc, X1_loc, X2_loc;          double dXdv10 = 0;
267  #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)          double dXdv20 = 0;
268          for (e = 0; e < numElements; e++) {          double dXdv01 = 0;
269              dXdv00 = 0;          double dXdv11 = 0;
270              dXdv10 = 0;          double dXdv21 = 0;
271              dXdv20 = 0;          for (int s = 0; s < numShape; s++) {
272              dXdv01 = 0;              const double X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];
273              dXdv11 = 0;              const double X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];
274              dXdv21 = 0;              const double X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];
275              for (s = 0; s < numShape; s++) {              dXdv00 += X0_loc * DTDV[s][0];
276                  X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv10 += X1_loc * DTDV[s][0];
277                  X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv20 += X2_loc * DTDV[s][0];
278                  X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];              dXdv01 += X0_loc * DTDV[s][1];
279                  dXdv00 += X0_loc * DTDV[s][0];              dXdv11 += X1_loc * DTDV[s][1];
280                  dXdv10 += X1_loc * DTDV[s][0];              dXdv21 += X2_loc * DTDV[s][1];
281                  dXdv20 += X2_loc * DTDV[s][0];          }
282                  dXdv01 += X0_loc * DTDV[s][1];          const double m00 = dXdv00 * dXdv00 + dXdv10 * dXdv10 + dXdv20 * dXdv20;
283                  dXdv11 += X1_loc * DTDV[s][1];          const double m01 = dXdv00 * dXdv01 + dXdv10 * dXdv11 + dXdv20 * dXdv21;
284                  dXdv21 += X2_loc * DTDV[s][1];          const double m11 = dXdv01 * dXdv01 + dXdv11 * dXdv11 + dXdv21 * dXdv21;
285              }          const double D = m00 * m11 - m01 * m01;
286              m00 = dXdv00 * dXdv00 + dXdv10 * dXdv10 + dXdv20 * dXdv20;          absD[e] = sqrt(D);
287              m01 = dXdv00 * dXdv01 + dXdv10 * dXdv11 + dXdv20 * dXdv21;          if (D == 0.) {
288              m11 = dXdv01 * dXdv01 + dXdv11 * dXdv11 + dXdv21 * dXdv21;              std::stringstream ss;
289              D = m00 * m11 - m01 * m01;              ss << "Assemble_jacobians_3D_M2D: element " << e
290              absD[e] = sqrt(D);                  << " (id " << elementId[e] << ") has area zero.";
291              if (D == 0.) {              throw DudleyException(ss.str());
292                  std::stringstream ss;          } else {
293                  ss << "Assemble_jacobians_3D_M2D: element " << e              const double invD = 1. / D;
294                      << " (id " << element_id[e] << ") has area zero.";              const double dvdX00 = (m00 * dXdv00 - m01 * dXdv01) * invD;
295                  const std::string msg(ss.str());              const double dvdX01 = (m00 * dXdv10 - m01 * dXdv11) * invD;
296                  throw DudleyException(msg);              const double dvdX02 = (m00 * dXdv20 - m01 * dXdv21) * invD;
297              } else {              const double dvdX10 = (-m01 * dXdv00 + m11 * dXdv01) * invD;
298                  invD = 1. / D;              const double dvdX11 = (-m01 * dXdv10 + m11 * dXdv11) * invD;
299                  dvdX00 = (m00 * dXdv00 - m01 * dXdv01) * invD;              const double dvdX12 = (-m01 * dXdv20 + m11 * dXdv21) * invD;
300                  dvdX01 = (m00 * dXdv10 - m01 * dXdv11) * invD;              for (int q = 0; q < numQuad; q++) {
301                  dvdX02 = (m00 * dXdv20 - m01 * dXdv21) * invD;                  for (int s = 0; s < numTest; s++) {
302                  dvdX10 = (-m01 * dXdv00 + m11 * dXdv01) * invD;                      dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] =
303                  dvdX11 = (-m01 * dXdv10 + m11 * dXdv11) * invD;                          DTDV[s][0] * dvdX00 + DTDV[s][1] * dvdX10;
304                  dvdX12 = (-m01 * dXdv20 + m11 * dXdv21) * invD;                      dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] =
305                  for (q = 0; q < numQuad; q++) {                          DTDV[s][0] * dvdX01 + DTDV[s][1] * dvdX11;
306                      for (s = 0; s < numTest; s++) {                      dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] =
307                          dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] = DTDV[s][0] * dvdX00 + DTDV[s][1] * dvdX10;                          DTDV[s][0] * dvdX02 + DTDV[s][1] * dvdX12;
dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] = DTDV[s][0] * dvdX01 + DTDV[s][1] * dvdX11;
dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] = DTDV[s][0] * dvdX02 + DTDV[s][1] * dvdX12;
}
308                  }                  }
309              }              }
310          }          }
311      } // end parallel section      } // end parallel for
312  }  }
313
314  } // namespace dudley  } // namespace dudley

Legend:
 Removed from v.6078 changed lines Added in v.6079