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

Revision 3198 - (show annotations)
Thu Sep 23 05:47:27 2010 UTC (9 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 3164 byte(s)
```Phew. Moving on.
```
 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 bool_t getQuadShape(dim_t dim, bool_t reduced, const double** shapearr) 23 { 24 #define _dudley_s_alpha 0.58541019662496852 25 #define _dudley_s_beta 0.1381966011250105 26 27 28 // {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints 29 static const double _dudley_V[3*2][12] = { 30 {0.5,0,0,0,0,0,0,0,0,0,0,0}, // Line single 31 {(1.-.577350269189626)/2.,(1.+.577350269189626)/2.,0,0,0,0,0,0,0,0,0,0}, // Line 2 points 32 {1/3.,1/3.,0,0,0,0,0,0,0,0,0,0}, // Tri single 33 {0.5,0,0,0.5,0.5,0.5,0,0,0,0,0,0}, // Tri 3 points 34 {0.25,0.25,0.25,0,0,0,0,0,0,0,0,0}, // Tet single 35 {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta, 36 _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta, 37 _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta, 38 _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha 39 } // 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 arr=malloc(8*sizeof(double*)); // point occupies two slots to make things simpler 50 arr[0]=malloc(1*sizeof(double)); 51 arr[0][0]=1.; // point 52 arr[1]=arr[0]; 53 arr[2]=malloc(4*sizeof(double)); // Line Single 54 arr[3]=malloc(4*sizeof(double)); // Line 2 55 56 57 for (int i=0;i<2;++i) 58 { 59 arr[2][2*i]=1-_dudley_V[0][i]; 60 arr[3][2*i]=1-_dudley_V[1][i]; 61 62 arr[2][2*i+1]=_dudley_V[0][i]; 63 arr[3][2*i+1]=_dudley_V[1][i]; 64 } 65 66 arr[4]=malloc(3*sizeof(double)); // Tri single 67 arr[4][0]=1.-_dudley_V[2][0]-_dudley_V[2][1]; 68 arr[4][1]=_dudley_V[2][0]; 69 arr[4][2]=_dudley_V[2][1]; 70 71 arr[5]=malloc(9*sizeof(double)); // Tri 3 72 for (int i=0;i<3;++i) 73 { 74 arr[5][3*i]=1-_dudley_V[3][2*i] -_dudley_V[3][2*i+1]; 75 arr[5][3*i+1]=_dudley_V[3][2*i]; 76 arr[5][3*i+2]=_dudley_V[3][2*i+1]; 77 78 79 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]); 80 } 81 arr[6]=malloc(4*sizeof(double)); // Tet single 82 arr[6][0]=1-_dudley_V[4][0] -_dudley_V[4][1] -_dudley_V[4][2]; 83 arr[6][1]=_dudley_V[4][0]; 84 arr[6][2]= _dudley_V[4][1]; 85 arr[6][3]=_dudley_V[4][2]; 86 87 arr[7]=malloc(16*sizeof(double)); // Tet 4 88 for (int i=0;i<4;++i) 89 { 90 double x=_dudley_V[5][3*i]; 91 double y=_dudley_V[5][3*i+1]; 92 double z=_dudley_V[5][3*i+2]; 93 arr[7][4*i]=1-x-y-z; 94 arr[7][4*i+1]=x; 95 arr[7][4*i+2]=y; 96 arr[7][4*i+3]=z; 97 } 98 } // end if 99 100 if ((dim>-1) && (dim<4)) 101 { 102 *shapearr=arr[(!reduced)?(2*dim+1):(2*dim)]; 103 return 1; 104 } 105 *shapearr=0; 106 return 0; 107 }