/[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 3196 - (hide 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 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     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