/[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 3210 - (show annotations)
Mon Sep 27 04:13:19 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3487 byte(s)


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 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 }
108
109
110 const char* getElementName(ElementTypeId id)
111 {
112 switch(id)
113 {
114 case Point1: return "Point1";
115 case Line2: return "Line2";
116 case Tri3: return "Tri3";
117 case Tet4: return "Tet4";
118 case Line2Face: return "Line2Face";
119 case Tri3Face: return "Tri3Face";
120 case Tet4Face: return "Tet4Face";
121 default:
122 return "noElement";
123 }
124 }
125
126
127

  ViewVC Help
Powered by ViewVC 1.1.26