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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3254 - (hide annotations)
Thu Oct 7 06:33:09 2010 UTC (8 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 3825 byte(s)
Back to C90 since M$ still won't support mixed code and declarations.
This should be a no-op.

1 jfenwick 3196 /*******************************************************
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 jfenwick 3227 #include "esysUtils/mem.h"
15 jfenwick 3196 #include <stdlib.h>
16    
17 jfenwick 3221 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files
18 jfenwick 3196
19 jfenwick 3221 This method is not threadsafe unless the initial call has completed
20     Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions
21     The dim argument is the dimension of the element not the dimension of the embedding space.
22     the reduced arg is whether the elements are reduced or not
23     */
24 jfenwick 3224 bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
25 jfenwick 3196 {
26     #define _dudley_s_alpha 0.58541019662496852
27     #define _dudley_s_beta 0.1381966011250105
28    
29 caltinay 3254 /* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */
30 jfenwick 3224 static const double _dudley_V[3 * 2][12] = {
31 caltinay 3254 {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 jfenwick 3224 {_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 caltinay 3254 _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha} /* Tet 4 points */
40 jfenwick 3224 };
41 jfenwick 3196
42     #undef _dudley_s_alpha
43     #undef _dudley_s_beta
44    
45 jfenwick 3224 static double **arr = 0;
46 jfenwick 3196
47 jfenwick 3224 if (arr == 0)
48     {
49 caltinay 3254 int i;
50     arr = MEMALLOC(8,double*); /* point occupies two slots to make things simpler */
51 jfenwick 3227 arr[0] = MEMALLOC(1,double);
52 caltinay 3254 arr[0][0] = 1.; /* point */
53 jfenwick 3224 arr[1] = arr[0];
54 caltinay 3254 arr[2] = MEMALLOC(4,double); /* Line Single */
55     arr[3] = MEMALLOC(4,double); /* Line 2 */
56 jfenwick 3196
57 caltinay 3254 for (i = 0; i < 2; ++i)
58 jfenwick 3224 {
59     arr[2][2 * i] = 1 - _dudley_V[0][i];
60     arr[3][2 * i] = 1 - _dudley_V[1][i];
61 jfenwick 3198
62 jfenwick 3224 arr[2][2 * i + 1] = _dudley_V[0][i];
63     arr[3][2 * i + 1] = _dudley_V[1][i];
64     }
65 jfenwick 3196
66 caltinay 3254 arr[4] = MEMALLOC(3,double); /* Tri single */
67 jfenwick 3224 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 jfenwick 3196
71 caltinay 3254 arr[5] = MEMALLOC(9,double); /* Tri 3 */
72     for (i = 0; i < 3; ++i)
73 jfenwick 3224 {
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 caltinay 3254 arr[6] = MEMALLOC(4, double); /* Tet single */
79 jfenwick 3224 arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
80     arr[6][1] = _dudley_V[4][0];
81     arr[6][2] = _dudley_V[4][1];
82     arr[6][3] = _dudley_V[4][2];
83 jfenwick 3198
84 caltinay 3254 arr[7] = MEMALLOC(16,double); /* Tet 4 */
85     for (i = 0; i < 4; ++i)
86 jfenwick 3224 {
87     double x = _dudley_V[5][3 * i];
88     double y = _dudley_V[5][3 * i + 1];
89     double z = _dudley_V[5][3 * i + 2];
90     arr[7][4 * i] = 1 - x - y - z;
91     arr[7][4 * i + 1] = x;
92     arr[7][4 * i + 2] = y;
93     arr[7][4 * i + 3] = z;
94     }
95 caltinay 3254 } /* end if */
96 jfenwick 3196
97 jfenwick 3224 if ((dim > -1) && (dim < 4))
98 jfenwick 3196 {
99 jfenwick 3224 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
100 jfenwick 3196 return 1;
101     }
102 jfenwick 3224 *shapearr = 0;
103 jfenwick 3196 return 0;
104 jfenwick 3198 }
105 jfenwick 3210
106 caltinay 3247 const char *getElementName(Dudley_ElementTypeId id)
107 jfenwick 3210 {
108 jfenwick 3224 switch (id)
109     {
110 caltinay 3247 case Dudley_Point1:
111 jfenwick 3224 return "Point1";
112 caltinay 3247 case Dudley_Line2:
113 jfenwick 3224 return "Line2";
114 caltinay 3247 case Dudley_Tri3:
115 jfenwick 3224 return "Tri3";
116 caltinay 3247 case Dudley_Tet4:
117 jfenwick 3224 return "Tet4";
118 caltinay 3247 case Dudley_Line2Face:
119 jfenwick 3224 return "Line2Face";
120 caltinay 3247 case Dudley_Tri3Face:
121 jfenwick 3224 return "Tri3Face";
122 caltinay 3247 case Dudley_Tet4Face:
123 jfenwick 3224 return "Tet4Face";
124     default:
125 jfenwick 3210 return "noElement";
126 jfenwick 3224 }
127 jfenwick 3210 }

  ViewVC Help
Powered by ViewVC 1.1.26