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 (2006) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27
28 /** \file LinearHexahedron.cpp
29 * \brief Implement shape function for linear hexahedron
30 * \author Jason Kraftcheck
31 */
32
33 #include "Mesquite.hpp"
34 #include "MsqError.hpp"
35 #include "LinearHexahedron.hpp"
36
37 namespace MBMesquite {
38
39 static const char* nonlinear_error
40 = "Attempt to use LinearHexahedron mapping function for a nonlinear element\n";
41
element_topology() const42 EntityTopology LinearHexahedron::element_topology() const
43 { return HEXAHEDRON; }
44
num_nodes() const45 int LinearHexahedron::num_nodes() const
46 { return 8; }
47
coeff_xi_sign(unsigned coeff)48 static inline int coeff_xi_sign( unsigned coeff )
49 { return 2*(((coeff+1)/2)%2) - 1; }
coeff_eta_sign(unsigned coeff)50 static inline int coeff_eta_sign( unsigned coeff )
51 { return 2*((coeff/2)%2) - 1; }
coeff_zeta_sign(unsigned coeff)52 static inline int coeff_zeta_sign( unsigned coeff )
53 { return 2*(coeff/4) - 1; }
54
coefficients_at_corner(unsigned corner,double * coeff_out,size_t * indices_out,size_t & num_coeff)55 static void coefficients_at_corner( unsigned corner,
56 double* coeff_out,
57 size_t* indices_out,
58 size_t& num_coeff )
59 {
60 num_coeff = 1;
61 coeff_out[0] = 1.0;
62 indices_out[0] = corner;
63 }
64
65 const unsigned xi = 0, eta = 1, zeta = 2;
66 const int edge_dir[] = { xi, eta, xi, eta, zeta, zeta, zeta, zeta, xi, eta, xi, eta };
67 const int edge_beg[] = { 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7 }; // start vertex by edge number
68 const int edge_end[] = { 1, 2, 3, 0, 4, 5, 6, 7, 5, 6, 7, 4 }; // end vetex by edge number
69 const int edge_opposite[] = {10,11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1 }; // opposite edge
70 const int edge_beg_orth1[] = { 3, 5, 1, 7, 1, 0, 3, 2, 7, 1, 5, 3 }; // vtx adjacent to edge start in edge_dir[e]+1 direction
71 const int edge_beg_orth2[] = { 4, 0, 6, 2, 3, 2, 1, 0, 0, 4, 2, 6 }; // vtx adjacent to edge start in edge_dir[e]+2 direction
72 const int edge_end_orth1[] = { 2, 6, 0, 4, 5, 4, 7, 6, 6, 2, 4, 0 }; // vtx adjacent to edge end in edge_dir[e]+1 direction
73 const int edge_end_orth2[] = { 5, 3, 7, 1, 7, 6, 5, 4, 1, 7, 3, 5 }; // vtx adjacent to edge end in edge_dir[e]+2 direction
74
coefficients_at_mid_edge(unsigned edge,double * coeff_out,size_t * indices_out,size_t & num_coeff)75 static void coefficients_at_mid_edge( unsigned edge,
76 double* coeff_out,
77 size_t* indices_out,
78 size_t& num_coeff )
79 {
80 num_coeff = 2;
81 coeff_out[0] = 0.5;
82 coeff_out[1] = 0.5;
83 indices_out[0] = edge_beg[edge];
84 indices_out[1] = edge_end[edge];
85 }
86
87
88 const int face_vtx[6][4] = { { 0, 1, 4, 5 }, // face 0 vertices
89 { 1, 2, 5, 6 }, // face 1 vertices
90 { 2, 3, 6, 7 }, // face 2
91 { 0, 3, 4, 7 }, // face 3
92 { 0, 1, 2, 3 }, // face 4
93 { 4, 5, 6, 7 } };// face 5
94 const int face_opp[6] = { 2, 3, 0, 1, 5, 4 }; // opposite faces on hex
95 const int face_dir[6] = { eta, xi, eta, xi, zeta, zeta }; // normal direction
96
coefficients_at_mid_face(unsigned face,double * coeff_out,size_t * indices_out,size_t & num_coeff)97 static void coefficients_at_mid_face( unsigned face,
98 double* coeff_out,
99 size_t* indices_out,
100 size_t& num_coeff )
101 {
102 num_coeff = 4;
103 coeff_out[0] = 0.25;
104 coeff_out[1] = 0.25;
105 coeff_out[2] = 0.25;
106 coeff_out[3] = 0.25;
107 indices_out[0] = face_vtx[face][0];
108 indices_out[1] = face_vtx[face][1];
109 indices_out[2] = face_vtx[face][2];
110 indices_out[3] = face_vtx[face][3];
111 }
112
113
114
coefficients_at_mid_elem(double * coeff_out,size_t * indices_out,size_t & num_coeff)115 static void coefficients_at_mid_elem( double* coeff_out,
116 size_t* indices_out,
117 size_t& num_coeff )
118 {
119 num_coeff = 8;
120 coeff_out[0] = 0.125;
121 coeff_out[1] = 0.125;
122 coeff_out[2] = 0.125;
123 coeff_out[3] = 0.125;
124 coeff_out[4] = 0.125;
125 coeff_out[5] = 0.125;
126 coeff_out[6] = 0.125;
127 coeff_out[7] = 0.125;
128 indices_out[0] = 0;
129 indices_out[1] = 1;
130 indices_out[2] = 2;
131 indices_out[3] = 3;
132 indices_out[4] = 4;
133 indices_out[5] = 5;
134 indices_out[6] = 6;
135 indices_out[7] = 7;
136 }
137
138
coefficients(Sample loc,NodeSet nodeset,double * coeff_out,size_t * indices_out,size_t & num_coeff,MsqError & err) const139 void LinearHexahedron::coefficients( Sample loc,
140 NodeSet nodeset,
141 double* coeff_out,
142 size_t* indices_out,
143 size_t& num_coeff,
144 MsqError& err ) const
145 {
146 if (nodeset.have_any_mid_node()) {
147 MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
148 return;
149 }
150
151 switch (loc.dimension) {
152 case 0:
153 coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
154 break;
155 case 1:
156 coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
157 break;
158 case 2:
159 coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
160 break;
161 case 3:
162 coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
163 break;
164 default:
165 MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
166 }
167 }
168
169
derivatives_at_corner(unsigned corner,size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx)170 static void derivatives_at_corner( unsigned corner,
171 size_t* vertex_indices_out,
172 MsqVector<3>* d_coeff_d_xi_out,
173 size_t& num_vtx )
174 {
175 const int xi_sign = coeff_xi_sign(corner);
176 const int eta_sign = coeff_eta_sign(corner);
177 const int zeta_sign = coeff_zeta_sign(corner);
178 const int offset = 4*(corner/4);
179 //const int adj_in_xi = offset + (9 - corner)%4;
180 const int adj_in_xi = corner + 1 - 2*(corner%2);
181 const int adj_in_eta = 3 + 2*offset - (int)corner;
182 const int adj_in_zeta = corner - zeta_sign*4;
183
184 num_vtx = 4;
185 vertex_indices_out[0] = corner;
186 vertex_indices_out[1] = adj_in_xi;
187 vertex_indices_out[2] = adj_in_eta;
188 vertex_indices_out[3] = adj_in_zeta;
189
190 d_coeff_d_xi_out[0][0] = xi_sign;
191 d_coeff_d_xi_out[0][1] = eta_sign;
192 d_coeff_d_xi_out[0][2] = zeta_sign;
193
194 d_coeff_d_xi_out[1][0] = -xi_sign ;
195 d_coeff_d_xi_out[1][1] = 0.0;
196 d_coeff_d_xi_out[1][2] = 0.0;
197
198 d_coeff_d_xi_out[2][0] = 0.0;
199 d_coeff_d_xi_out[2][1] = -eta_sign;
200 d_coeff_d_xi_out[2][2] = 0.0;
201
202 d_coeff_d_xi_out[3][0] = 0.0;
203 d_coeff_d_xi_out[3][1] = 0.0;
204 d_coeff_d_xi_out[3][2] =-zeta_sign;
205 }
206
207
derivatives_at_mid_edge(unsigned edge,size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx)208 static void derivatives_at_mid_edge( unsigned edge,
209 size_t* vertex_indices_out,
210 MsqVector<3>* d_coeff_d_xi_out,
211 size_t& num_vtx )
212 {
213 const int direction = edge_dir[edge];
214 const int ortho1 = (direction+1)%3;
215 const int ortho2 = (direction+2)%3;
216 const int vtx = edge_beg[edge];
217 const int signs[] = { coeff_xi_sign(vtx), coeff_eta_sign(vtx), coeff_zeta_sign(vtx) };
218 const int sign_dir = signs[direction];
219 const int sign_or1 = signs[ortho1 ];
220 const int sign_or2 = signs[ortho2 ];
221
222 num_vtx = 6;
223 vertex_indices_out[0] = edge_beg[edge];
224 vertex_indices_out[1] = edge_end[edge];
225 vertex_indices_out[2] = edge_beg_orth1[edge];
226 vertex_indices_out[3] = edge_end_orth1[edge];
227 vertex_indices_out[4] = edge_beg_orth2[edge];
228 vertex_indices_out[5] = edge_end_orth2[edge];
229
230 d_coeff_d_xi_out[0][direction] = sign_dir;
231 d_coeff_d_xi_out[0][ortho1 ] = sign_or1 * 0.5;
232 d_coeff_d_xi_out[0][ortho2 ] = sign_or2 * 0.5;
233
234 d_coeff_d_xi_out[1][direction] = -sign_dir;
235 d_coeff_d_xi_out[1][ortho1 ] = sign_or1 * 0.5;
236 d_coeff_d_xi_out[1][ortho2 ] = sign_or2 * 0.5;
237
238 d_coeff_d_xi_out[2][direction] = 0.0;
239 d_coeff_d_xi_out[2][ortho1 ] = -sign_or1 * 0.5;
240 d_coeff_d_xi_out[2][ortho2 ] = 0.0;
241
242 d_coeff_d_xi_out[3][direction] = 0.0;
243 d_coeff_d_xi_out[3][ortho1 ] = -sign_or1 * 0.5;
244 d_coeff_d_xi_out[3][ortho2 ] = 0.0;
245
246 d_coeff_d_xi_out[4][direction] = 0.0;
247 d_coeff_d_xi_out[4][ortho1 ] = 0.0;
248 d_coeff_d_xi_out[4][ortho2 ] = -sign_or2 * 0.5;
249
250 d_coeff_d_xi_out[5][direction] = 0.0;
251 d_coeff_d_xi_out[5][ortho1 ] = 0.0;
252 d_coeff_d_xi_out[5][ortho2 ] = -sign_or2 * 0.5;
253 }
254
255
derivatives_at_mid_face(unsigned face,size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx)256 static void derivatives_at_mid_face( unsigned face,
257 size_t* vertex_indices_out,
258 MsqVector<3>* d_coeff_d_xi_out,
259 size_t& num_vtx )
260 {
261 const int vtx_signs[4][3] = { { coeff_xi_sign (face_vtx[face][0]),
262 coeff_eta_sign (face_vtx[face][0]),
263 coeff_zeta_sign(face_vtx[face][0]) },
264 { coeff_xi_sign (face_vtx[face][1]),
265 coeff_eta_sign (face_vtx[face][1]),
266 coeff_zeta_sign(face_vtx[face][1]) },
267 { coeff_xi_sign (face_vtx[face][2]),
268 coeff_eta_sign (face_vtx[face][2]),
269 coeff_zeta_sign(face_vtx[face][2]) },
270 { coeff_xi_sign (face_vtx[face][3]),
271 coeff_eta_sign (face_vtx[face][3]),
272 coeff_zeta_sign(face_vtx[face][3]) } };
273 const int ortho_dir = face_dir[face];
274 const int face_dir1 = (ortho_dir+1) % 3;
275 const int face_dir2 = (ortho_dir+2) % 3;
276 const int ortho_sign = vtx_signs[0][ortho_dir]; // same for all 4 verts
277
278 num_vtx = 8;
279 vertex_indices_out[0] = face_vtx[face][0];
280 vertex_indices_out[1] = face_vtx[face][1];
281 vertex_indices_out[2] = face_vtx[face][2];
282 vertex_indices_out[3] = face_vtx[face][3];
283 vertex_indices_out[4] = face_vtx[face_opp[face]][0];
284 vertex_indices_out[5] = face_vtx[face_opp[face]][1];
285 vertex_indices_out[6] = face_vtx[face_opp[face]][2];
286 vertex_indices_out[7] = face_vtx[face_opp[face]][3];
287
288 d_coeff_d_xi_out[0][ortho_dir] = ortho_sign * 0.25;
289 d_coeff_d_xi_out[0][face_dir1] = vtx_signs[0][face_dir1] * 0.5;
290 d_coeff_d_xi_out[0][face_dir2] = vtx_signs[0][face_dir2] * 0.5;
291
292 d_coeff_d_xi_out[1][ortho_dir] = ortho_sign * 0.25;
293 d_coeff_d_xi_out[1][face_dir1] = vtx_signs[1][face_dir1] * 0.5;
294 d_coeff_d_xi_out[1][face_dir2] = vtx_signs[1][face_dir2] * 0.5;
295
296 d_coeff_d_xi_out[2][ortho_dir] = ortho_sign * 0.25;
297 d_coeff_d_xi_out[2][face_dir1] = vtx_signs[2][face_dir1] * 0.5;
298 d_coeff_d_xi_out[2][face_dir2] = vtx_signs[2][face_dir2] * 0.5;
299
300 d_coeff_d_xi_out[3][ortho_dir] = ortho_sign * 0.25;
301 d_coeff_d_xi_out[3][face_dir1] = vtx_signs[3][face_dir1] * 0.5;
302 d_coeff_d_xi_out[3][face_dir2] = vtx_signs[3][face_dir2] * 0.5;
303
304 d_coeff_d_xi_out[4][ortho_dir] = -ortho_sign * 0.25;
305 d_coeff_d_xi_out[4][face_dir1] = 0.0;
306 d_coeff_d_xi_out[4][face_dir2] = 0.0;
307
308 d_coeff_d_xi_out[5][ortho_dir] = -ortho_sign * 0.25;
309 d_coeff_d_xi_out[5][face_dir1] = 0.0;
310 d_coeff_d_xi_out[5][face_dir2] = 0.0;
311
312 d_coeff_d_xi_out[6][ortho_dir] = -ortho_sign * 0.25;
313 d_coeff_d_xi_out[6][face_dir1] = 0.0;
314 d_coeff_d_xi_out[6][face_dir2] = 0.0;
315
316 d_coeff_d_xi_out[7][ortho_dir] = -ortho_sign * 0.25;
317 d_coeff_d_xi_out[7][face_dir1] = 0.0;
318 d_coeff_d_xi_out[7][face_dir2] = 0.0;
319 }
320
derivatives_at_mid_elem(size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx)321 static void derivatives_at_mid_elem( size_t* vertex_indices_out,
322 MsqVector<3>* d_coeff_d_xi_out,
323 size_t& num_vtx )
324
325 {
326 num_vtx = 8;
327 vertex_indices_out[0] = 0;
328 vertex_indices_out[1] = 1;
329 vertex_indices_out[2] = 2;
330 vertex_indices_out[3] = 3;
331 vertex_indices_out[4] = 4;
332 vertex_indices_out[5] = 5;
333 vertex_indices_out[6] = 6;
334 vertex_indices_out[7] = 7;
335
336 d_coeff_d_xi_out[0][0] = -0.25;
337 d_coeff_d_xi_out[0][1] = -0.25;
338 d_coeff_d_xi_out[0][2] = -0.25;
339
340 d_coeff_d_xi_out[1][0] = 0.25;
341 d_coeff_d_xi_out[1][1] = -0.25;
342 d_coeff_d_xi_out[1][2] = -0.25;
343
344 d_coeff_d_xi_out[2][0] = 0.25;
345 d_coeff_d_xi_out[2][1] = 0.25;
346 d_coeff_d_xi_out[2][2] = -0.25;
347
348 d_coeff_d_xi_out[3][0] = -0.25;
349 d_coeff_d_xi_out[3][1] = 0.25;
350 d_coeff_d_xi_out[3][2] = -0.25;
351
352 d_coeff_d_xi_out[4][0] = -0.25;
353 d_coeff_d_xi_out[4][1] = -0.25;
354 d_coeff_d_xi_out[4][2] = 0.25;
355
356 d_coeff_d_xi_out[5][0] = 0.25;
357 d_coeff_d_xi_out[5][1] = -0.25;
358 d_coeff_d_xi_out[5][2] = 0.25;
359
360 d_coeff_d_xi_out[6][0] = 0.25;
361 d_coeff_d_xi_out[6][1] = 0.25;
362 d_coeff_d_xi_out[6][2] = 0.25;
363
364 d_coeff_d_xi_out[7][0] = -0.25;
365 d_coeff_d_xi_out[7][1] = 0.25;
366 d_coeff_d_xi_out[7][2] = 0.25;
367 }
368
369
derivatives(Sample loc,NodeSet nodeset,size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx,MsqError & err) const370 void LinearHexahedron::derivatives( Sample loc,
371 NodeSet nodeset,
372 size_t* vertex_indices_out,
373 MsqVector<3>* d_coeff_d_xi_out,
374 size_t& num_vtx,
375 MsqError& err ) const
376 {
377 if (nodeset.have_any_mid_node()) {
378 MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
379 return;
380 }
381
382 switch (loc.dimension) {
383 case 0:
384 derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
385 break;
386 case 1:
387 derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
388 break;
389 case 2:
390 derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
391 break;
392 case 3:
393 derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
394 break;
395 default:
396 MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
397 }
398 }
399
ideal(Sample,MsqMatrix<3,3> & J,MsqError &) const400 void LinearHexahedron::ideal( Sample ,
401 MsqMatrix<3,3>& J,
402 MsqError& ) const
403 {
404 J(0,0) = J(1,1) = J(2,2) = 1.0;
405 J(1,0) = J(0,1) = J(0,2) = 0.0;
406 J(2,0) = J(2,1) = J(1,2) = 0.0;
407 }
408
409 } // namespace MBMesquite
410