/[escript]/trunk/dudley/src/ShapeTable.c
ViewVC logotype

Contents of /trunk/dudley/src/ShapeTable.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3264 - (show annotations)
Tue Oct 12 04:40:06 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 4064 byte(s)
Minor changes to avoid optimiser bug in intel-10

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 "esysUtils/mem.h"
15 #include <stdlib.h>
16
17 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files
18
19 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 bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
25 {
26 #define _dudley_s_alpha 0.58541019662496852
27 #define _dudley_s_beta 0.1381966011250105
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} /* 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 int i;
50 arr = MEMALLOC(8,double*); /* point occupies two slots to make things simpler */
51 arr[0] = MEMALLOC(1,double);
52 arr[0][0] = 1.; /* point */
53 arr[1] = arr[0];
54 arr[2] = MEMALLOC(4,double); /* Line Single */
55 arr[3] = MEMALLOC(4,double); /* Line 2 */
56
57 /*
58 for (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
68 for (i = 0; i < 2; ++i)
69 {
70 arr[2][2 * i] = 1 - _dudley_V[0][i];
71 arr[2][2 * i + 1] = _dudley_V[0][i];
72 }
73 for (i = 0; i < 2; ++i)
74 {
75 arr[3][2 * i] = 1 - _dudley_V[1][i];
76 arr[3][2 * i + 1] = _dudley_V[1][i];
77 }
78
79
80
81 arr[4] = MEMALLOC(3,double); /* Tri single */
82 arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
83 arr[4][1] = _dudley_V[2][0];
84 arr[4][2] = _dudley_V[2][1];
85
86 arr[5] = MEMALLOC(9,double); /* Tri 3 */
87 for (i = 0; i < 3; ++i)
88 {
89 arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
90 arr[5][3 * i + 1] = _dudley_V[3][2 * i];
91 arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
92 }
93 arr[6] = MEMALLOC(4, double); /* Tet single */
94 arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
95 arr[6][1] = _dudley_V[4][0];
96 arr[6][2] = _dudley_V[4][1];
97 arr[6][3] = _dudley_V[4][2];
98
99 arr[7] = MEMALLOC(16,double); /* Tet 4 */
100 for (i = 0; i < 4; ++i)
101 {
102 double x = _dudley_V[5][3 * i];
103 double y = _dudley_V[5][3 * i + 1];
104 double z = _dudley_V[5][3 * i + 2];
105 arr[7][4 * i] = 1 - x - y - z;
106 arr[7][4 * i + 1] = x;
107 arr[7][4 * i + 2] = y;
108 arr[7][4 * i + 3] = z;
109 }
110 } /* end if */
111
112 if ((dim > -1) && (dim < 4))
113 {
114 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
115 return 1;
116 }
117 *shapearr = 0;
118 return 0;
119 }
120
121 const char *getElementName(Dudley_ElementTypeId id)
122 {
123 switch (id)
124 {
125 case Dudley_Point1:
126 return "Point1";
127 case Dudley_Line2:
128 return "Line2";
129 case Dudley_Tri3:
130 return "Tri3";
131 case Dudley_Tet4:
132 return "Tet4";
133 case Dudley_Line2Face:
134 return "Line2Face";
135 case Dudley_Tri3Face:
136 return "Tri3Face";
137 case Dudley_Tet4Face:
138 return "Tet4Face";
139 default:
140 return "noElement";
141 }
142 }

  ViewVC Help
Powered by ViewVC 1.1.26