Contents of /branches/domexper/dudley/src/ShapeTable.c

Revision 3221 - (show annotations)
Wed Sep 29 01:00:21 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3478 byte(s)
```Comment stripping

```
 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 15 16 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files 17 18 This method is not threadsafe unless the initial call has completed 19 Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions 20 The dim argument is the dimension of the element not the dimension of the embedding space. 21 the reduced arg is whether the elements are reduced or not 22 */ 23 bool_t getQuadShape(dim_t dim, bool_t reduced, const double** shapearr) 24 { 25 #define _dudley_s_alpha 0.58541019662496852 26 #define _dudley_s_beta 0.1381966011250105 27 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 40 } // Tet 4 points 41 }; 42 43 #undef _dudley_s_alpha 44 #undef _dudley_s_beta 45 46 static double** arr=0; 47 48 if (arr==0) 49 { 50 arr=malloc(8*sizeof(double*)); // point occupies two slots to make things simpler 51 arr[0]=malloc(1*sizeof(double)); 52 arr[0][0]=1.; // point 53 arr[1]=arr[0]; 54 arr[2]=malloc(4*sizeof(double)); // Line Single 55 arr[3]=malloc(4*sizeof(double)); // Line 2 56 57 58 for (int 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 arr[4]=malloc(3*sizeof(double)); // Tri single 68 arr[4][0]=1.-_dudley_V[2][0]-_dudley_V[2][1]; 69 arr[4][1]=_dudley_V[2][0]; 70 arr[4][2]=_dudley_V[2][1]; 71 72 arr[5]=malloc(9*sizeof(double)); // Tri 3 73 for (int i=0;i<3;++i) 74 { 75 arr[5][3*i]=1-_dudley_V[3][2*i] -_dudley_V[3][2*i+1]; 76 arr[5][3*i+1]=_dudley_V[3][2*i]; 77 arr[5][3*i+2]=_dudley_V[3][2*i+1]; 78 79 80 fprintf(stderr, "\nGG 3 (%d)%f (%d)%f (%d)%f\n", 3*i,arr[5][3*i], 3*i+1, arr[5][3*i+1], 3*i+2, arr[5][3*i+2]); 81 } 82 arr[6]=malloc(4*sizeof(double)); // Tet single 83 arr[6][0]=1-_dudley_V[4][0] -_dudley_V[4][1] -_dudley_V[4][2]; 84 arr[6][1]=_dudley_V[4][0]; 85 arr[6][2]= _dudley_V[4][1]; 86 arr[6][3]=_dudley_V[4][2]; 87 88 arr[7]=malloc(16*sizeof(double)); // Tet 4 89 for (int i=0;i<4;++i) 90 { 91 double x=_dudley_V[5][3*i]; 92 double y=_dudley_V[5][3*i+1]; 93 double z=_dudley_V[5][3*i+2]; 94 arr[7][4*i]=1-x-y-z; 95 arr[7][4*i+1]=x; 96 arr[7][4*i+2]=y; 97 arr[7][4*i+3]=z; 98 } 99 } // end if 100 101 if ((dim>-1) && (dim<4)) 102 { 103 *shapearr=arr[(!reduced)?(2*dim+1):(2*dim)]; 104 return 1; 105 } 106 *shapearr=0; 107 return 0; 108 } 109 110 111 const char* getElementName(ElementTypeId id) 112 { 113 switch(id) 114 { 115 case Point1: return "Point1"; 116 case Line2: return "Line2"; 117 case Tri3: return "Tri3"; 118 case Tet4: return "Tet4"; 119 case Line2Face: return "Line2Face"; 120 case Tri3Face: return "Tri3Face"; 121 case Tet4Face: return "Tet4Face"; 122 default: 123 return "noElement"; 124 } 125 } 126 127 128