← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/ShapeTable.c

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/ShapeTable.c
 
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 "esysUtils/mem.h"
 
15
#include <stdlib.h>
 
16
 
 
17
/* Joel Fenwick - derived from info in Finley's Quadrature and shape files
 
18
 
 
19
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
bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
 
25
{
 
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
    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
 
 
42
#undef _dudley_s_alpha
 
43
#undef _dudley_s_beta
 
44
 
 
45
    static double **arr = 0;
 
46
 
 
47
    if (arr == 0)
 
48
    {
 
49
        int i;
 
50
        arr = MEMALLOC(8,double*);      /* point occupies two slots to make things simpler */
 
51
        arr[0] = MEMALLOC(1,double);
 
52
        arr[0][0] = 1.;         /* point */
 
53
        arr[1] = arr[0];
 
54
        arr[2] = MEMALLOC(4,double);    /* Line Single */
 
55
        arr[3] = MEMALLOC(4,double);    /* Line 2 */
 
56
 
 
57
        for (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
        arr[4] = MEMALLOC(3,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
        arr[5] = MEMALLOC(9,double);    /* Tri 3 */
 
72
        for (i = 0; i < 3; ++i)
 
73
        {
 
74
            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
        arr[6] = MEMALLOC(4, double);   /* Tet single */
 
79
        arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
 
80
        arr[6][1] = _dudley_V[4][0];
 
81
        arr[6][2] = _dudley_V[4][1];
 
82
        arr[6][3] = _dudley_V[4][2];
 
83
 
 
84
        arr[7] = MEMALLOC(16,double);   /* Tet 4 */
 
85
        for (i = 0; i < 4; ++i)
 
86
        {
 
87
            double x = _dudley_V[5][3 * i];
 
88
            double y = _dudley_V[5][3 * i + 1];
 
89
            double z = _dudley_V[5][3 * i + 2];
 
90
            arr[7][4 * i] = 1 - x - y - z;
 
91
            arr[7][4 * i + 1] = x;
 
92
            arr[7][4 * i + 2] = y;
 
93
            arr[7][4 * i + 3] = z;
 
94
        }
 
95
    }                           /* end if */
 
96
 
 
97
    if ((dim > -1) && (dim < 4))
 
98
    {
 
99
        *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
 
100
        return 1;
 
101
    }
 
102
    *shapearr = 0;
 
103
    return 0;
 
104
}
 
105
 
 
106
const char *getElementName(Dudley_ElementTypeId id)
 
107
{
 
108
    switch (id)
 
109
    {
 
110
    case Dudley_Point1:
 
111
        return "Point1";
 
112
    case Dudley_Line2:
 
113
        return "Line2";
 
114
    case Dudley_Tri3:
 
115
        return "Tri3";
 
116
    case Dudley_Tet4:
 
117
        return "Tet4";
 
118
    case Dudley_Line2Face:
 
119
        return "Line2Face";
 
120
    case Dudley_Tri3Face:
 
121
        return "Tri3Face";
 
122
    case Dudley_Tet4Face:
 
123
        return "Tet4Face";
 
124
    default:
 
125
        return "noElement";
 
126
    }
 
127
}