1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2006 Lawrence Livermore National Laboratory. Under
5 the terms of Contract B545069 with the University of Wisconsin --
6 Madison, Lawrence Livermore National Laboratory retains certain
7 rights in 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 kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26 /** \file HexLagrangeShape.cpp
27 * \author Nicholas Voshell
28 */
29
30 #include "HexLagrangeShape.hpp"
31 #include "MsqError.hpp"
32
33 namespace MBMesquite {
34
element_topology() const35 EntityTopology HexLagrangeShape::element_topology() const
36 { return HEXAHEDRON; }
37
num_nodes() const38 int HexLagrangeShape::num_nodes() const
39 { return 27; }
40
coefficients(Sample loc,NodeSet nodeset,double * coeff_out,size_t * indices_out,size_t & num_coeff,MsqError & err) const41 void HexLagrangeShape::coefficients( Sample loc,
42 NodeSet nodeset,
43 double* coeff_out,
44 size_t* indices_out,
45 size_t& num_coeff,
46 MsqError& err ) const
47 {
48 if (nodeset.num_nodes() != 27) {
49 MSQ_SETERR(err)("Mapping function supports only 27-node hexahedra with no slaved nodes.",
50 MsqError::UNSUPPORTED_ELEMENT);
51 return;
52 }
53
54 switch (loc.dimension) {
55 case 0: //Corner sample point - assume that it is there always
56 num_coeff = 1;
57 indices_out[0] = loc.number;
58 coeff_out[0] = 1.0;
59 break;
60 case 1: //Line sample point - check if it is there,
61 num_coeff = 1;
62 indices_out[0] = loc.number+8;
63 coeff_out[0] = 1.0;
64 break;
65 case 2: //Face sample point - check if it is there,
66 num_coeff = 1;
67 indices_out[0] = loc.number+20;
68 coeff_out[0] = 1.0;
69 break;
70 case 3: //Center sample point - check if it is there,
71 num_coeff = 1;
72 indices_out[0] = 26;
73 coeff_out[0] = 1.0;
74 break;
75 default:
76 MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT,
77 "Request for dimension %d mapping function value"
78 "for a quadratic hexahedral element", loc.dimension);
79 }
80
81 }
82
83
derivatives(Sample loc,NodeSet nodeset,size_t * vertices,MsqVector<3> * derivs,size_t & num_vtx,MsqError & err) const84 void HexLagrangeShape::derivatives( Sample loc,
85 NodeSet nodeset,
86 size_t* vertices,
87 MsqVector<3>* derivs,
88 size_t& num_vtx,
89 MsqError& err ) const
90 {
91 if (nodeset.num_nodes() != 27) {
92 MSQ_SETERR(err)("Mapping function supports only 27-node hexahedra with no slaved nodes.",
93 MsqError::UNSUPPORTED_ELEMENT);
94 return;
95 }
96
97 //r coordinate
98 const int HLS1[]= {1, 3, 3, 1, 1, 3, 3, 1, 2, 3, 2, 1, 1, 3, 3, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 2, 2};
99 //s coordinate
100 const int HLS2[]= {1, 1, 3, 3, 1, 1, 3, 3, 1, 2, 3, 2, 1, 1, 3, 3, 1, 2, 3, 2, 1, 2, 3, 2, 2, 2, 2};
101 //t coordinate
102 const int HLS3[]= {1, 1, 1, 1, 3, 3, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 2, 2, 2, 2, 1, 3, 2};
103
104 const int HLSup1[]= {12, 13, 14, 15, -1, -1, -1, -1, 20, 21, 22, 23, 4, 5, 6, 7, -1, -1, -1, -1, 16, 17, 18, 19, 26, -1, 25};
105 const int HLSup2[]= {4, 5, 6, 7, -1, -1, -1, -1, 16, 17, 18, 19, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 25, -1, -1};
106 const int HLSdn1[]= {-1, -1, -1, -1, 12, 13, 14, 15, -1, -1, -1, -1, 0, 1, 2, 3, 20, 21, 22, 23, 8, 9, 10, 11, -1, 26, 24};
107 const int HLSdn2[]= {-1, -1, -1, -1, 0, 1, 2, 3, -1, -1, -1, -1, -1, -1, -1, -1, 8, 9, 10, 11, -1, -1, -1, -1, -1, 24, -1};
108
109 const int HLSlt1[]= {-1, 8, 10, -1, -1, 16, 18, -1, 0, 24, 3, -1, -1, 20, 22, -1, 4, 25, 7, -1, 12, 26, 15, -1, 11, 19, 23};
110 const int HLSlt2[]= {-1, 0, 3, -1, -1, 4, 7, -1, -1, 11, -1, -1, -1, 12, 15, -1, -1, 19, -1, -1, -1, 23, -1, -1, -1, -1, -1};
111 const int HLSrt1[]= {8, -1, -1, 10, 16, -1, -1, 18, 1, -1, 2, 24, 20, -1, -1, 22, 5, -1, 6, 25, 13, -1, 14, 26, 9, 17, 21};
112 const int HLSrt2[]= {1, -1, -1, 2, 5, -1, -1, 6, -1, -1, -1, 9, 13, -1, -1, 14, -1, -1, -1, 17, -1, -1, -1, 21, -1, -1, -1};
113
114 const int HLSft1[]= {-1, -1, 9, 11, -1, -1, 17, 19, -1, 1, 24, 0, -1, -1, 21, 23, -1, 5, 25, 4, -1, 13, 26, 12, 8, 16, 20};
115 const int HLSft2[]= {-1, -1, 1, 0, -1, -1, 5, 4, -1, -1, 8, -1, -1, -1, 13, 12, -1, -1, 16, -1, -1, -1, 20, -1, -1, -1, -1};
116 const int HLSbk1[]= {11, 9, -1, -1, 19, 17, -1, -1, 24, 2, -1, 3, 23, 21, -1, -1, 25, 6, -1, 7, 26, 14, -1, 15, 10, 18, 22};
117 const int HLSbk2[]= {3, 2, -1, -1, 7, 6, -1, -1, 10, -1, -1, -1, 15, 14, -1, -1, 18, -1, -1, -1, 22, -1, -1, -1, -1, -1, -1};
118
119 double entries[] = {0, -3, 0, 3};
120
121 int location=loc.number;
122
123 switch (loc.dimension) {
124 case 1:
125 location+=8;
126 break;
127 case 2:
128 location+=20;
129 break;
130 case 3:
131 location+=26;
132 break;
133 }
134
135 num_vtx=0;
136
137 if (location != 26) {
138 vertices[num_vtx]=location;
139 derivs[num_vtx][0]=entries[HLS1[location]];
140 derivs[num_vtx][1]=entries[HLS2[location]];
141 derivs[num_vtx][2]=entries[HLS3[location]];
142 num_vtx++;
143 }
144
145 if(HLSup1[location]!=-1){
146 vertices[num_vtx]=HLSup1[location];
147 if(HLS3[location]==2) {
148 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=1;
149 } else {
150 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=4;
151 }
152 num_vtx++;
153 }
154
155 if(HLSdn1[location]!=-1){
156 vertices[num_vtx]=HLSdn1[location];
157 if(HLS3[location]==2) {
158 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=-1;
159 } else {
160 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=-4;
161 }
162 num_vtx++;
163 }
164
165 if(HLSup2[location]!=-1){
166 vertices[num_vtx]=HLSup2[location];
167 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=-1;
168 num_vtx++;
169 }
170
171 if(HLSdn2[location]!=-1){
172 vertices[num_vtx]=HLSdn2[location];
173 derivs[num_vtx][0]=0; derivs[num_vtx][1]=0; derivs[num_vtx][2]=1;
174 num_vtx++;
175 }
176
177 if(HLSlt1[location]!=-1){
178 vertices[num_vtx]=HLSlt1[location];
179 if(HLS1[location]==2) {
180 derivs[num_vtx][0]=-1; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
181 } else {
182 derivs[num_vtx][0]=-4; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
183 }
184 num_vtx++;
185 }
186
187 if(HLSrt1[location]!=-1){
188 vertices[num_vtx]=HLSrt1[location];
189 if(HLS1[location]==2) {
190 derivs[num_vtx][0]=1; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
191 } else {
192 derivs[num_vtx][0]=4; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
193 }
194 num_vtx++;
195 }
196
197 if(HLSlt2[location]!=-1){
198 vertices[num_vtx]=HLSlt2[location];
199 derivs[num_vtx][0]=1; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
200 num_vtx++;
201 }
202
203 if(HLSrt2[location]!=-1){
204 vertices[num_vtx]=HLSrt2[location];
205 derivs[num_vtx][0]=-1; derivs[num_vtx][1]=0; derivs[num_vtx][2]=0;
206 num_vtx++;
207 }
208
209 if(HLSft1[location]!=-1){
210 vertices[num_vtx]=HLSft1[location];
211 if(HLS2[location]==2) {
212 derivs[num_vtx][0]=0; derivs[num_vtx][1]=-1; derivs[num_vtx][2]=0;
213 } else {
214 derivs[num_vtx][0]=0; derivs[num_vtx][1]=-4; derivs[num_vtx][2]=0;
215 }
216 num_vtx++;
217 }
218
219 if(HLSbk1[location]!=-1){
220 vertices[num_vtx]=HLSbk1[location];
221 if(HLS2[location]==2) {
222 derivs[num_vtx][0]=0; derivs[num_vtx][1]=1; derivs[num_vtx][2]=0;
223 } else {
224 derivs[num_vtx][0]=0; derivs[num_vtx][1]=4; derivs[num_vtx][2]=0;
225 }
226 num_vtx++;
227 }
228
229 if(HLSft2[location]!=-1){
230 vertices[num_vtx]=HLSft2[location];
231 derivs[num_vtx][0]=0; derivs[num_vtx][1]=1; derivs[num_vtx][2]=0;
232 num_vtx++;
233 }
234
235 if(HLSbk2[location]!=-1){
236 vertices[num_vtx]=HLSbk2[location];
237 derivs[num_vtx][0]=0; derivs[num_vtx][1]=-1; derivs[num_vtx][2]=0;
238 num_vtx++;
239 }
240
241 }
242
243
ideal(Sample,MsqMatrix<3,3> & J,MsqError &) const244 void HexLagrangeShape::ideal( Sample ,
245 MsqMatrix<3,3>& J,
246 MsqError& ) const
247 {
248 J = MsqMatrix<3,3>(1.0);
249 }
250
251 } // namespace MBMesquite
252