# Contents of /trunk/dudley/src/ShapeTable.c

Revision 3264 - (show annotations)
Tue Oct 12 04:40:06 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 4064 byte(s)
```Minor changes to avoid optimiser bug in intel-10

```
 1 /******************************************************* 2 * 3 * Copyright (c) 2003-2010 by University of Queensland 4 * Earth Systems Science Computational Center (ESSCC) 5 * http://www.uq.edu.au/esscc 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 * http://www.opensource.org/licenses/osl-3.0.php 10 * 11 *******************************************************/ 12 13 #include "ShapeTable.h" 14 #include "esysUtils/mem.h" 15 #include 16 17 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files 18 19 This method is not threadsafe unless the initial call has completed 20 Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions 21 The dim argument is the dimension of the element not the dimension of the embedding space. 22 the reduced arg is whether the elements are reduced or not 23 */ 24 bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr) 25 { 26 #define _dudley_s_alpha 0.58541019662496852 27 #define _dudley_s_beta 0.1381966011250105 28 29 /* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */ 30 static const double _dudley_V[3 * 2][12] = { 31 {0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line single */ 32 {(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line 2 points */ 33 {1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tri single */ 34 {0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0}, /* Tri 3 points */ 35 {0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tet single */ 36 {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta, 37 _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta, 38 _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta, 39 _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha} /* Tet 4 points */ 40 }; 41 42 #undef _dudley_s_alpha 43 #undef _dudley_s_beta 44 45 static double **arr = 0; 46 47 if (arr == 0) 48 { 49 int i; 50 arr = MEMALLOC(8,double*); /* point occupies two slots to make things simpler */ 51 arr[0] = MEMALLOC(1,double); 52 arr[0][0] = 1.; /* point */ 53 arr[1] = arr[0]; 54 arr[2] = MEMALLOC(4,double); /* Line Single */ 55 arr[3] = MEMALLOC(4,double); /* Line 2 */ 56 57 /* 58 for (i = 0; i < 2; ++i) 59 { 60 arr[2][2 * i] = 1 - _dudley_V[0][i]; 61 arr[3][2 * i] = 1 - _dudley_V[1][i]; 62 63 arr[2][2 * i + 1] = _dudley_V[0][i]; 64 arr[3][2 * i + 1] = _dudley_V[1][i]; 65 } 66 */ 67 68 for (i = 0; i < 2; ++i) 69 { 70 arr[2][2 * i] = 1 - _dudley_V[0][i]; 71 arr[2][2 * i + 1] = _dudley_V[0][i]; 72 } 73 for (i = 0; i < 2; ++i) 74 { 75 arr[3][2 * i] = 1 - _dudley_V[1][i]; 76 arr[3][2 * i + 1] = _dudley_V[1][i]; 77 } 78 79 80 81 arr[4] = MEMALLOC(3,double); /* Tri single */ 82 arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1]; 83 arr[4][1] = _dudley_V[2][0]; 84 arr[4][2] = _dudley_V[2][1]; 85 86 arr[5] = MEMALLOC(9,double); /* Tri 3 */ 87 for (i = 0; i < 3; ++i) 88 { 89 arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1]; 90 arr[5][3 * i + 1] = _dudley_V[3][2 * i]; 91 arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1]; 92 } 93 arr[6] = MEMALLOC(4, double); /* Tet single */ 94 arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2]; 95 arr[6][1] = _dudley_V[4][0]; 96 arr[6][2] = _dudley_V[4][1]; 97 arr[6][3] = _dudley_V[4][2]; 98 99 arr[7] = MEMALLOC(16,double); /* Tet 4 */ 100 for (i = 0; i < 4; ++i) 101 { 102 double x = _dudley_V[5][3 * i]; 103 double y = _dudley_V[5][3 * i + 1]; 104 double z = _dudley_V[5][3 * i + 2]; 105 arr[7][4 * i] = 1 - x - y - z; 106 arr[7][4 * i + 1] = x; 107 arr[7][4 * i + 2] = y; 108 arr[7][4 * i + 3] = z; 109 } 110 } /* end if */ 111 112 if ((dim > -1) && (dim < 4)) 113 { 114 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)]; 115 return 1; 116 } 117 *shapearr = 0; 118 return 0; 119 } 120 121 const char *getElementName(Dudley_ElementTypeId id) 122 { 123 switch (id) 124 { 125 case Dudley_Point1: 126 return "Point1"; 127 case Dudley_Line2: 128 return "Line2"; 129 case Dudley_Tri3: 130 return "Tri3"; 131 case Dudley_Tet4: 132 return "Tet4"; 133 case Dudley_Line2Face: 134 return "Line2Face"; 135 case Dudley_Tri3Face: 136 return "Tri3Face"; 137 case Dudley_Tet4Face: 138 return "Tet4Face"; 139 default: 140 return "noElement"; 141 } 142 }