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