1 /* *****************************************************************
2     MESQUITE -- The Mesh Quality Improvement Toolkit
3 
4     Copyright 2006 Sandia National Laboratories.  Developed at the
5     University of Wisconsin--Madison under SNL contract number
6     624796.  The U.S. Government and the University of Wisconsin
7     retain certain rights to this software.
8 
9     This library is free software; you can redistribute it and/or
10     modify it under the terms of the GNU Lesser General Public
11     License as published by the Free Software Foundation; either
12     version 2.1 of the License, or (at your option) any later version.
13 
14     This library is distributed in the hope that it will be useful,
15     but WITHOUT ANY WARRANTY; without even the implied warranty of
16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17     Lesser General Public License for more details.
18 
19     You should have received a copy of the GNU Lesser General Public License
20     (lgpl.txt) along with this library; if not, write to the Free Software
21     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22 
23     (2006) kraftche@cae.wisc.edu
24 
25   ***************************************************************** */
26 
27 
28 /** \file IdealElements.cpp
29  *  \brief
30  *  \author Jason Kraftcheck
31  */
32 
33 #include "Mesquite.hpp"
34 #include "IdealElements.hpp"
35 #include "Vector3D.hpp"
36 
37 namespace MBMesquite {
38 
39 static Vector3D unit_quad[4] = { Vector3D( -0.5, -0.5, 0.0 ),
40                                  Vector3D(  0.5, -0.5, 0.0 ),
41                                  Vector3D(  0.5,  0.5, 0.0 ),
42                                  Vector3D( -0.5,  0.5, 0.0 ) };
43 
44 static Vector3D unit_hex[8] = { Vector3D(  0.5, -0.5, -0.5 ),
45                                 Vector3D(  0.5,  0.5, -0.5 ),
46                                 Vector3D( -0.5,  0.5, -0.5 ),
47                                 Vector3D( -0.5, -0.5, -0.5 ),
48                                 Vector3D(  0.5, -0.5,  0.5 ),
49                                 Vector3D(  0.5,  0.5,  0.5 ),
50                                 Vector3D( -0.5,  0.5,  0.5 ),
51                                 Vector3D( -0.5, -0.5,  0.5 ) };
52 
53 static Vector3D unit_edge_tri[3];
54 static Vector3D unit_edge_tet[4];
55 static Vector3D unit_edge_pyr[5];
56 static Vector3D unit_edge_wdg[6];
57 static Vector3D unit_height_pyr[5];
58 
59 static Vector3D unit_tri[3];
60 static Vector3D unit_tet[4];
61 static Vector3D unit_pyr[5];
62 static Vector3D unit_wdg[6];
63 static Vector3D unit_hex_pyr[5];
64 
65 static void init_tri( Vector3D* coords, double side );
66 static void init_tet( Vector3D* coords, double side );
67 static void init_pyr( Vector3D* coords, double side );
68 static void init_wdg( Vector3D* coords, double side );
69 static void init_hex_pyr( Vector3D* coords, double height );
70 
71 static const Vector3D* const* init_unit_edge( Vector3D** );
72 static const Vector3D* const* init_unit_elem( Vector3D** );
73 
unit_edge_element(EntityTopology ptype,bool punit_pyr)74 const Vector3D* unit_edge_element( EntityTopology ptype, bool punit_pyr )
75 {
76   static Vector3D* values[MIXED+1];
77   static const Vector3D* const* data = init_unit_edge( values );
78   return (ptype == PYRAMID && punit_pyr) ? data[MIXED] : data[ptype];
79 }
80 
unit_element(EntityTopology ptype,bool punit_pyr)81 const Vector3D* unit_element( EntityTopology ptype, bool punit_pyr )
82 {
83   static Vector3D* values[MIXED+1];
84   static const Vector3D* const* data = init_unit_elem( values );
85   return (ptype == PYRAMID && punit_pyr) ? data[MIXED] : data[ptype];
86 }
87 
init_unit_edge(Vector3D ** ptr)88 static const Vector3D* const* init_unit_edge( Vector3D** ptr )
89 {
90   for (unsigned i = 0; i < MIXED; ++i)
91     ptr[i] = 0;
92 
93   init_tri( unit_edge_tri, 1.0 );
94   init_tet( unit_edge_tet, 1.0 );
95   init_pyr( unit_edge_pyr, 1.0 );
96   init_wdg( unit_edge_wdg, 1.0 );
97   init_hex_pyr( unit_height_pyr, 1.0 );
98 
99   ptr[TRIANGLE] = unit_edge_tri;
100   ptr[QUADRILATERAL] = unit_quad;
101   ptr[TETRAHEDRON] = unit_edge_tet;
102   ptr[PYRAMID] = unit_edge_pyr;
103   ptr[PRISM] = unit_edge_wdg;
104   ptr[HEXAHEDRON] = unit_hex;
105   ptr[MIXED] = unit_height_pyr;
106   return ptr;
107 }
108 
init_unit_elem(Vector3D ** ptr)109 static const Vector3D* const* init_unit_elem( Vector3D** ptr )
110 {
111   for (unsigned i = 0; i < MIXED; ++i)
112     ptr[i] = 0;
113 
114   init_tri( unit_tri, 2.0 * pow( 3.0, -0.25 ) );
115   init_tet( unit_tet, MBMesquite::cbrt( 3.0 ) * sqrt(2.0) );
116   init_pyr( unit_pyr, pow( 18.0, 1.0/6.0 ) );
117   init_wdg( unit_wdg, MBMesquite::cbrt( 4.0 ) * pow( 3.0, -1.0/6.0 ) );
118   init_hex_pyr( unit_hex_pyr, MBMesquite::cbrt( 3.0 ) );
119 
120   ptr[TRIANGLE] = unit_tri;
121   ptr[QUADRILATERAL] = unit_quad;
122   ptr[TETRAHEDRON] = unit_tet;
123   ptr[PYRAMID] = unit_pyr;
124   ptr[PRISM] = unit_wdg;
125   ptr[HEXAHEDRON] = unit_hex;
126   ptr[MIXED] = unit_hex_pyr;
127   return ptr;
128 }
129 
init_tri(Vector3D * coords,double side)130 static void init_tri( Vector3D* coords, double side )
131 {
132   const double third_height = side * sqrt(3.0) / 6.0;
133   coords[1] = Vector3D( -0.5*side,  -third_height, 0.0 );
134   coords[2] = Vector3D(  0.5*side,  -third_height, 0.0 );
135   coords[0] = Vector3D(  0.0,      2*third_height, 0.0 );
136 }
137 
init_tet(Vector3D * coords,double side)138 static void init_tet( Vector3D* coords, double side )
139 {
140   const double height = side * sqrt(2.0/3.0);
141   const double third_base = side * sqrt(3.0) / 6.0;
142   coords[0] = Vector3D( -0.5*side,  -third_base, -0.25*height );
143   coords[1] = Vector3D(  0.5*side,  -third_base, -0.25*height );
144   coords[2] = Vector3D(  0.0,      2*third_base, -0.25*height );
145   coords[3] = Vector3D(  0.0,      0.0,           0.75*height );
146 }
147 
init_pyr(Vector3D * coords,double side)148 static void init_pyr( Vector3D* coords, double side )
149 {
150   const double height = side * sqrt(2.0) * 0.5;
151   coords[0] = Vector3D(  0.5*side, -0.5*side, -0.25*height );
152   coords[1] = Vector3D(  0.5*side,  0.5*side, -0.25*height );
153   coords[2] = Vector3D( -0.5*side,  0.5*side, -0.25*height );
154   coords[3] = Vector3D( -0.5*side, -0.5*side, -0.25*height );
155   coords[4] = Vector3D(  0.0,       0.0,       0.75*height );
156 }
157 
init_wdg(Vector3D * coords,double side)158 static void init_wdg( Vector3D* coords, double side )
159 {
160   const double third_height = side * sqrt(3.0) / 6.0;
161   coords[0] = Vector3D( -0.5*side,  -third_height, -0.5*side );
162   coords[1] = Vector3D(  0.5*side,  -third_height, -0.5*side );
163   coords[2] = Vector3D(  0.0,      2*third_height, -0.5*side );
164   coords[3] = Vector3D( -0.5*side,  -third_height,  0.5*side );
165   coords[4] = Vector3D(  0.5*side,  -third_height,  0.5*side );
166   coords[5] = Vector3D(  0.0,      2*third_height,  0.5*side );
167 }
168 
init_hex_pyr(Vector3D * coords,double side)169 static void init_hex_pyr( Vector3D* coords, double side )
170 {
171   coords[0] = Vector3D(  0.5*side, -0.5*side, -0.25*side );
172   coords[1] = Vector3D(  0.5*side,  0.5*side, -0.25*side );
173   coords[2] = Vector3D( -0.5*side,  0.5*side, -0.25*side );
174   coords[3] = Vector3D( -0.5*side, -0.5*side, -0.25*side );
175   coords[4] = Vector3D(  0.0,       0.0,       0.75*side );
176 }
177 
178 } // namespace MBMesquite
179