Branch data Line data Source code
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) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TetLagrangeShape.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TetLagrangeShape.hpp"
34 : : #include "MsqError.hpp"
35 : : #include <assert.h>
36 : :
37 : : namespace MBMesquite
38 : : {
39 : :
40 : 4228329 : EntityTopology TetLagrangeShape::element_topology() const
41 : : {
42 : 4228329 : return TETRAHEDRON;
43 : : }
44 : :
45 : 4228329 : int TetLagrangeShape::num_nodes() const
46 : : {
47 : 4228329 : return 10;
48 : : }
49 : :
50 : 4488804 : NodeSet TetLagrangeShape::sample_points( NodeSet ns ) const
51 : : {
52 [ - + ]: 4488804 : if( ns.have_any_mid_node() ) { ns.set_all_corner_nodes( TETRAHEDRON ); }
53 : : else
54 : : {
55 : 4488804 : ns.clear();
56 : 4488804 : ns.set_mid_region_node();
57 : : }
58 : 4488804 : return ns;
59 : : }
60 : :
61 : 0 : static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff )
62 : : {
63 : 0 : num_coeff = 1;
64 : 0 : indices_out[0] = corner;
65 : 0 : coeff_out[0] = 1.0;
66 : 0 : }
67 : :
68 : 0 : static void coefficients_at_mid_edge( unsigned edge, NodeSet nodeset, double* coeff_out, size_t* indices_out,
69 : : size_t& num_coeff )
70 : : {
71 [ # # ]: 0 : if( nodeset.mid_edge_node( edge ) )
72 : : { // if mid-edge node is present
73 : 0 : num_coeff = 1;
74 : 0 : indices_out[0] = 4 + edge;
75 : 0 : coeff_out[0] = 1.0;
76 : : }
77 : : else
78 : : { // no mid node on edge
79 : 0 : num_coeff = 2;
80 : 0 : coeff_out[0] = coeff_out[1] = 0.5;
81 [ # # ]: 0 : if( edge < 3 )
82 : : {
83 : 0 : indices_out[0] = edge;
84 : 0 : indices_out[1] = ( edge + 1 ) % 3;
85 : : }
86 : : else
87 : : {
88 : 0 : indices_out[0] = edge - 3;
89 : 0 : indices_out[1] = 3;
90 : : }
91 : : }
92 : 0 : }
93 : :
94 : 0 : static void coefficients_at_mid_face( unsigned face, NodeSet nodeset, double* coeff_out, size_t* indices_out,
95 : : size_t& num_coeff )
96 : : {
97 : 0 : const double one_ninth = 1.0 / 9.0;
98 : 0 : const double two_ninth = 2.0 / 9.0;
99 : 0 : const double four_ninth = 4.0 / 9.0;
100 : :
101 [ # # ]: 0 : if( face < 3 )
102 : : {
103 : 0 : const int next = ( face + 1 ) % 3;
104 : 0 : indices_out[0] = face;
105 : 0 : indices_out[1] = next;
106 : 0 : indices_out[2] = 3;
107 : 0 : coeff_out[0] = -one_ninth;
108 : 0 : coeff_out[1] = -one_ninth;
109 : 0 : coeff_out[2] = -one_ninth;
110 : 0 : num_coeff = 3;
111 [ # # ]: 0 : if( nodeset.mid_edge_node( face ) )
112 : : {
113 : 0 : indices_out[num_coeff] = 4 + face;
114 : 0 : coeff_out[num_coeff] = four_ninth;
115 : 0 : ++num_coeff;
116 : : }
117 : : else
118 : : {
119 : 0 : coeff_out[0] += two_ninth;
120 : 0 : coeff_out[1] += two_ninth;
121 : : }
122 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 + next ) )
123 : : {
124 : 0 : indices_out[num_coeff] = 7 + next;
125 : 0 : coeff_out[num_coeff] = four_ninth;
126 : 0 : ++num_coeff;
127 : : }
128 : : else
129 : : {
130 : 0 : coeff_out[1] += two_ninth;
131 : 0 : coeff_out[2] += two_ninth;
132 : : }
133 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 + face ) )
134 : : {
135 : 0 : indices_out[num_coeff] = 7 + face;
136 : 0 : coeff_out[num_coeff] = four_ninth;
137 : 0 : ++num_coeff;
138 : : }
139 : : else
140 : : {
141 : 0 : coeff_out[0] += two_ninth;
142 : 0 : coeff_out[2] += two_ninth;
143 : : }
144 : : }
145 : : else
146 : : {
147 [ # # ]: 0 : assert( face == 3 );
148 : 0 : indices_out[0] = 0;
149 : 0 : indices_out[1] = 1;
150 : 0 : indices_out[2] = 2;
151 : 0 : coeff_out[0] = -one_ninth;
152 : 0 : coeff_out[1] = -one_ninth;
153 : 0 : coeff_out[2] = -one_ninth;
154 : 0 : num_coeff = 3;
155 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
156 : : {
157 : 0 : indices_out[num_coeff] = 4;
158 : 0 : coeff_out[num_coeff] = four_ninth;
159 : 0 : ++num_coeff;
160 : : }
161 : : else
162 : : {
163 : 0 : coeff_out[0] += two_ninth;
164 : 0 : coeff_out[1] += two_ninth;
165 : : }
166 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
167 : : {
168 : 0 : indices_out[num_coeff] = 5;
169 : 0 : coeff_out[num_coeff] = four_ninth;
170 : 0 : ++num_coeff;
171 : : }
172 : : else
173 : : {
174 : 0 : coeff_out[1] += two_ninth;
175 : 0 : coeff_out[2] += two_ninth;
176 : : }
177 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
178 : : {
179 : 0 : indices_out[num_coeff] = 6;
180 : 0 : coeff_out[num_coeff] = four_ninth;
181 : 0 : ++num_coeff;
182 : : }
183 : : else
184 : : {
185 : 0 : coeff_out[2] += two_ninth;
186 : 0 : coeff_out[0] += two_ninth;
187 : : }
188 : : }
189 : 0 : }
190 : :
191 : 0 : static void coefficients_at_mid_elem( NodeSet nodeset, double* coeff_out, size_t* indices_out, size_t& num_coeff )
192 : : {
193 : 0 : num_coeff = 4;
194 : 0 : indices_out[0] = 0;
195 : 0 : indices_out[1] = 1;
196 : 0 : indices_out[2] = 2;
197 : 0 : indices_out[3] = 3;
198 : 0 : coeff_out[0] = -0.125;
199 : 0 : coeff_out[1] = -0.125;
200 : 0 : coeff_out[2] = -0.125;
201 : 0 : coeff_out[3] = -0.125;
202 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
203 : : {
204 : 0 : indices_out[num_coeff] = 4;
205 : 0 : coeff_out[num_coeff] = 0.25;
206 : 0 : ++num_coeff;
207 : : }
208 : : else
209 : : {
210 : 0 : coeff_out[0] += 0.125;
211 : 0 : coeff_out[1] += 0.125;
212 : : }
213 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
214 : : {
215 : 0 : indices_out[num_coeff] = 5;
216 : 0 : coeff_out[num_coeff] = 0.25;
217 : 0 : ++num_coeff;
218 : : }
219 : : else
220 : : {
221 : 0 : coeff_out[1] += 0.125;
222 : 0 : coeff_out[2] += 0.125;
223 : : }
224 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
225 : : {
226 : 0 : indices_out[num_coeff] = 6;
227 : 0 : coeff_out[num_coeff] = 0.25;
228 : 0 : ++num_coeff;
229 : : }
230 : : else
231 : : {
232 : 0 : coeff_out[2] += 0.125;
233 : 0 : coeff_out[0] += 0.125;
234 : : }
235 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
236 : : {
237 : 0 : indices_out[num_coeff] = 7;
238 : 0 : coeff_out[num_coeff] = 0.25;
239 : 0 : ++num_coeff;
240 : : }
241 : : else
242 : : {
243 : 0 : coeff_out[0] += 0.125;
244 : 0 : coeff_out[3] += 0.125;
245 : : }
246 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
247 : : {
248 : 0 : indices_out[num_coeff] = 8;
249 : 0 : coeff_out[num_coeff] = 0.25;
250 : 0 : ++num_coeff;
251 : : }
252 : : else
253 : : {
254 : 0 : coeff_out[1] += 0.125;
255 : 0 : coeff_out[3] += 0.125;
256 : : }
257 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
258 : : {
259 : 0 : indices_out[num_coeff] = 9;
260 : 0 : coeff_out[num_coeff] = 0.25;
261 : 0 : ++num_coeff;
262 : : }
263 : : else
264 : : {
265 : 0 : coeff_out[2] += 0.125;
266 : 0 : coeff_out[3] += 0.125;
267 : : }
268 : 0 : }
269 : :
270 : 0 : void TetLagrangeShape::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out,
271 : : size_t& num_coeff, MsqError& err ) const
272 : : {
273 [ # # ]: 0 : if( nodeset.have_any_mid_face_node() | nodeset.have_any_mid_region_node() )
274 : : {
275 : : MSQ_SETERR( err )
276 [ # # ]: 0 : ( "TetLagrangeShape does not support mid-face/mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
277 : 0 : return;
278 : : }
279 : :
280 [ # # # # : 0 : switch( loc.dimension )
# ]
281 : : {
282 : : case 0:
283 : 0 : coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
284 : 0 : break;
285 : : case 1:
286 : 0 : coefficients_at_mid_edge( loc.number, nodeset, coeff_out, indices_out, num_coeff );
287 : 0 : break;
288 : : case 2:
289 : 0 : coefficients_at_mid_face( loc.number, nodeset, coeff_out, indices_out, num_coeff );
290 : 0 : break;
291 : : case 3:
292 : 0 : coefficients_at_mid_elem( nodeset, coeff_out, indices_out, num_coeff );
293 : 0 : break;
294 : : default:
295 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
296 : : }
297 : : }
298 : :
299 : 4565348 : static void get_linear_derivatives( size_t* vertices, MsqVector< 3 >* derivs )
300 : : {
301 : 4565348 : vertices[0] = 0;
302 : 4565348 : derivs[0][0] = -1.0;
303 : 4565348 : derivs[0][1] = -1.0;
304 : 4565348 : derivs[0][2] = -1.0;
305 : :
306 : 4565348 : vertices[1] = 1;
307 : 4565348 : derivs[1][0] = 1.0;
308 : 4565348 : derivs[1][1] = 0.0;
309 : 4565348 : derivs[1][2] = 0.0;
310 : :
311 : 4565348 : vertices[2] = 2;
312 : 4565348 : derivs[2][0] = 0.0;
313 : 4565348 : derivs[2][1] = 1.0;
314 : 4565348 : derivs[2][2] = 0.0;
315 : :
316 : 4565348 : vertices[3] = 3;
317 : 4565348 : derivs[3][0] = 0.0;
318 : 4565348 : derivs[3][1] = 0.0;
319 : 4565348 : derivs[3][2] = 1.0;
320 : 4565348 : }
321 : :
322 : : static const unsigned edges[][2] = { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 3 }, { 2, 3 } };
323 : :
324 : 0 : static void derivatives_at_corner( unsigned corner, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs,
325 : : size_t& num_vtx )
326 : : {
327 : : // begin with derivatives for linear tetrahedron
328 : 0 : num_vtx = 4;
329 : 0 : get_linear_derivatives( vertices, derivs );
330 : :
331 : : // adjust for the presence of mid-edge nodes
332 [ # # # # : 0 : switch( corner )
# ]
333 : : {
334 : : case 0:
335 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
336 : : {
337 : 0 : vertices[num_vtx] = 4;
338 : 0 : derivs[num_vtx][0] = 4.0;
339 : 0 : derivs[num_vtx][1] = 0.0;
340 : 0 : derivs[num_vtx][2] = 0.0;
341 : 0 : derivs[0][0] -= 2.0;
342 : 0 : derivs[1][0] -= 2.0;
343 : 0 : ++num_vtx;
344 : : }
345 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
346 : : {
347 : 0 : vertices[num_vtx] = 6;
348 : 0 : derivs[num_vtx][0] = 0.0;
349 : 0 : derivs[num_vtx][1] = 4.0;
350 : 0 : derivs[num_vtx][2] = 0.0;
351 : 0 : derivs[0][1] -= 2.0;
352 : 0 : derivs[2][1] -= 2.0;
353 : 0 : ++num_vtx;
354 : : }
355 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
356 : : {
357 : 0 : vertices[num_vtx] = 7;
358 : 0 : derivs[num_vtx][0] = 0.0;
359 : 0 : derivs[num_vtx][1] = 0.0;
360 : 0 : derivs[num_vtx][2] = 4.0;
361 : 0 : derivs[0][2] -= 2.0;
362 : 0 : derivs[3][2] -= 2.0;
363 : 0 : ++num_vtx;
364 : : }
365 : 0 : break;
366 : :
367 : : case 1:
368 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
369 : : {
370 : 0 : vertices[num_vtx] = 4;
371 : 0 : derivs[num_vtx][0] = -4.0;
372 : 0 : derivs[num_vtx][1] = -4.0;
373 : 0 : derivs[num_vtx][2] = -4.0;
374 : 0 : derivs[0][0] += 2.0;
375 : 0 : derivs[0][1] += 2.0;
376 : 0 : derivs[0][2] += 2.0;
377 : 0 : derivs[1][0] += 2.0;
378 : 0 : derivs[1][1] += 2.0;
379 : 0 : derivs[1][2] += 2.0;
380 : 0 : ++num_vtx;
381 : : }
382 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
383 : : {
384 : 0 : vertices[num_vtx] = 5;
385 : 0 : derivs[num_vtx][0] = 0.0;
386 : 0 : derivs[num_vtx][1] = 4.0;
387 : 0 : derivs[num_vtx][2] = 0.0;
388 : 0 : derivs[1][1] -= 2.0;
389 : 0 : derivs[2][1] -= 2.0;
390 : 0 : ++num_vtx;
391 : : }
392 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
393 : : {
394 : 0 : vertices[num_vtx] = 8;
395 : 0 : derivs[num_vtx][0] = 0.0;
396 : 0 : derivs[num_vtx][1] = 0.0;
397 : 0 : derivs[num_vtx][2] = 4.0;
398 : 0 : derivs[1][2] -= 2.0;
399 : 0 : derivs[3][2] -= 2.0;
400 : 0 : ++num_vtx;
401 : : }
402 : 0 : break;
403 : :
404 : : case 2:
405 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
406 : : {
407 : 0 : vertices[num_vtx] = 5;
408 : 0 : derivs[num_vtx][0] = 4.0;
409 : 0 : derivs[num_vtx][1] = 0.0;
410 : 0 : derivs[num_vtx][2] = 0.0;
411 : 0 : derivs[1][0] -= 2.0;
412 : 0 : derivs[2][0] -= 2.0;
413 : 0 : ++num_vtx;
414 : : }
415 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
416 : : {
417 : 0 : vertices[num_vtx] = 6;
418 : 0 : derivs[num_vtx][0] = -4.0;
419 : 0 : derivs[num_vtx][1] = -4.0;
420 : 0 : derivs[num_vtx][2] = -4.0;
421 : 0 : derivs[0][0] += 2.0;
422 : 0 : derivs[0][1] += 2.0;
423 : 0 : derivs[0][2] += 2.0;
424 : 0 : derivs[2][0] += 2.0;
425 : 0 : derivs[2][1] += 2.0;
426 : 0 : derivs[2][2] += 2.0;
427 : 0 : ++num_vtx;
428 : : }
429 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
430 : : {
431 : 0 : vertices[num_vtx] = 9;
432 : 0 : derivs[num_vtx][0] = 0.0;
433 : 0 : derivs[num_vtx][1] = 0.0;
434 : 0 : derivs[num_vtx][2] = 4.0;
435 : 0 : derivs[2][2] -= 2.0;
436 : 0 : derivs[3][2] -= 2.0;
437 : 0 : ++num_vtx;
438 : : }
439 : 0 : break;
440 : :
441 : : case 3:
442 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
443 : : {
444 : 0 : vertices[num_vtx] = 7;
445 : 0 : derivs[num_vtx][0] = -4.0;
446 : 0 : derivs[num_vtx][1] = -4.0;
447 : 0 : derivs[num_vtx][2] = -4.0;
448 : 0 : derivs[0][0] += 2.0;
449 : 0 : derivs[0][1] += 2.0;
450 : 0 : derivs[0][2] += 2.0;
451 : 0 : derivs[3][0] += 2.0;
452 : 0 : derivs[3][1] += 2.0;
453 : 0 : derivs[3][2] += 2.0;
454 : 0 : ++num_vtx;
455 : : }
456 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
457 : : {
458 : 0 : vertices[num_vtx] = 8;
459 : 0 : derivs[num_vtx][0] = 4.0;
460 : 0 : derivs[num_vtx][1] = 0.0;
461 : 0 : derivs[num_vtx][2] = 0.0;
462 : 0 : derivs[1][0] -= 2.0;
463 : 0 : derivs[3][0] -= 2.0;
464 : 0 : ++num_vtx;
465 : : }
466 : :
467 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
468 : : {
469 : 0 : vertices[num_vtx] = 9;
470 : 0 : derivs[num_vtx][0] = 0.0;
471 : 0 : derivs[num_vtx][1] = 4.0;
472 : 0 : derivs[num_vtx][2] = 0.0;
473 : 0 : derivs[2][1] -= 2.0;
474 : 0 : derivs[3][1] -= 2.0;
475 : 0 : ++num_vtx;
476 : : }
477 : 0 : break;
478 : : }
479 : 0 : }
480 : :
481 : 0 : static void derivatives_at_mid_edge( unsigned edge, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs,
482 : : size_t& num_vtx )
483 : : {
484 : : int sign;
485 : 0 : num_vtx = 2;
486 [ # # # # : 0 : switch( edge )
# # # ]
487 : : {
488 : : case 0:
489 : 0 : vertices[0] = 0;
490 : 0 : derivs[0][0] = -1.0;
491 : 0 : derivs[0][1] = -1.0;
492 : 0 : derivs[0][2] = -1.0;
493 : :
494 : 0 : vertices[1] = 1;
495 : 0 : derivs[1][0] = 1.0;
496 : 0 : derivs[1][1] = 0.0;
497 : 0 : derivs[1][2] = 0.0;
498 : :
499 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 2 ) )
500 : : {
501 : 0 : vertices[num_vtx] = 2;
502 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 1 );
503 : 0 : derivs[num_vtx][0] = 0.0;
504 : 0 : derivs[num_vtx][1] = sign;
505 : 0 : derivs[num_vtx][2] = 0.0;
506 : 0 : ++num_vtx;
507 : : }
508 : :
509 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) == nodeset.mid_edge_node( 4 ) )
510 : : {
511 : 0 : vertices[num_vtx] = 3;
512 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 3 );
513 : 0 : derivs[num_vtx][0] = 0.0;
514 : 0 : derivs[num_vtx][1] = 0.0;
515 : 0 : derivs[num_vtx][2] = sign;
516 : 0 : ++num_vtx;
517 : : }
518 : :
519 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
520 : : {
521 : 0 : vertices[num_vtx] = 4;
522 : 0 : derivs[num_vtx][0] = 0.0;
523 : 0 : derivs[num_vtx][1] = -2.0;
524 : 0 : derivs[num_vtx][2] = -2.0;
525 : 0 : derivs[0][1] += 1.0;
526 : 0 : derivs[0][2] += 1.0;
527 : 0 : derivs[1][1] += 1.0;
528 : 0 : derivs[1][2] += 1.0;
529 : 0 : ++num_vtx;
530 : : }
531 : :
532 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
533 : : {
534 : 0 : vertices[num_vtx] = 5;
535 : 0 : derivs[num_vtx][0] = 0.0;
536 : 0 : derivs[num_vtx][1] = 2.0;
537 : 0 : derivs[num_vtx][2] = 0.0;
538 : 0 : derivs[1][1] -= 1.0;
539 : 0 : ++num_vtx;
540 : : }
541 : :
542 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
543 : : {
544 : 0 : vertices[num_vtx] = 6;
545 : 0 : derivs[num_vtx][0] = 0.0;
546 : 0 : derivs[num_vtx][1] = 2.0;
547 : 0 : derivs[num_vtx][2] = 0.0;
548 : 0 : derivs[0][1] -= 1.0;
549 : 0 : ++num_vtx;
550 : : }
551 : :
552 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
553 : : {
554 : 0 : vertices[num_vtx] = 7;
555 : 0 : derivs[num_vtx][0] = 0.0;
556 : 0 : derivs[num_vtx][1] = 0.0;
557 : 0 : derivs[num_vtx][2] = 2.0;
558 : 0 : derivs[0][2] -= 1.0;
559 : 0 : ++num_vtx;
560 : : }
561 : :
562 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
563 : : {
564 : 0 : vertices[num_vtx] = 8;
565 : 0 : derivs[num_vtx][0] = 0.0;
566 : 0 : derivs[num_vtx][1] = 0.0;
567 : 0 : derivs[num_vtx][2] = 2.0;
568 : 0 : derivs[1][2] -= 1.0;
569 : 0 : ++num_vtx;
570 : : }
571 : 0 : break;
572 : :
573 : : case 1:
574 : 0 : vertices[0] = 1;
575 : 0 : derivs[0][0] = 1.0;
576 : 0 : derivs[0][1] = 0.0;
577 : 0 : derivs[0][2] = 0.0;
578 : :
579 : 0 : vertices[1] = 2;
580 : 0 : derivs[1][0] = 0.0;
581 : 0 : derivs[1][1] = 1.0;
582 : 0 : derivs[1][2] = 0.0;
583 : :
584 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 2 ) )
585 : : {
586 : 0 : vertices[num_vtx] = 0;
587 : 0 : sign = 2 * nodeset.mid_edge_node( 0 ) - 1;
588 : 0 : derivs[num_vtx][0] = sign;
589 : 0 : derivs[num_vtx][1] = sign;
590 : 0 : derivs[num_vtx][2] = sign;
591 : 0 : ++num_vtx;
592 : : }
593 : :
594 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) == nodeset.mid_edge_node( 5 ) )
595 : : {
596 : 0 : vertices[num_vtx] = 3;
597 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 4 );
598 : 0 : derivs[num_vtx][0] = 0.0;
599 : 0 : derivs[num_vtx][1] = 0.0;
600 : 0 : derivs[num_vtx][2] = sign;
601 : 0 : ++num_vtx;
602 : : }
603 : :
604 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
605 : : {
606 : 0 : vertices[num_vtx] = 4;
607 : 0 : derivs[num_vtx][0] = -2.0;
608 : 0 : derivs[num_vtx][1] = -2.0;
609 : 0 : derivs[num_vtx][2] = -2.0;
610 : 0 : ++num_vtx;
611 : 0 : derivs[0][0] += 1.0;
612 : 0 : derivs[0][1] += 1.0;
613 : 0 : derivs[0][2] += 1.0;
614 : : }
615 : :
616 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
617 : : {
618 : 0 : vertices[num_vtx] = 5;
619 : 0 : derivs[num_vtx][0] = 2.0;
620 : 0 : derivs[num_vtx][1] = 2.0;
621 : 0 : derivs[num_vtx][2] = 0.0;
622 : 0 : ++num_vtx;
623 : 0 : derivs[0][0] -= 1.0;
624 : 0 : derivs[0][1] -= 1.0;
625 : 0 : derivs[1][0] -= 1.0;
626 : 0 : derivs[1][1] -= 1.0;
627 : : }
628 : :
629 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
630 : : {
631 : 0 : vertices[num_vtx] = 6;
632 : 0 : derivs[num_vtx][0] = -2.0;
633 : 0 : derivs[num_vtx][1] = -2.0;
634 : 0 : derivs[num_vtx][2] = -2.0;
635 : 0 : ++num_vtx;
636 : 0 : derivs[1][0] += 1.0;
637 : 0 : derivs[1][1] += 1.0;
638 : 0 : derivs[1][2] += 1.0;
639 : : }
640 : :
641 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
642 : : {
643 : 0 : vertices[num_vtx] = 8;
644 : 0 : derivs[num_vtx][0] = 0.0;
645 : 0 : derivs[num_vtx][1] = 0.0;
646 : 0 : derivs[num_vtx][2] = 2.0;
647 : 0 : ++num_vtx;
648 : 0 : derivs[0][2] -= 1.0;
649 : : }
650 : :
651 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
652 : : {
653 : 0 : vertices[num_vtx] = 9;
654 : 0 : derivs[num_vtx][0] = 0.0;
655 : 0 : derivs[num_vtx][1] = 0.0;
656 : 0 : derivs[num_vtx][2] = 2.0;
657 : 0 : ++num_vtx;
658 : 0 : derivs[1][2] -= 1.0;
659 : : }
660 : 0 : break;
661 : :
662 : : case 2:
663 : 0 : vertices[0] = 0;
664 : 0 : derivs[0][0] = -1.0;
665 : 0 : derivs[0][1] = -1.0;
666 : 0 : derivs[0][2] = -1.0;
667 : :
668 : 0 : vertices[1] = 2;
669 : 0 : derivs[1][0] = 0.0;
670 : 0 : derivs[1][1] = 1.0;
671 : 0 : derivs[1][2] = 0.0;
672 : :
673 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 1 ) )
674 : : {
675 : 0 : vertices[num_vtx] = 1;
676 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 0 );
677 : 0 : derivs[num_vtx][0] = sign;
678 : 0 : derivs[num_vtx][1] = 0.0;
679 : 0 : derivs[num_vtx][2] = 0.0;
680 : 0 : ++num_vtx;
681 : : }
682 : :
683 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) == nodeset.mid_edge_node( 5 ) )
684 : : {
685 : 0 : vertices[num_vtx] = 3;
686 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 3 );
687 : 0 : derivs[num_vtx][0] = 0.0;
688 : 0 : derivs[num_vtx][1] = 0.0;
689 : 0 : derivs[num_vtx][2] = sign;
690 : 0 : ++num_vtx;
691 : : }
692 : :
693 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
694 : : {
695 : 0 : vertices[num_vtx] = 4;
696 : 0 : derivs[num_vtx][0] = 2.0;
697 : 0 : derivs[num_vtx][1] = 0.0;
698 : 0 : derivs[num_vtx][2] = 0.0;
699 : 0 : ++num_vtx;
700 : 0 : derivs[0][0] -= 1.0;
701 : : }
702 : :
703 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
704 : : {
705 : 0 : vertices[num_vtx] = 5;
706 : 0 : derivs[num_vtx][0] = 2.0;
707 : 0 : derivs[num_vtx][1] = 0.0;
708 : 0 : derivs[num_vtx][2] = 0.0;
709 : 0 : ++num_vtx;
710 : 0 : derivs[1][0] -= 1.0;
711 : : }
712 : :
713 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
714 : : {
715 : 0 : vertices[num_vtx] = 6;
716 : 0 : derivs[num_vtx][0] = -2.0;
717 : 0 : derivs[num_vtx][1] = 0.0;
718 : 0 : derivs[num_vtx][2] = -2.0;
719 : 0 : ++num_vtx;
720 : 0 : derivs[0][0] += 1.0;
721 : 0 : derivs[0][2] += 1.0;
722 : 0 : derivs[1][0] += 1.0;
723 : 0 : derivs[1][2] += 1.0;
724 : : }
725 : :
726 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
727 : : {
728 : 0 : vertices[num_vtx] = 7;
729 : 0 : derivs[num_vtx][0] = 0.0;
730 : 0 : derivs[num_vtx][1] = 0.0;
731 : 0 : derivs[num_vtx][2] = 2.0;
732 : 0 : ++num_vtx;
733 : 0 : derivs[0][2] -= 1.0;
734 : : }
735 : :
736 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
737 : : {
738 : 0 : vertices[num_vtx] = 9;
739 : 0 : derivs[num_vtx][0] = 0.0;
740 : 0 : derivs[num_vtx][1] = 0.0;
741 : 0 : derivs[num_vtx][2] = 2.0;
742 : 0 : ++num_vtx;
743 : 0 : derivs[1][2] -= 1.0;
744 : : }
745 : 0 : break;
746 : :
747 : : case 3:
748 : 0 : vertices[0] = 0;
749 : 0 : derivs[0][0] = -1.0;
750 : 0 : derivs[0][1] = -1.0;
751 : 0 : derivs[0][2] = -1.0;
752 : :
753 : 0 : vertices[1] = 3;
754 : 0 : derivs[1][0] = 0.0;
755 : 0 : derivs[1][1] = 0.0;
756 : 0 : derivs[1][2] = 1.0;
757 : :
758 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 4 ) )
759 : : {
760 : 0 : vertices[num_vtx] = 1;
761 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 0 );
762 : 0 : derivs[num_vtx][0] = sign;
763 : 0 : derivs[num_vtx][1] = 0.0;
764 : 0 : derivs[num_vtx][2] = 0.0;
765 : 0 : ++num_vtx;
766 : : }
767 : :
768 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) == nodeset.mid_edge_node( 5 ) )
769 : : {
770 : 0 : vertices[num_vtx] = 2;
771 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 2 );
772 : 0 : derivs[num_vtx][0] = 0.0;
773 : 0 : derivs[num_vtx][1] = sign;
774 : 0 : derivs[num_vtx][2] = 0.0;
775 : 0 : ++num_vtx;
776 : : }
777 : :
778 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
779 : : {
780 : 0 : vertices[num_vtx] = 4;
781 : 0 : derivs[num_vtx][0] = 2.0;
782 : 0 : derivs[num_vtx][1] = 0.0;
783 : 0 : derivs[num_vtx][2] = 0.0;
784 : 0 : ++num_vtx;
785 : 0 : derivs[0][0] -= 1.0;
786 : : }
787 : :
788 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
789 : : {
790 : 0 : vertices[num_vtx] = 6;
791 : 0 : derivs[num_vtx][0] = 0.0;
792 : 0 : derivs[num_vtx][1] = 2.0;
793 : 0 : derivs[num_vtx][2] = 0.0;
794 : 0 : ++num_vtx;
795 : 0 : derivs[0][1] -= 1.0;
796 : : }
797 : :
798 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
799 : : {
800 : 0 : vertices[num_vtx] = 7;
801 : 0 : derivs[num_vtx][0] = -2.0;
802 : 0 : derivs[num_vtx][1] = -2.0;
803 : 0 : derivs[num_vtx][2] = 0.0;
804 : 0 : ++num_vtx;
805 : 0 : derivs[0][0] += 1.0;
806 : 0 : derivs[0][1] += 1.0;
807 : 0 : derivs[1][0] += 1.0;
808 : 0 : derivs[1][1] += 1.0;
809 : : }
810 : :
811 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
812 : : {
813 : 0 : vertices[num_vtx] = 8;
814 : 0 : derivs[num_vtx][0] = 2.0;
815 : 0 : derivs[num_vtx][1] = 0.0;
816 : 0 : derivs[num_vtx][2] = 0.0;
817 : 0 : ++num_vtx;
818 : 0 : derivs[1][0] -= 1.0;
819 : : }
820 : :
821 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
822 : : {
823 : 0 : vertices[num_vtx] = 9;
824 : 0 : derivs[num_vtx][0] = 0.0;
825 : 0 : derivs[num_vtx][1] = 2.0;
826 : 0 : derivs[num_vtx][2] = 0.0;
827 : 0 : ++num_vtx;
828 : 0 : derivs[1][1] -= 1.0;
829 : : }
830 : 0 : break;
831 : :
832 : : case 4:
833 : 0 : vertices[0] = 1;
834 : 0 : derivs[0][0] = 1.0;
835 : 0 : derivs[0][1] = 0.0;
836 : 0 : derivs[0][2] = 0.0;
837 : :
838 : 0 : vertices[1] = 3;
839 : 0 : derivs[1][0] = 0.0;
840 : 0 : derivs[1][1] = 0.0;
841 : 0 : derivs[1][2] = 1.0;
842 : :
843 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 3 ) )
844 : : {
845 : 0 : vertices[num_vtx] = 0;
846 : 0 : sign = 2 * nodeset.mid_edge_node( 0 ) - 1;
847 : 0 : derivs[num_vtx][0] = sign;
848 : 0 : derivs[num_vtx][1] = sign;
849 : 0 : derivs[num_vtx][2] = sign;
850 : 0 : ++num_vtx;
851 : : }
852 : :
853 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 5 ) )
854 : : {
855 : 0 : vertices[num_vtx] = 2;
856 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 1 );
857 : 0 : derivs[num_vtx][0] = 0.0;
858 : 0 : derivs[num_vtx][1] = sign;
859 : 0 : derivs[num_vtx][2] = 0.0;
860 : 0 : ++num_vtx;
861 : : }
862 : :
863 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
864 : : {
865 : 0 : vertices[num_vtx] = 4;
866 : 0 : derivs[num_vtx][0] = -2.0;
867 : 0 : derivs[num_vtx][1] = -2.0;
868 : 0 : derivs[num_vtx][2] = -2.0;
869 : 0 : ++num_vtx;
870 : 0 : derivs[0][0] += 1.0;
871 : 0 : derivs[0][1] += 1.0;
872 : 0 : derivs[0][2] += 1.0;
873 : : }
874 : :
875 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
876 : : {
877 : 0 : vertices[num_vtx] = 5;
878 : 0 : derivs[num_vtx][0] = 0.0;
879 : 0 : derivs[num_vtx][1] = 2.0;
880 : 0 : derivs[num_vtx][2] = 0.0;
881 : 0 : ++num_vtx;
882 : 0 : derivs[0][1] -= 1.0;
883 : : }
884 : :
885 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
886 : : {
887 : 0 : vertices[num_vtx] = 7;
888 : 0 : derivs[num_vtx][0] = -2.0;
889 : 0 : derivs[num_vtx][1] = -2.0;
890 : 0 : derivs[num_vtx][2] = -2.0;
891 : 0 : ++num_vtx;
892 : 0 : derivs[1][0] += 1.0;
893 : 0 : derivs[1][1] += 1.0;
894 : 0 : derivs[1][2] += 1.0;
895 : : }
896 : :
897 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
898 : : {
899 : 0 : vertices[num_vtx] = 8;
900 : 0 : derivs[num_vtx][0] = 2.0;
901 : 0 : derivs[num_vtx][1] = 0.0;
902 : 0 : derivs[num_vtx][2] = 2.0;
903 : 0 : ++num_vtx;
904 : 0 : derivs[0][0] -= 1.0;
905 : 0 : derivs[0][2] -= 1.0;
906 : 0 : derivs[1][0] -= 1.0;
907 : 0 : derivs[1][2] -= 1.0;
908 : : }
909 : :
910 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
911 : : {
912 : 0 : vertices[num_vtx] = 9;
913 : 0 : derivs[num_vtx][0] = 0.0;
914 : 0 : derivs[num_vtx][1] = 2.0;
915 : 0 : derivs[num_vtx][2] = 0.0;
916 : 0 : ++num_vtx;
917 : 0 : derivs[1][1] -= 1.0;
918 : : }
919 : 0 : break;
920 : :
921 : : case 5:
922 : 0 : vertices[0] = 2;
923 : 0 : derivs[0][0] = 0.0;
924 : 0 : derivs[0][1] = 1.0;
925 : 0 : derivs[0][2] = 0.0;
926 : :
927 : 0 : vertices[1] = 3;
928 : 0 : derivs[1][0] = 0.0;
929 : 0 : derivs[1][1] = 0.0;
930 : 0 : derivs[1][2] = 1.0;
931 : :
932 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) == nodeset.mid_edge_node( 3 ) )
933 : : {
934 : 0 : vertices[num_vtx] = 0;
935 : 0 : sign = 2 * nodeset.mid_edge_node( 2 ) - 1;
936 : 0 : derivs[num_vtx][0] = sign;
937 : 0 : derivs[num_vtx][1] = sign;
938 : 0 : derivs[num_vtx][2] = sign;
939 : 0 : ++num_vtx;
940 : : }
941 : :
942 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 4 ) )
943 : : {
944 : 0 : vertices[num_vtx] = 1;
945 : 0 : sign = 1 - 2 * nodeset.mid_edge_node( 1 );
946 : 0 : derivs[num_vtx][0] = sign;
947 : 0 : derivs[num_vtx][1] = 0.0;
948 : 0 : derivs[num_vtx][2] = 0.0;
949 : 0 : ++num_vtx;
950 : : }
951 : :
952 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
953 : : {
954 : 0 : vertices[num_vtx] = 5;
955 : 0 : derivs[num_vtx][0] = 2.0;
956 : 0 : derivs[num_vtx][1] = 0.0;
957 : 0 : derivs[num_vtx][2] = 0.0;
958 : 0 : ++num_vtx;
959 : 0 : derivs[0][0] -= 1.0;
960 : : }
961 : :
962 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
963 : : {
964 : 0 : vertices[num_vtx] = 6;
965 : 0 : derivs[num_vtx][0] = -2.0;
966 : 0 : derivs[num_vtx][1] = -2.0;
967 : 0 : derivs[num_vtx][2] = -2.0;
968 : 0 : ++num_vtx;
969 : 0 : derivs[0][0] += 1.0;
970 : 0 : derivs[0][1] += 1.0;
971 : 0 : derivs[0][2] += 1.0;
972 : : }
973 : :
974 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
975 : : {
976 : 0 : vertices[num_vtx] = 7;
977 : 0 : derivs[num_vtx][0] = -2.0;
978 : 0 : derivs[num_vtx][1] = -2.0;
979 : 0 : derivs[num_vtx][2] = -2.0;
980 : 0 : ++num_vtx;
981 : 0 : derivs[1][0] += 1.0;
982 : 0 : derivs[1][1] += 1.0;
983 : 0 : derivs[1][2] += 1.0;
984 : : }
985 : :
986 [ # # ]: 0 : if( nodeset.mid_edge_node( 4 ) )
987 : : {
988 : 0 : vertices[num_vtx] = 8;
989 : 0 : derivs[num_vtx][0] = 2.0;
990 : 0 : derivs[num_vtx][1] = 0.0;
991 : 0 : derivs[num_vtx][2] = 0.0;
992 : 0 : ++num_vtx;
993 : 0 : derivs[1][0] -= 1.0;
994 : : }
995 : :
996 [ # # ]: 0 : if( nodeset.mid_edge_node( 5 ) )
997 : : {
998 : 0 : vertices[num_vtx] = 9;
999 : 0 : derivs[num_vtx][0] = 0.0;
1000 : 0 : derivs[num_vtx][1] = 2.0;
1001 : 0 : derivs[num_vtx][2] = 2.0;
1002 : 0 : ++num_vtx;
1003 : 0 : derivs[0][1] -= 1.0;
1004 : 0 : derivs[0][2] -= 1.0;
1005 : 0 : derivs[1][1] -= 1.0;
1006 : 0 : derivs[1][2] -= 1.0;
1007 : : }
1008 : 0 : break;
1009 : : }
1010 : 0 : }
1011 : :
1012 : : // Derivatives of coefficients for mid-edge nodes
1013 : :
1014 : : const double ft = 4.0 / 3.0;
1015 : :
1016 : : const double ho_dr[6][4] = { { 0., -ft, ft, 0. }, { 0., ft, ft, ft }, { 0., -ft, -ft, -ft },
1017 : 31 : { -ft, -ft, -ft, 0. }, { ft, ft, ft, 0. }, { 0., 0., 0., 0. } };
1018 : :
1019 : : const double ho_ds[6][4] = { { -ft, -ft, 0., -ft }, { ft, ft, 0., ft }, { ft, -ft, 0., 0. },
1020 : 31 : { -ft, -ft, -ft, 0. }, { 0., 0., 0., 0. }, { ft, ft, ft, 0. } };
1021 : :
1022 : : const double ho_dt[6][4] = { { -ft, -ft, 0., -ft }, { 0., 0., 0., 0. }, { 0., -ft, -ft, -ft },
1023 : 31 : { 0., -ft, 0., ft }, { ft, ft, 0., ft }, { 0., ft, ft, ft } };
1024 : :
1025 : 0 : static void derivatives_at_mid_face( unsigned face, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs,
1026 : : size_t& num_vtx )
1027 : : {
1028 : : // begin with derivatives for linear tetrahedron
1029 : 0 : num_vtx = 4;
1030 : 0 : get_linear_derivatives( vertices, derivs );
1031 : :
1032 [ # # ]: 0 : for( unsigned i = 0; i < 6; ++i )
1033 [ # # ]: 0 : if( nodeset.mid_edge_node( i ) )
1034 : : {
1035 : 0 : vertices[num_vtx] = i + 4;
1036 : 0 : derivs[num_vtx][0] = ho_dr[i][face];
1037 : 0 : derivs[num_vtx][1] = ho_ds[i][face];
1038 : 0 : derivs[num_vtx][2] = ho_dt[i][face];
1039 : 0 : ++num_vtx;
1040 : 0 : int j = edges[i][0];
1041 : 0 : derivs[j][0] -= 0.5 * ho_dr[i][face];
1042 : 0 : derivs[j][1] -= 0.5 * ho_ds[i][face];
1043 : 0 : derivs[j][2] -= 0.5 * ho_dt[i][face];
1044 : 0 : j = edges[i][1];
1045 : 0 : derivs[j][0] -= 0.5 * ho_dr[i][face];
1046 : 0 : derivs[j][1] -= 0.5 * ho_ds[i][face];
1047 : 0 : derivs[j][2] -= 0.5 * ho_dt[i][face];
1048 : : }
1049 : 0 : }
1050 : :
1051 : : // position (0->r, 1->s, 2->t) of zero-valued term for mid-edge node
1052 : : static const int zeros[6] = { 0, 2, 1, 2, 1, 0 };
1053 : : // value of mid-edge terms
1054 : : static const int signs[6] = { -1, 1, -1, -1, 1, 1 };
1055 : :
1056 : 0 : static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs, size_t& num_vtx )
1057 : : {
1058 : :
1059 : 0 : bool corners[4] = { false, false, false, false };
1060 : 0 : double corner_vals[4][3] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
1061 : :
1062 : 0 : num_vtx = 0;
1063 [ # # ]: 0 : for( unsigned i = 4; i < 10; ++i )
1064 : : {
1065 : 0 : int sign = signs[i - 4];
1066 : 0 : int zero = zeros[i - 4];
1067 : :
1068 [ # # ][ # # ]: 0 : if( nodeset.mid_edge_node( i - 4 ) )
1069 : : {
1070 : 0 : vertices[num_vtx] = i;
1071 [ # # ]: 0 : derivs[num_vtx][0] = (double)sign;
1072 [ # # ]: 0 : derivs[num_vtx][1] = (double)sign;
1073 [ # # ]: 0 : derivs[num_vtx][2] = (double)sign;
1074 [ # # ]: 0 : derivs[num_vtx][zero] = 0.0;
1075 : 0 : ++num_vtx;
1076 : : }
1077 : : else
1078 : : {
1079 [ # # ]: 0 : for( unsigned j = 0; j < 2; ++j )
1080 : : {
1081 : 0 : int corner = edges[i - 4][j];
1082 : 0 : int v1 = ( zero + 1 ) % 3;
1083 : 0 : int v2 = ( zero + 2 ) % 3;
1084 : 0 : corners[corner] = true;
1085 : 0 : corner_vals[corner][v1] += 0.5 * sign;
1086 : 0 : corner_vals[corner][v2] += 0.5 * sign;
1087 : : }
1088 : : }
1089 : : }
1090 : :
1091 [ # # ]: 0 : for( unsigned i = 0; i < 4; ++i )
1092 [ # # ]: 0 : if( corners[i] )
1093 : : {
1094 : 0 : vertices[num_vtx] = i;
1095 [ # # ]: 0 : derivs[num_vtx][0] = corner_vals[i][0];
1096 [ # # ]: 0 : derivs[num_vtx][1] = corner_vals[i][1];
1097 [ # # ]: 0 : derivs[num_vtx][2] = corner_vals[i][2];
1098 : 0 : ++num_vtx;
1099 : : }
1100 : 0 : }
1101 : :
1102 : 4565348 : void TetLagrangeShape::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
1103 : : MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
1104 : : {
1105 [ + - ]: 4565348 : if( !nodeset.have_any_mid_node() )
1106 : : {
1107 : 4565348 : num_vtx = 4;
1108 : 4565348 : get_linear_derivatives( vertex_indices_out, d_coeff_d_xi_out );
1109 : 4565348 : return;
1110 : : }
1111 : :
1112 [ # # ]: 0 : if( nodeset.have_any_mid_face_node() | nodeset.have_any_mid_region_node() )
1113 : : {
1114 : : MSQ_SETERR( err )
1115 [ # # ]: 0 : ( "TetLagrangeShape does not support mid-face/mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
1116 : 0 : return;
1117 : : }
1118 : :
1119 [ # # # # : 0 : switch( loc.dimension )
# ]
1120 : : {
1121 : : case 0:
1122 : 0 : derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
1123 : 0 : break;
1124 : : case 1:
1125 : 0 : derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
1126 : 0 : break;
1127 : : case 2:
1128 : 0 : derivatives_at_mid_face( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
1129 : 0 : break;
1130 : : case 3:
1131 : 0 : derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
1132 : 0 : break;
1133 : : default:
1134 [ # # ]: 4565348 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
1135 : : }
1136 : : }
1137 : :
1138 : 0 : void TetLagrangeShape::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const
1139 : : {
1140 : 0 : const double a = 1.122462048309373; // 6th root of 2
1141 : 0 : const double b = 1.7320508075688772; // sqrt(3)
1142 : 0 : const double c = 1.5874010519681994; // 2 to the 2/3
1143 : 0 : J( 0, 0 ) = a;
1144 : 0 : J( 0, 1 ) = 0.5 * a;
1145 : 0 : J( 0, 2 ) = 0.5 * a;
1146 : 0 : J( 1, 0 ) = 0.0;
1147 : 0 : J( 1, 1 ) = 0.5 * a * b;
1148 : 0 : J( 1, 2 ) = 0.5 * a / b;
1149 : 0 : J( 2, 0 ) = 0.0;
1150 : 0 : J( 2, 1 ) = 0.0;
1151 : 0 : J( 2, 2 ) = c / b;
1152 : 0 : }
1153 : :
1154 [ + - ][ + - ]: 124 : } // namespace MBMesquite
|
