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 LinearPyramid.cpp
28 : : * \brief LinearPyramid implementation
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "LinearPyramid.hpp"
34 : : #include "MsqError.hpp"
35 : :
36 : : namespace MBMesquite
37 : : {
38 : :
39 : : static const char* nonlinear_error = "Attempt to use LinearTriangle mapping function for a nonlinear element\n";
40 : :
41 : 0 : static inline void set_equal_derivatives( double value, size_t* indices, MsqVector< 3 >* derivs, size_t& num_vtx )
42 : : {
43 : 0 : num_vtx = 5;
44 : 0 : indices[0] = 0;
45 : 0 : indices[1] = 1;
46 : 0 : indices[2] = 2;
47 : 0 : indices[3] = 3;
48 : 0 : indices[4] = 4;
49 : :
50 : 0 : derivs[0][0] = -value;
51 : 0 : derivs[0][1] = -value;
52 : 0 : derivs[0][2] = -0.25;
53 : :
54 : 0 : derivs[1][0] = value;
55 : 0 : derivs[1][1] = -value;
56 : 0 : derivs[1][2] = -0.25;
57 : :
58 : 0 : derivs[2][0] = value;
59 : 0 : derivs[2][1] = value;
60 : 0 : derivs[2][2] = -0.25;
61 : :
62 : 0 : derivs[3][0] = -value;
63 : 0 : derivs[3][1] = value;
64 : 0 : derivs[3][2] = -0.25;
65 : :
66 : 0 : derivs[4][0] = 0.0;
67 : 0 : derivs[4][1] = 0.0;
68 : 0 : derivs[4][2] = 1.0;
69 : 0 : }
70 : :
71 : 0 : static inline void set_edge_derivatives( unsigned base_corner, double value, size_t* indices, MsqVector< 3 >* derivs,
72 : : size_t& num_vtx )
73 : : {
74 : 0 : const int direction = base_corner % 2;
75 : 0 : const int edge_beg = base_corner;
76 : 0 : const int edge_end = ( base_corner + 1 ) % 4;
77 : 0 : const int adj_end = ( base_corner + 2 ) % 4;
78 : 0 : const int adj_beg = ( base_corner + 3 ) % 4;
79 : 0 : const int dir_sign = 2 * ( edge_beg / 2 ) - 1;
80 : 0 : const int oth_sign = 2 * ( ( edge_beg + 1 ) / 2 % 2 ) - 1;
81 : :
82 : 0 : num_vtx = 5;
83 : 0 : indices[0] = edge_beg;
84 : 0 : indices[1] = edge_end;
85 : 0 : indices[2] = adj_end;
86 : 0 : indices[3] = adj_beg;
87 : 0 : indices[4] = 4;
88 : :
89 : 0 : derivs[0][direction] = 2 * dir_sign * value;
90 : 0 : derivs[0][1 - direction] = oth_sign * value;
91 : 0 : derivs[0][2] = -0.5;
92 : :
93 : 0 : derivs[1][direction] = -2 * dir_sign * value;
94 : 0 : derivs[1][1 - direction] = oth_sign * value;
95 : 0 : derivs[1][2] = -0.5;
96 : :
97 : 0 : derivs[2][direction] = 0.0;
98 : 0 : derivs[2][1 - direction] = -oth_sign * value;
99 : 0 : derivs[2][2] = 0.0;
100 : :
101 : 0 : derivs[3][direction] = 0.0;
102 : 0 : derivs[3][1 - direction] = -oth_sign * value;
103 : 0 : derivs[3][2] = 0.0;
104 : :
105 : 0 : derivs[4][0] = 0.0;
106 : 0 : derivs[4][1] = 0.0;
107 : 0 : derivs[4][2] = 1.0;
108 : 0 : }
109 : :
110 : 2280 : static inline void set_corner_derivatives( unsigned corner, double value, size_t* indices, MsqVector< 3 >* derivs,
111 : : size_t& num_vtx )
112 : : {
113 : 2280 : const unsigned adj_in_xi = ( 5 - corner ) % 4;
114 : 2280 : const unsigned adj_in_eta = 3 - corner;
115 : :
116 : 2280 : const int dxi_sign = 2 * ( ( corner + 1 ) / 2 % 2 ) - 1;
117 : 2280 : const int deta_sign = 2 * ( corner / 2 ) - 1;
118 : 2280 : const double dxi_value = dxi_sign * value;
119 : 2280 : const double deta_value = deta_sign * value;
120 : :
121 : 2280 : num_vtx = 4;
122 : 2280 : indices[0] = corner;
123 : 2280 : indices[1] = adj_in_xi;
124 : 2280 : indices[2] = adj_in_eta;
125 : 2280 : indices[3] = 4;
126 : :
127 : 2280 : derivs[0][0] = dxi_value;
128 : 2280 : derivs[0][1] = deta_value;
129 : 2280 : derivs[0][2] = -1.0;
130 : :
131 : 2280 : derivs[1][0] = -dxi_value;
132 : 2280 : derivs[1][1] = 0.0;
133 : 2280 : derivs[1][2] = 0.0;
134 : :
135 : 2280 : derivs[2][0] = 0.0;
136 : 2280 : derivs[2][1] = -deta_value;
137 : 2280 : derivs[2][2] = 0.0;
138 : :
139 : 2280 : derivs[3][0] = 0.0;
140 : 2280 : derivs[3][1] = 0.0;
141 : 2280 : derivs[3][2] = 1.0;
142 : 2280 : }
143 : :
144 : 2280 : EntityTopology LinearPyramid::element_topology() const
145 : : {
146 : 2280 : return PYRAMID;
147 : : }
148 : :
149 : 2280 : int LinearPyramid::num_nodes() const
150 : : {
151 : 2280 : return 5;
152 : : }
153 : :
154 : 570 : NodeSet LinearPyramid::sample_points( NodeSet ) const
155 : : {
156 [ + - ]: 570 : NodeSet result;
157 [ + - ]: 570 : result.set_all_corner_nodes( PYRAMID );
158 [ + - ]: 570 : result.clear_corner_node( 4 );
159 : 570 : return result;
160 : : }
161 : :
162 : 0 : static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff )
163 : : {
164 : 0 : num_coeff = 1;
165 : 0 : indices_out[0] = corner;
166 : 0 : coeff_out[0] = 1.0;
167 : 0 : }
168 : :
169 : 0 : static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff )
170 : : {
171 : 0 : num_coeff = 2;
172 : 0 : coeff_out[0] = 0.5;
173 : 0 : coeff_out[1] = 0.5;
174 : :
175 [ # # ]: 0 : if( edge < 4 )
176 : : {
177 : 0 : indices_out[0] = edge;
178 : 0 : indices_out[1] = ( edge + 1 ) % 4;
179 : : }
180 : : else
181 : : {
182 : 0 : indices_out[0] = edge - 4;
183 : 0 : indices_out[1] = 4;
184 : : }
185 : 0 : }
186 : :
187 : 0 : static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff )
188 : : {
189 [ # # ]: 0 : if( face == 4 )
190 : : {
191 : 0 : num_coeff = 4;
192 : 0 : coeff_out[0] = 0.25;
193 : 0 : coeff_out[1] = 0.25;
194 : 0 : coeff_out[2] = 0.25;
195 : 0 : coeff_out[3] = 0.25;
196 : 0 : indices_out[0] = 0;
197 : 0 : indices_out[1] = 1;
198 : 0 : indices_out[2] = 2;
199 : 0 : indices_out[3] = 3;
200 : : }
201 : : else
202 : : {
203 : 0 : num_coeff = 3;
204 : 0 : indices_out[0] = face;
205 : 0 : indices_out[1] = ( face + 1 ) % 4;
206 : 0 : indices_out[2] = 4;
207 : 0 : coeff_out[0] = 0.25;
208 : 0 : coeff_out[1] = 0.25;
209 : 0 : coeff_out[2] = 0.50;
210 : : }
211 : 0 : }
212 : :
213 : 0 : static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff )
214 : : {
215 : 0 : num_coeff = 5;
216 : 0 : coeff_out[0] = 0.125;
217 : 0 : coeff_out[1] = 0.125;
218 : 0 : coeff_out[2] = 0.125;
219 : 0 : coeff_out[3] = 0.125;
220 : 0 : coeff_out[4] = 0.500;
221 : 0 : indices_out[0] = 0;
222 : 0 : indices_out[1] = 1;
223 : 0 : indices_out[2] = 2;
224 : 0 : indices_out[3] = 3;
225 : 0 : indices_out[4] = 4;
226 : 0 : }
227 : :
228 : 0 : void LinearPyramid::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out,
229 : : size_t& num_coeff, MsqError& err ) const
230 : : {
231 [ # # ]: 0 : if( nodeset.have_any_mid_node() )
232 : : {
233 [ # # ]: 0 : MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
234 : 0 : return;
235 : : }
236 : :
237 [ # # # # : 0 : switch( loc.dimension )
# ]
238 : : {
239 : : case 0:
240 : 0 : coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
241 : 0 : break;
242 : : case 1:
243 : 0 : coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
244 : 0 : break;
245 : : case 2:
246 : 0 : coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
247 : 0 : break;
248 : : case 3:
249 : 0 : coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
250 : 0 : break;
251 : : default:
252 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
253 : : }
254 : : }
255 : :
256 : 2280 : void LinearPyramid::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
257 : : MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
258 : : {
259 [ - + ]: 2280 : if( nodeset.have_any_mid_node() )
260 : : {
261 [ # # ]: 0 : MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
262 : 2280 : return;
263 : : }
264 : :
265 [ + - - - : 2280 : switch( loc.dimension )
- ]
266 : : {
267 : : case 0:
268 [ - + ]: 2280 : if( loc.number == 4 ) { set_equal_derivatives( 0.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
269 : : else
270 : : {
271 : 2280 : set_corner_derivatives( loc.number, 1.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
272 : : }
273 : 2280 : break;
274 : : case 1:
275 [ # # ]: 0 : if( loc.number < 4 )
276 : 0 : { set_edge_derivatives( loc.number, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
277 : : else
278 : : {
279 : 0 : set_corner_derivatives( loc.number - 4, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
280 : : }
281 : 0 : break;
282 : : case 2:
283 [ # # ]: 0 : if( loc.number == 4 ) { set_equal_derivatives( 0.5, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
284 : : else
285 : : {
286 : 0 : set_edge_derivatives( loc.number, 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
287 : : }
288 : 0 : break;
289 : : case 3:
290 : 0 : set_equal_derivatives( 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
291 : 0 : break;
292 : : default:
293 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
294 : : }
295 : : }
296 : :
297 : 0 : void LinearPyramid::ideal( Sample location, MsqMatrix< 3, 3 >& J, MsqError& ) const
298 : : {
299 : : // For an ideal element with unit edge length at the base and unit
300 : : // height, the Jacobian matrix is:
301 : : // | 1-zeta 0 1/2 - xi |
302 : : // | 0 1-zeta 1/2 - eta |
303 : : // | 0 0 1 |
304 : : //
305 : : // The coefficient to produce a unit determinant
306 : : // is therefore (1-zeta)^(-2/3).
307 : : //
308 : : // Thus the unit-determinate ideal element Jacobian
309 : : // is, given alpha = (1-zeta)^(-1/3):
310 : : //
311 : : // | 1/alpha 0 alpha^2 (1/2 - xi) |
312 : : // | 0 1/alpha alpha^2 (1/2 - eta)|
313 : : // | 0 0 alpha^2 |
314 : : //
315 : : // There are only three zeta values of interest:
316 : : // zeta = 1 : the degenerate case
317 : : // zeta = 0 : both 1/alpha and alpha^2 are 1.0
318 : : // zeta = 1/2 : 1/alpha = 1/cbrt(2.0) and alpha^2 = 2*(1/alpha)
319 : :
320 : : // special case for apex
321 [ # # ][ # # ]: 0 : if( location.dimension == 0 && location.number == 4 )
322 : : {
323 : 0 : J = MsqMatrix< 3, 3 >( 0.0 );
324 : 0 : return;
325 : : }
326 : :
327 : : // These are always zero
328 : 0 : J( 0, 1 ) = J( 1, 0 ) = J( 2, 0 ) = J( 2, 1 ) = 0.0;
329 : :
330 : : // Set diagonal terms and magnitude of terms in 3rd column based on zeta
331 : :
332 : : // All of the zeta=0 locations
333 : : double f;
334 [ # # ][ # # ]: 0 : if( location.dimension == 0 || ( location.dimension == 1 && location.number < 4 ) ||
[ # # ][ # # ]
335 [ # # ]: 0 : ( location.dimension == 2 && location.number == 4 ) )
336 : : {
337 : 0 : J( 0, 0 ) = J( 1, 1 ) = J( 2, 2 ) = 1.0;
338 : 0 : f = 0.5;
339 : : }
340 : : // all of the zeta=1/2 locations
341 : : else
342 : : {
343 : 0 : f = J( 0, 0 ) = J( 1, 1 ) = 0.79370052598409979;
344 : 0 : J( 2, 2 ) = 2.0 * f;
345 : : }
346 : :
347 : : // Set terms in 3rd column based on xi,eta
348 : :
349 : : // The xi = eta = 0.5 locations (mid-element in xi and eta)
350 [ # # ][ # # ]: 0 : if( location.dimension == 3 || ( location.dimension == 2 && location.number == 4 ) )
[ # # ]
351 : 0 : { J( 0, 2 ) = J( 1, 2 ) = 0.0; }
352 : : // The corner locations
353 [ # # ][ # # ]: 0 : else if( location.dimension == 0 || ( location.dimension == 1 && location.number >= 4 ) )
[ # # ]
354 : : {
355 [ # # # # : 0 : switch( location.number % 4 )
# ]
356 : : {
357 : : case 0:
358 : 0 : J( 0, 2 ) = f;
359 : 0 : J( 1, 2 ) = f;
360 : 0 : break;
361 : : case 1:
362 : 0 : J( 0, 2 ) = -f;
363 : 0 : J( 1, 2 ) = f;
364 : 0 : break;
365 : : case 2:
366 : 0 : J( 0, 2 ) = -f;
367 : 0 : J( 1, 2 ) = -f;
368 : 0 : break;
369 : : case 3:
370 : 0 : J( 0, 2 ) = f;
371 : 0 : J( 1, 2 ) = -f;
372 : 0 : break;
373 : : }
374 : : }
375 : : // The mid-edge locations
376 : : else
377 : : {
378 [ # # # # : 0 : switch( location.number )
# ]
379 : : {
380 : : case 0:
381 : 0 : J( 0, 2 ) = 0;
382 : 0 : J( 1, 2 ) = f;
383 : 0 : break;
384 : : case 1:
385 : 0 : J( 0, 2 ) = -f;
386 : 0 : J( 1, 2 ) = 0;
387 : 0 : break;
388 : : case 2:
389 : 0 : J( 0, 2 ) = 0;
390 : 0 : J( 1, 2 ) = -f;
391 : 0 : break;
392 : : case 3:
393 : 0 : J( 0, 2 ) = f;
394 : 0 : J( 1, 2 ) = 0;
395 : 0 : break;
396 : : }
397 : : }
398 : : }
399 : :
400 : : } // namespace MBMesquite
|
