/[escript]/branches/domexper/dudley/src/ShapeTable.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 <stdlib.h>
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

  ViewVC Help
Powered by ViewVC 1.1.26