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

Revision 3196 - (show annotations)
Wed Sep 22 01:18:52 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 2990 byte(s)
moving slowly
 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 arr[4]=malloc(4*sizeof(double)); // Tri single 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 arr[4][2*i]=1-_dudley_V[2][i]; 62 63 arr[2][2*i+1]=_dudley_V[0][i]; 64 arr[3][2*i+1]=_dudley_V[1][i]; 65 arr[4][2*i+1]=_dudley_V[2][i]; 66 } 67 68 arr[5]=malloc(9*sizeof(double)); // Tri 3 69 for (int i=0;i<3;++i) 70 { 71 arr[5][3*i]=1-_dudley_V[3][3*i] -_dudley_V[2][3*i+1]; 72 arr[5][3*i+1]=_dudley_V[3][3*i]; 73 arr[5][3*i+2]=_dudley_V[3][3*i+1]; 74 } 75 arr[6]=malloc(3*sizeof(double)); // Tet single 76 arr[6][0]=1-_dudley_V[4][0] -_dudley_V[2][1] -_dudley_V[2][2]; 77 arr[6][1]=_dudley_V[4][0]; 78 arr[6][2]= _dudley_V[4][1]; 79 arr[6][3]=_dudley_V[4][2]; 80 81 arr[7]=malloc(16*sizeof(double)); 82 for (int i=0;i<4;++i) 83 { 84 double x=_dudley_V[5][0]; 85 double y=_dudley_V[5][1]; 86 double z=_dudley_V[5][2]; 87 arr[7][4*i]=1-x-y-z; 88 arr[7][4*i+1]=x; 89 arr[7][4*i+2]=y; 90 arr[7][4*i+3]=z; 91 } 92 } // end if 93 94 if ((dim>-1) && (dim<4)) 95 { 96 *shapearr=arr[(!reduced)?(2*dim+1):(2*dim)]; 97 return 1; 98 } 99 *shapearr=0; 100 return 0; 101 }