/[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 3227 - (hide annotations)
Thu Sep 30 06:07:08 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3743 byte(s)
Pass1 or moving MPI stuff out of paso

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     // {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     {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 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 jfenwick 3227 arr = MEMALLOC(8,double*); // point occupies two slots to make things simpler
50     arr[0] = MEMALLOC(1,double);
51 jfenwick 3224 arr[0][0] = 1.; // point
52     arr[1] = arr[0];
53 jfenwick 3227 arr[2] = MEMALLOC(4,double); // Line Single
54     arr[3] = MEMALLOC(4,double); // Line 2
55 jfenwick 3196
56 jfenwick 3224 for (int i = 0; i < 2; ++i)
57     {
58     arr[2][2 * i] = 1 - _dudley_V[0][i];
59     arr[3][2 * i] = 1 - _dudley_V[1][i];
60 jfenwick 3198
61 jfenwick 3224 arr[2][2 * i + 1] = _dudley_V[0][i];
62     arr[3][2 * i + 1] = _dudley_V[1][i];
63     }
64 jfenwick 3196
65 jfenwick 3227 arr[4] = MEMALLOC(3,double); // Tri single
66 jfenwick 3224 arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
67     arr[4][1] = _dudley_V[2][0];
68     arr[4][2] = _dudley_V[2][1];
69 jfenwick 3196
70 jfenwick 3227 arr[5] = MEMALLOC(9,double); // Tri 3
71 jfenwick 3224 for (int i = 0; i < 3; ++i)
72     {
73     arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
74     arr[5][3 * i + 1] = _dudley_V[3][2 * i];
75     arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
76     }
77 jfenwick 3227 arr[6] = MEMALLOC(4, double); // Tet single
78 jfenwick 3224 arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
79     arr[6][1] = _dudley_V[4][0];
80     arr[6][2] = _dudley_V[4][1];
81     arr[6][3] = _dudley_V[4][2];
82 jfenwick 3198
83 jfenwick 3227 arr[7] = MEMALLOC(16,double); // Tet 4
84 jfenwick 3224 for (int i = 0; i < 4; ++i)
85     {
86     double x = _dudley_V[5][3 * i];
87     double y = _dudley_V[5][3 * i + 1];
88     double z = _dudley_V[5][3 * i + 2];
89     arr[7][4 * i] = 1 - x - y - z;
90     arr[7][4 * i + 1] = x;
91     arr[7][4 * i + 2] = y;
92     arr[7][4 * i + 3] = z;
93     }
94     } // end if
95 jfenwick 3196
96 jfenwick 3224 if ((dim > -1) && (dim < 4))
97 jfenwick 3196 {
98 jfenwick 3224 *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
99 jfenwick 3196 return 1;
100     }
101 jfenwick 3224 *shapearr = 0;
102 jfenwick 3196 return 0;
103 jfenwick 3198 }
104 jfenwick 3210
105 jfenwick 3224 const char *getElementName(ElementTypeId id)
106 jfenwick 3210 {
107 jfenwick 3224 switch (id)
108     {
109     case Point1:
110     return "Point1";
111     case Line2:
112     return "Line2";
113     case Tri3:
114     return "Tri3";
115     case Tet4:
116     return "Tet4";
117     case Line2Face:
118     return "Line2Face";
119     case Tri3Face:
120     return "Tri3Face";
121     case Tet4Face:
122     return "Tet4Face";
123     default:
124 jfenwick 3210 return "noElement";
125 jfenwick 3224 }
126 jfenwick 3210 }

  ViewVC Help
Powered by ViewVC 1.1.26