/[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 3210 - (hide annotations)
Mon Sep 27 04:13:19 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3487 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.26