dudley/src/ShapeTable.cpp

/*****************************************************************************
*
* Copyright (c) 2003-2013 by University of Queensland
* http://www.uq.edu.au
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
* Development since 2012 by School of Earth Sciences
*
*****************************************************************************/

#include "ShapeTable.h"
#include "esysUtils/mem.h"
#include <stdlib.h>

/* Joel Fenwick - derived from info in Finley's Quadrature and shape files

This method is not threadsafe unless the initial call has completed
Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions
The dim argument is the dimension of the element not the dimension of the embedding space.
the reduced arg is whether the elements are reduced or not
*/
bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
{
#define _dudley_s_alpha 0.58541019662496852
#define _dudley_s_beta  0.1381966011250105

/* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */
    static const double _dudley_V[3 * 2][12] = {
    {0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line single */
    {(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line 2 points */
    {1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tri single */
    {0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0},   /* Tri 3 points */
    {0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0},  /* Tet single */
    {_dudley_s_beta, _dudley_s_beta, _dudley_s_beta,
     _dudley_s_alpha, _dudley_s_beta, _dudley_s_beta,
     _dudley_s_beta, _dudley_s_alpha, _dudley_s_beta,
     _dudley_s_beta, _dudley_s_beta, _dudley_s_alpha}   /* Tet 4 points */
    };

#undef _dudley_s_alpha
#undef _dudley_s_beta

    static double **arr = 0;

    if (arr == 0)
    {
    int i;
    arr = new double*[8];   /* point occupies two slots to make things simpler */
    arr[0] = new double[1];
    arr[0][0] = 1.;     /* point */
    arr[1] = arr[0];
    arr[2] = new double[4]; /* Line Single */
    arr[3] = new double[4]; /* Line 2 */

/*
    for (i = 0; i < 2; ++i)
    {
        arr[2][2 * i] = 1 - _dudley_V[0][i];
        arr[3][2 * i] = 1 - _dudley_V[1][i];

        arr[2][2 * i + 1] = _dudley_V[0][i];
        arr[3][2 * i + 1] = _dudley_V[1][i];
    }
*/

    for (i = 0; i < 2; ++i)
    {
        arr[2][2 * i] = 1 - _dudley_V[0][i];
        arr[2][2 * i + 1] = _dudley_V[0][i];
    }
    for (i = 0; i < 2; ++i)
    {
        arr[3][2 * i] = 1 - _dudley_V[1][i];
        arr[3][2 * i + 1] = _dudley_V[1][i];
    }


    arr[4] = new double[3]; /* Tri single */
    arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
    arr[4][1] = _dudley_V[2][0];
    arr[4][2] = _dudley_V[2][1];

    arr[5] = new double[9]; /* Tri 3 */
    for (i = 0; i < 3; ++i)
    {
        arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
        arr[5][3 * i + 1] = _dudley_V[3][2 * i];
        arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
    }
    arr[6] = new  double[4];    /* Tet single */
    arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
    arr[6][1] = _dudley_V[4][0];
    arr[6][2] = _dudley_V[4][1];
    arr[6][3] = _dudley_V[4][2];

    arr[7] = new double[16];    /* Tet 4 */
    for (i = 0; i < 4; ++i)
    {
        double x = _dudley_V[5][3 * i];
        double y = _dudley_V[5][3 * i + 1];
        double z = _dudley_V[5][3 * i + 2];
        arr[7][4 * i] = 1 - x - y - z;
        arr[7][4 * i + 1] = x;
        arr[7][4 * i + 2] = y;
        arr[7][4 * i + 3] = z;
    }
    }               /* end if */

    if ((dim > -1) && (dim < 4))
    {
    *shapearr = arr[(!reduced) ? (2 * dim + 1) : (2 * dim)];
    return 1;
    }
    *shapearr = 0;
    return 0;
}

const char *getElementName(Dudley_ElementTypeId id)
{
    switch (id)
    {
    case Dudley_Point1:
    return "Point1";
    case Dudley_Line2:
    return "Line2";
    case Dudley_Tri3:
    return "Tri3";
    case Dudley_Tet4:
    return "Tet4";
    case Dudley_Line2Face:
    return "Line2Face";
    case Dudley_Tri3Face:
    return "Tri3Face";
    case Dudley_Tet4Face:
    return "Tet4Face";
    default:
    return "noElement";
    }
}
1	/*****************************************************************************
2	*
3	* Copyright (c) 2003-2013 by University of Queensland
4	* http://www.uq.edu.au
5	*
6	* Primary Business: Queensland, Australia
7	* Licensed under the Open Software License version 3.0
8	* http://www.opensource.org/licenses/osl-3.0.php
9	*
10	* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11	* Development since 2012 by School of Earth Sciences
12	*
13	*****************************************************************************/
14
15	#include "ShapeTable.h"
16	#include "esysUtils/mem.h"
17	#include <stdlib.h>
18
19	/* Joel Fenwick - derived from info in Finley's Quadrature and shape files
20
21	This method is not threadsafe unless the initial call has completed
22	Evaluates the shape functions at nodes (This is the S value from the finley ShapeFunctions
23	The dim argument is the dimension of the element not the dimension of the embedding space.
24	the reduced arg is whether the elements are reduced or not
25	*/
26	bool_t getQuadShape(dim_t dim, bool_t reduced, const double **shapearr)
27	{
28	#define _dudley_s_alpha 0.58541019662496852
29	#define _dudley_s_beta 0.1381966011250105
30
31	/* {Line, TRI, TET} X {single_quad_point, more} X max number of quadpoints */
32	static const double _dudley_V[3 * 2][12] = {
33	{0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line single */
34	{(1. - .577350269189626) / 2., (1. + .577350269189626) / 2., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Line 2 points */
35	{1 / 3., 1 / 3., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tri single */
36	{0.5, 0, 0, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0}, /* Tri 3 points */
37	{0.25, 0.25, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0}, /* Tet single */
38	{_dudley_s_beta, _dudley_s_beta, _dudley_s_beta,
39	_dudley_s_alpha, _dudley_s_beta, _dudley_s_beta,
40	_dudley_s_beta, _dudley_s_alpha, _dudley_s_beta,
41	_dudley_s_beta, _dudley_s_beta, _dudley_s_alpha} /* Tet 4 points */
42	};
43
44	#undef _dudley_s_alpha
45	#undef _dudley_s_beta
46
47	static double **arr = 0;
48
49	if (arr == 0)
50	{
51	int i;
52	arr = new double[8]; / point occupies two slots to make things simpler */
53	arr[0] = new double[1];
54	arr[0][0] = 1.; /* point */
55	arr[1] = arr[0];
56	arr[2] = new double[4]; /* Line Single */
57	arr[3] = new double[4]; /* Line 2 */
58
59	/*
60	for (i = 0; i < 2; ++i)
61	{
62	arr[2][2 * i] = 1 - _dudley_V[0][i];
63	arr[3][2 * i] = 1 - _dudley_V[1][i];
64
65	arr[2][2 * i + 1] = _dudley_V[0][i];
66	arr[3][2 * i + 1] = _dudley_V[1][i];
67	}
68	*/
69
70	for (i = 0; i < 2; ++i)
71	{
72	arr[2][2 * i] = 1 - _dudley_V[0][i];
73	arr[2][2 * i + 1] = _dudley_V[0][i];
74	}
75	for (i = 0; i < 2; ++i)
76	{
77	arr[3][2 * i] = 1 - _dudley_V[1][i];
78	arr[3][2 * i + 1] = _dudley_V[1][i];
79	}
80
81
82
83	arr[4] = new double[3]; /* Tri single */
84	arr[4][0] = 1. - _dudley_V[2][0] - _dudley_V[2][1];
85	arr[4][1] = _dudley_V[2][0];
86	arr[4][2] = _dudley_V[2][1];
87
88	arr[5] = new double[9]; /* Tri 3 */
89	for (i = 0; i < 3; ++i)
90	{
91	arr[5][3 * i] = 1 - _dudley_V[3][2 * i] - _dudley_V[3][2 * i + 1];
92	arr[5][3 * i + 1] = _dudley_V[3][2 * i];
93	arr[5][3 * i + 2] = _dudley_V[3][2 * i + 1];
94	}
95	arr[6] = new double[4]; /* Tet single */
96	arr[6][0] = 1 - _dudley_V[4][0] - _dudley_V[4][1] - _dudley_V[4][2];
97	arr[6][1] = _dudley_V[4][0];
98	arr[6][2] = _dudley_V[4][1];
99	arr[6][3] = _dudley_V[4][2];
100
101	arr[7] = new double[16]; /* Tet 4 */
102	for (i = 0; i < 4; ++i)
103	{
104	double x = _dudley_V[5][3 * i];
105	double y = _dudley_V[5][3 * i + 1];
106	double z = _dudley_V[5][3 * i + 2];
107	arr[7][4 * i] = 1 - x - y - z;
108	arr[7][4 * i + 1] = x;
109	arr[7][4 * i + 2] = y;
110	arr[7][4 * i + 3] = z;
111	}
112	} /* end if */
113
114	if ((dim > -1) && (dim < 4))
115	{
116	shapearr = arr[(!reduced) ? (2 dim + 1) : (2 * dim)];
117	return 1;
118	}
119	*shapearr = 0;
120	return 0;
121	}
122
123	const char *getElementName(Dudley_ElementTypeId id)
124	{
125	switch (id)
126	{
127	case Dudley_Point1:
128	return "Point1";
129	case Dudley_Line2:
130	return "Line2";
131	case Dudley_Tri3:
132	return "Tri3";
133	case Dudley_Tet4:
134	return "Tet4";
135	case Dudley_Line2Face:
136	return "Line2Face";
137	case Dudley_Tri3Face:
138	return "Tri3Face";
139	case Dudley_Tet4Face:
140	return "Tet4Face";
141	default:
142	return "noElement";
143	}
144	}