/[escript]/branches/doubleplusgood/dudley/src/ShapeTable.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/ShapeTable.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 4143 byte(s)
like sand though the hourglass
1 /*****************************************************************************
2 *
3 * Copyright (c) 2003-2013 by University of Queensland
4 * http://www.uq.edu.au
5 *
6 * Primary Business: Queensland, Australia
7 * Licensed under the Open Software License version 3.0
8 * http://www.opensource.org/licenses/osl-3.0.php
9 *
10 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 * Development since 2012 by School of Earth Sciences
12 *
13 *****************************************************************************/
14
15 #include "ShapeTable.h"
16 #include "esysUtils/mem.h"
17 #include <stdlib.h>
18
19 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files
20
21 This method is not threadsafe unless the initial call has completed
22 Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions
23 The dim argument is the dimension of the element not the dimension of the embedding space.
24 the reduced arg is whether the elements are reduced or not
25 */
26 bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
27 {
28 #define _dudley_s_alpha 0.58541019662496852
29 #define _dudley_s_beta 0.1381966011250105
30
31 /* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */
32 static const double _dudley_V[3 * 2][12] = {
33 {0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line single */
34 {(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line 2 points */
35 {1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tri single */
36 {0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0}, /* Tri 3 points */
37 {0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tet single */
38 {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta,
39 _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta,
40 _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta,
41 _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha} /* Tet 4 points */
42 };
43
44 #undef _dudley_s_alpha
45 #undef _dudley_s_beta
46
47 static double **arr = 0;
48
49 if (arr == 0)
50 {
51 int i;
52 arr = new double*[8]; /* point occupies two slots to make things simpler */
53 arr[0] = new double[1];
54 arr[0][0] = 1.; /* point */
55 arr[1] = arr[0];
56 arr[2] = new double[4]; /* Line Single */
57 arr[3] = new double[4]; /* Line 2 */
58
59 /*
60 for (i = 0; i < 2; ++i)
61 {
62 arr[2][2 * i] = 1 - _dudley_V[0][i];
63 arr[3][2 * i] = 1 - _dudley_V[1][i];
64
65 arr[2][2 * i + 1] = _dudley_V[0][i];
66 arr[3][2 * i + 1] = _dudley_V[1][i];
67 }
68 */
69
70 for (i = 0; i < 2; ++i)
71 {
72 arr[2][2 * i] = 1 - _dudley_V[0][i];
73 arr[2][2 * i + 1] = _dudley_V[0][i];
74 }
75 for (i = 0; i < 2; ++i)
76 {
77 arr[3][2 * i] = 1 - _dudley_V[1][i];
78 arr[3][2 * i + 1] = _dudley_V[1][i];
79 }
80
81
82
83 arr[4] = new double[3]; /* Tri single */
84 arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
85 arr[4][1] = _dudley_V[2][0];
86 arr[4][2] = _dudley_V[2][1];
87
88 arr[5] = new double[9]; /* Tri 3 */
89 for (i = 0; i < 3; ++i)
90 {
91 arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
92 arr[5][3 * i + 1] = _dudley_V[3][2 * i];
93 arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
94 }
95 arr[6] = new double[4]; /* Tet single */
96 arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
97 arr[6][1] = _dudley_V[4][0];
98 arr[6][2] = _dudley_V[4][1];
99 arr[6][3] = _dudley_V[4][2];
100
101 arr[7] = new double[16]; /* Tet 4 */
102 for (i = 0; i < 4; ++i)
103 {
104 double x = _dudley_V[5][3 * i];
105 double y = _dudley_V[5][3 * i + 1];
106 double z = _dudley_V[5][3 * i + 2];
107 arr[7][4 * i] = 1 - x - y - z;
108 arr[7][4 * i + 1] = x;
109 arr[7][4 * i + 2] = y;
110 arr[7][4 * i + 3] = z;
111 }
112 } /* end if */
113
114 if ((dim > -1) && (dim < 4))
115 {
116 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
117 return 1;
118 }
119 *shapearr = 0;
120 return 0;
121 }
122
123 const char *getElementName(Dudley_ElementTypeId id)
124 {
125 switch (id)
126 {
127 case Dudley_Point1:
128 return "Point1";
129 case Dudley_Line2:
130 return "Line2";
131 case Dudley_Tri3:
132 return "Tri3";
133 case Dudley_Tet4:
134 return "Tet4";
135 case Dudley_Line2Face:
136 return "Line2Face";
137 case Dudley_Tri3Face:
138 return "Tri3Face";
139 case Dudley_Tet4Face:
140 return "Tet4Face";
141 default:
142 return "noElement";
143 }
144 }

  ViewVC Help
Powered by ViewVC 1.1.26