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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3196 - (show annotations)
Wed Sep 22 01:18:52 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 2990 byte(s)
moving slowly
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 <stdlib.h>
15
16 // Joel Fenwick - derived from info in Finley's Quadrature and shape files
17
18 // 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 bool_t getQuadShape(dim_t dim, bool_t reduced, const double** shapearr)
23 {
24 #define _dudley_s_alpha 0.58541019662496852
25 #define _dudley_s_beta 0.1381966011250105
26
27
28 // {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints
29 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
39 } // 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 arr=malloc(8*sizeof(double*)); // point occupies two slots to make things simpler
50 arr[0]=malloc(1*sizeof(double));
51 arr[0][0]=1.; // point
52 arr[1]=arr[0];
53 arr[2]=malloc(4*sizeof(double)); // Line Single
54 arr[3]=malloc(4*sizeof(double)); // Line 2
55 arr[4]=malloc(4*sizeof(double)); // Tri single
56
57 for (int i=0;i<2;++i)
58 {
59 arr[2][2*i]=1-_dudley_V[0][i];
60 arr[3][2*i]=1-_dudley_V[1][i];
61 arr[4][2*i]=1-_dudley_V[2][i];
62
63 arr[2][2*i+1]=_dudley_V[0][i];
64 arr[3][2*i+1]=_dudley_V[1][i];
65 arr[4][2*i+1]=_dudley_V[2][i];
66 }
67
68 arr[5]=malloc(9*sizeof(double)); // Tri 3
69 for (int i=0;i<3;++i)
70 {
71 arr[5][3*i]=1-_dudley_V[3][3*i] -_dudley_V[2][3*i+1];
72 arr[5][3*i+1]=_dudley_V[3][3*i];
73 arr[5][3*i+2]=_dudley_V[3][3*i+1];
74 }
75 arr[6]=malloc(3*sizeof(double)); // Tet single
76 arr[6][0]=1-_dudley_V[4][0] -_dudley_V[2][1] -_dudley_V[2][2];
77 arr[6][1]=_dudley_V[4][0];
78 arr[6][2]= _dudley_V[4][1];
79 arr[6][3]=_dudley_V[4][2];
80
81 arr[7]=malloc(16*sizeof(double));
82 for (int i=0;i<4;++i)
83 {
84 double x=_dudley_V[5][0];
85 double y=_dudley_V[5][1];
86 double z=_dudley_V[5][2];
87 arr[7][4*i]=1-x-y-z;
88 arr[7][4*i+1]=x;
89 arr[7][4*i+2]=y;
90 arr[7][4*i+3]=z;
91 }
92 } // end if
93
94 if ((dim>-1) && (dim<4))
95 {
96 *shapearr=arr[(!reduced)?(2*dim+1):(2*dim)];
97 return 1;
98 }
99 *shapearr=0;
100 return 0;
101 }

  ViewVC Help
Powered by ViewVC 1.1.26