/[escript]/branches/trilinos_from_5897/dudley/src/ShapeTable.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/ShapeTable.cpp

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

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 15  Line 15 
15    
16  #include "ShapeTable.h"  #include "ShapeTable.h"
17    
 #include <cstdlib>  
   
18  namespace dudley {  namespace dudley {
19    
20  /* Joel Fenwick - derived from info in Finley's Quadrature and shape files  // Joel Fenwick - derived from info in Finley's Quadrature and shape files
21    
22  This method is not threadsafe unless the initial call has completed  // This method is not thread-safe unless the initial call has completed.
23  Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions  // Evaluates the shape functions at nodes (This is the S value from the finley
24  The dim argument is the dimension of the element not the dimension of the embedding space.  // ShapeFunctions). The dim argument is the dimension of the element not the
25  the reduced arg is whether the elements are reduced or not  // dimension of the embedding space. The reduced arg is whether the elements
26  */  // are reduced or not
27  bool getQuadShape(dim_t dim, bool reduced, const double **shapearr)  bool getQuadShape(dim_t dim, bool reduced, const double **shapearr)
28  {  {
29  #define _dudley_s_alpha 0.58541019662496852  #define _dudley_s_alpha 0.58541019662496852
30  #define _dudley_s_beta  0.1381966011250105  #define _dudley_s_beta  0.1381966011250105
31    
32  /* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */      // {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints
33      static const double _dudley_V[3 * 2][12] = {      static const double _dudley_V[3 * 2][12] = {
34          {0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line single */          {0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // Line single
35          {(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},     /* Line 2 points */          {(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},     // Line 2 points
36          {1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tri single */          {1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // Tri single
37          {0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0},   /* Tri 3 points */          {0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0},   // Tri 3 points
38          {0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0},  /* Tet single */          {0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0},  // Tet single
39          {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta,          {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta,
40           _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta,           _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta,
41           _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta,           _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta,
42           _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha}       /* Tet 4 points */           _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha} // Tet 4 points
43      };      };
44    
45  #undef _dudley_s_alpha  #undef _dudley_s_alpha
46  #undef _dudley_s_beta  #undef _dudley_s_beta
47    
48      static double **arr = 0;      static double **arr = NULL;
49    
50      if (arr == 0)      if (!arr) {
     {  
51          int i;          int i;
52          arr = new double*[8];   /* point occupies two slots to make things simpler */          arr = new double*[8];   // point occupies two slots to make things simpler
53          arr[0] = new double[1];          arr[0] = new double[1];
54          arr[0][0] = 1.;         /* point */          arr[0][0] = 1.;         // point
55          arr[1] = arr[0];          arr[1] = arr[0];
56          arr[2] = new double[4]; /* Line Single */          arr[2] = new double[4]; // Line Single
57          arr[3] = new double[4]; /* Line 2 */          arr[3] = new double[4]; // Line 2
58    
59  /*          for (i = 0; i < 2; ++i) {
         for (i = 0; i < 2; ++i)  
         {  
60              arr[2][2 * i] = 1 - _dudley_V[0][i];              arr[2][2 * i] = 1 - _dudley_V[0][i];
             arr[3][2 * i] = 1 - _dudley_V[1][i];  
   
61              arr[2][2 * i + 1] = _dudley_V[0][i];              arr[2][2 * i + 1] = _dudley_V[0][i];
             arr[3][2 * i + 1] = _dudley_V[1][i];  
         }  
 */  
   
         for (i = 0; i < 2; ++i)  
         {  
             arr[2][2 * i] = 1 - _dudley_V[0][i];  
             arr[2][2 * i + 1] = _dudley_V[0][i];  
         }  
         for (i = 0; i < 2; ++i)  
         {  
62              arr[3][2 * i] = 1 - _dudley_V[1][i];              arr[3][2 * i] = 1 - _dudley_V[1][i];
63              arr[3][2 * i + 1] = _dudley_V[1][i];              arr[3][2 * i + 1] = _dudley_V[1][i];
64          }          }
65    
66            arr[4] = new double[3]; // Tri single
   
         arr[4] = new double[3]; /* Tri single */  
67          arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];          arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
68          arr[4][1] = _dudley_V[2][0];          arr[4][1] = _dudley_V[2][0];
69          arr[4][2] = _dudley_V[2][1];          arr[4][2] = _dudley_V[2][1];
70    
71          arr[5] = new double[9]; /* Tri 3 */          arr[5] = new double[9]; // Tri 3
72          for (i = 0; i < 3; ++i)          for (i = 0; i < 3; ++i) {
         {  
73              arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];              arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
74              arr[5][3 * i + 1] = _dudley_V[3][2 * i];              arr[5][3 * i + 1] = _dudley_V[3][2 * i];
75              arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];              arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
76          }          }
77          arr[6] = new  double[4];        /* Tet single */          arr[6] = new double[4]; // Tet single
78          arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];          arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
79          arr[6][1] = _dudley_V[4][0];          arr[6][1] = _dudley_V[4][0];
80          arr[6][2] = _dudley_V[4][1];          arr[6][2] = _dudley_V[4][1];
81          arr[6][3] = _dudley_V[4][2];          arr[6][3] = _dudley_V[4][2];
82    
83          arr[7] = new double[16];        /* Tet 4 */          arr[7] = new double[16]; // Tet 4
84          for (i = 0; i < 4; ++i)          for (i = 0; i < 4; ++i) {
85          {              const double x = _dudley_V[5][3 * i];
86              double x = _dudley_V[5][3 * i];              const double y = _dudley_V[5][3 * i + 1];
87              double y = _dudley_V[5][3 * i + 1];              const double z = _dudley_V[5][3 * i + 2];
             double z = _dudley_V[5][3 * i + 2];  
88              arr[7][4 * i] = 1 - x - y - z;              arr[7][4 * i] = 1 - x - y - z;
89              arr[7][4 * i + 1] = x;              arr[7][4 * i + 1] = x;
90              arr[7][4 * i + 2] = y;              arr[7][4 * i + 2] = y;
91              arr[7][4 * i + 3] = z;              arr[7][4 * i + 3] = z;
92          }          }
93      }                           /* end if */      } // end if
94    
95      if ((dim > -1) && (dim < 4))      if (dim > -1 && dim < 4) {
     {  
96          *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];          *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
97          return 1;          return 1;
98      }      }
99      *shapearr = 0;      *shapearr = NULL;
100      return 0;      return 0;
101  }  }
102    
103  const char *getElementName(Dudley_ElementTypeId id)  const char *getElementName(ElementTypeId id)
104  {  {
105      switch (id)      switch (id) {
106      {          case Dudley_Point1:
107      case Dudley_Point1:              return "Point1";
108          return "Point1";          case Dudley_Line2:
109      case Dudley_Line2:              return "Line2";
110          return "Line2";          case Dudley_Tri3:
111      case Dudley_Tri3:              return "Tri3";
112          return "Tri3";          case Dudley_Tet4:
113      case Dudley_Tet4:              return "Tet4";
114          return "Tet4";          case Dudley_Line2Face:
115      case Dudley_Line2Face:              return "Line2Face";
116          return "Line2Face";          case Dudley_Tri3Face:
117      case Dudley_Tri3Face:              return "Tri3Face";
118          return "Tri3Face";          case Dudley_Tet4Face:
119      case Dudley_Tet4Face:              return "Tet4Face";
120          return "Tet4Face";          default:
121      default:              return "noElement";
         return "noElement";  
122      }      }
123  }  }
124    

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

  ViewVC Help
Powered by ViewVC 1.1.26