/[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 3224 - (hide annotations)
Wed Sep 29 05:19:37 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3780 byte(s)
indent -linux -nce -i4 -bl -bli0 -l120

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     #include <stdlib.h>
15    
16 jfenwick 3221 /* Joel Fenwick - derived from info in Finley's Quadrature and shape files
17 jfenwick 3196
18 jfenwick 3221 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 jfenwick 3224 bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
24 jfenwick 3196 {
25     #define _dudley_s_alpha 0.58541019662496852
26     #define _dudley_s_beta 0.1381966011250105
27    
28     // {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints
29 jfenwick 3224 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} // Tet 4 points
39     };
40 jfenwick 3196
41     #undef _dudley_s_alpha
42     #undef _dudley_s_beta
43    
44 jfenwick 3224 static double **arr = 0;
45 jfenwick 3196
46 jfenwick 3224 if (arr == 0)
47     {
48     arr = malloc(8 * sizeof(double *)); // point occupies two slots to make things simpler
49     arr[0] = malloc(1 * sizeof(double));
50     arr[0][0] = 1.; // point
51     arr[1] = arr[0];
52     arr[2] = malloc(4 * sizeof(double)); // Line Single
53     arr[3] = malloc(4 * sizeof(double)); // Line 2
54 jfenwick 3196
55 jfenwick 3224 for (int i = 0; i < 2; ++i)
56     {
57     arr[2][2 * i] = 1 - _dudley_V[0][i];
58     arr[3][2 * i] = 1 - _dudley_V[1][i];
59 jfenwick 3198
60 jfenwick 3224 arr[2][2 * i + 1] = _dudley_V[0][i];
61     arr[3][2 * i + 1] = _dudley_V[1][i];
62     }
63 jfenwick 3196
64 jfenwick 3224 arr[4] = malloc(3 * sizeof(double)); // Tri single
65     arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
66     arr[4][1] = _dudley_V[2][0];
67     arr[4][2] = _dudley_V[2][1];
68 jfenwick 3196
69 jfenwick 3224 arr[5] = malloc(9 * sizeof(double)); // Tri 3
70     for (int i = 0; i < 3; ++i)
71     {
72     arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
73     arr[5][3 * i + 1] = _dudley_V[3][2 * i];
74     arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
75     }
76     arr[6] = malloc(4 * sizeof(double)); // Tet single
77     arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
78     arr[6][1] = _dudley_V[4][0];
79     arr[6][2] = _dudley_V[4][1];
80     arr[6][3] = _dudley_V[4][2];
81 jfenwick 3198
82 jfenwick 3224 arr[7] = malloc(16 * sizeof(double)); // Tet 4
83     for (int i = 0; i < 4; ++i)
84     {
85     double x = _dudley_V[5][3 * i];
86     double y = _dudley_V[5][3 * i + 1];
87     double z = _dudley_V[5][3 * i + 2];
88     arr[7][4 * i] = 1 - x - y - z;
89     arr[7][4 * i + 1] = x;
90     arr[7][4 * i + 2] = y;
91     arr[7][4 * i + 3] = z;
92     }
93     } // end if
94 jfenwick 3196
95 jfenwick 3224 if ((dim > -1) && (dim < 4))
96 jfenwick 3196 {
97 jfenwick 3224 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
98 jfenwick 3196 return 1;
99     }
100 jfenwick 3224 *shapearr = 0;
101 jfenwick 3196 return 0;
102 jfenwick 3198 }
103 jfenwick 3210
104 jfenwick 3224 const char *getElementName(ElementTypeId id)
105 jfenwick 3210 {
106 jfenwick 3224 switch (id)
107     {
108     case Point1:
109     return "Point1";
110     case Line2:
111     return "Line2";
112     case Tri3:
113     return "Tri3";
114     case Tet4:
115     return "Tet4";
116     case Line2Face:
117     return "Line2Face";
118     case Tri3Face:
119     return "Tri3Face";
120     case Tet4Face:
121     return "Tet4Face";
122     default:
123 jfenwick 3210 return "noElement";
124 jfenwick 3224 }
125 jfenwick 3210 }

  ViewVC Help
Powered by ViewVC 1.1.26