MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2006 Sandia National Laboratories. Developed at the 00005 University of Wisconsin--Madison under SNL contract number 00006 624796. The U.S. Government and the University of Wisconsin 00007 retain certain rights to this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 (2006) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TetLagrangeShape.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TetLagrangeShape.hpp" 00034 #include "MsqError.hpp" 00035 #include <cassert> 00036 00037 namespace MBMesquite 00038 { 00039 00040 EntityTopology TetLagrangeShape::element_topology() const 00041 { 00042 return TETRAHEDRON; 00043 } 00044 00045 int TetLagrangeShape::num_nodes() const 00046 { 00047 return 10; 00048 } 00049 00050 NodeSet TetLagrangeShape::sample_points( NodeSet ns ) const 00051 { 00052 if( ns.have_any_mid_node() ) { ns.set_all_corner_nodes( TETRAHEDRON ); } 00053 else 00054 { 00055 ns.clear(); 00056 ns.set_mid_region_node(); 00057 } 00058 return ns; 00059 } 00060 00061 static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00062 { 00063 num_coeff = 1; 00064 indices_out[0] = corner; 00065 coeff_out[0] = 1.0; 00066 } 00067 00068 static void coefficients_at_mid_edge( unsigned edge, NodeSet nodeset, double* coeff_out, size_t* indices_out, 00069 size_t& num_coeff ) 00070 { 00071 if( nodeset.mid_edge_node( edge ) ) 00072 { // if mid-edge node is present 00073 num_coeff = 1; 00074 indices_out[0] = 4 + edge; 00075 coeff_out[0] = 1.0; 00076 } 00077 else 00078 { // no mid node on edge 00079 num_coeff = 2; 00080 coeff_out[0] = coeff_out[1] = 0.5; 00081 if( edge < 3 ) 00082 { 00083 indices_out[0] = edge; 00084 indices_out[1] = ( edge + 1 ) % 3; 00085 } 00086 else 00087 { 00088 indices_out[0] = edge - 3; 00089 indices_out[1] = 3; 00090 } 00091 } 00092 } 00093 00094 static void coefficients_at_mid_face( unsigned face, NodeSet nodeset, double* coeff_out, size_t* indices_out, 00095 size_t& num_coeff ) 00096 { 00097 const double one_ninth = 1.0 / 9.0; 00098 const double two_ninth = 2.0 / 9.0; 00099 const double four_ninth = 4.0 / 9.0; 00100 00101 if( face < 3 ) 00102 { 00103 const int next = ( face + 1 ) % 3; 00104 indices_out[0] = face; 00105 indices_out[1] = next; 00106 indices_out[2] = 3; 00107 coeff_out[0] = -one_ninth; 00108 coeff_out[1] = -one_ninth; 00109 coeff_out[2] = -one_ninth; 00110 num_coeff = 3; 00111 if( nodeset.mid_edge_node( face ) ) 00112 { 00113 indices_out[num_coeff] = 4 + face; 00114 coeff_out[num_coeff] = four_ninth; 00115 ++num_coeff; 00116 } 00117 else 00118 { 00119 coeff_out[0] += two_ninth; 00120 coeff_out[1] += two_ninth; 00121 } 00122 if( nodeset.mid_edge_node( 3 + next ) ) 00123 { 00124 indices_out[num_coeff] = 7 + next; 00125 coeff_out[num_coeff] = four_ninth; 00126 ++num_coeff; 00127 } 00128 else 00129 { 00130 coeff_out[1] += two_ninth; 00131 coeff_out[2] += two_ninth; 00132 } 00133 if( nodeset.mid_edge_node( 3 + face ) ) 00134 { 00135 indices_out[num_coeff] = 7 + face; 00136 coeff_out[num_coeff] = four_ninth; 00137 ++num_coeff; 00138 } 00139 else 00140 { 00141 coeff_out[0] += two_ninth; 00142 coeff_out[2] += two_ninth; 00143 } 00144 } 00145 else 00146 { 00147 assert( face == 3 ); 00148 indices_out[0] = 0; 00149 indices_out[1] = 1; 00150 indices_out[2] = 2; 00151 coeff_out[0] = -one_ninth; 00152 coeff_out[1] = -one_ninth; 00153 coeff_out[2] = -one_ninth; 00154 num_coeff = 3; 00155 if( nodeset.mid_edge_node( 0 ) ) 00156 { 00157 indices_out[num_coeff] = 4; 00158 coeff_out[num_coeff] = four_ninth; 00159 ++num_coeff; 00160 } 00161 else 00162 { 00163 coeff_out[0] += two_ninth; 00164 coeff_out[1] += two_ninth; 00165 } 00166 if( nodeset.mid_edge_node( 1 ) ) 00167 { 00168 indices_out[num_coeff] = 5; 00169 coeff_out[num_coeff] = four_ninth; 00170 ++num_coeff; 00171 } 00172 else 00173 { 00174 coeff_out[1] += two_ninth; 00175 coeff_out[2] += two_ninth; 00176 } 00177 if( nodeset.mid_edge_node( 2 ) ) 00178 { 00179 indices_out[num_coeff] = 6; 00180 coeff_out[num_coeff] = four_ninth; 00181 ++num_coeff; 00182 } 00183 else 00184 { 00185 coeff_out[2] += two_ninth; 00186 coeff_out[0] += two_ninth; 00187 } 00188 } 00189 } 00190 00191 static void coefficients_at_mid_elem( NodeSet nodeset, double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00192 { 00193 num_coeff = 4; 00194 indices_out[0] = 0; 00195 indices_out[1] = 1; 00196 indices_out[2] = 2; 00197 indices_out[3] = 3; 00198 coeff_out[0] = -0.125; 00199 coeff_out[1] = -0.125; 00200 coeff_out[2] = -0.125; 00201 coeff_out[3] = -0.125; 00202 if( nodeset.mid_edge_node( 0 ) ) 00203 { 00204 indices_out[num_coeff] = 4; 00205 coeff_out[num_coeff] = 0.25; 00206 ++num_coeff; 00207 } 00208 else 00209 { 00210 coeff_out[0] += 0.125; 00211 coeff_out[1] += 0.125; 00212 } 00213 if( nodeset.mid_edge_node( 1 ) ) 00214 { 00215 indices_out[num_coeff] = 5; 00216 coeff_out[num_coeff] = 0.25; 00217 ++num_coeff; 00218 } 00219 else 00220 { 00221 coeff_out[1] += 0.125; 00222 coeff_out[2] += 0.125; 00223 } 00224 if( nodeset.mid_edge_node( 2 ) ) 00225 { 00226 indices_out[num_coeff] = 6; 00227 coeff_out[num_coeff] = 0.25; 00228 ++num_coeff; 00229 } 00230 else 00231 { 00232 coeff_out[2] += 0.125; 00233 coeff_out[0] += 0.125; 00234 } 00235 if( nodeset.mid_edge_node( 3 ) ) 00236 { 00237 indices_out[num_coeff] = 7; 00238 coeff_out[num_coeff] = 0.25; 00239 ++num_coeff; 00240 } 00241 else 00242 { 00243 coeff_out[0] += 0.125; 00244 coeff_out[3] += 0.125; 00245 } 00246 if( nodeset.mid_edge_node( 4 ) ) 00247 { 00248 indices_out[num_coeff] = 8; 00249 coeff_out[num_coeff] = 0.25; 00250 ++num_coeff; 00251 } 00252 else 00253 { 00254 coeff_out[1] += 0.125; 00255 coeff_out[3] += 0.125; 00256 } 00257 if( nodeset.mid_edge_node( 5 ) ) 00258 { 00259 indices_out[num_coeff] = 9; 00260 coeff_out[num_coeff] = 0.25; 00261 ++num_coeff; 00262 } 00263 else 00264 { 00265 coeff_out[2] += 0.125; 00266 coeff_out[3] += 0.125; 00267 } 00268 } 00269 00270 void TetLagrangeShape::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out, 00271 size_t& num_coeff, MsqError& err ) const 00272 { 00273 if( nodeset.have_any_mid_face_node() | nodeset.have_any_mid_region_node() ) 00274 { 00275 MSQ_SETERR( err ) 00276 ( "TetLagrangeShape does not support mid-face/mid-element nodes", MsqError::UNSUPPORTED_ELEMENT ); 00277 return; 00278 } 00279 00280 switch( loc.dimension ) 00281 { 00282 case 0: 00283 coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff ); 00284 break; 00285 case 1: 00286 coefficients_at_mid_edge( loc.number, nodeset, coeff_out, indices_out, num_coeff ); 00287 break; 00288 case 2: 00289 coefficients_at_mid_face( loc.number, nodeset, coeff_out, indices_out, num_coeff ); 00290 break; 00291 case 3: 00292 coefficients_at_mid_elem( nodeset, coeff_out, indices_out, num_coeff ); 00293 break; 00294 default: 00295 MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG ); 00296 } 00297 } 00298 00299 static void get_linear_derivatives( size_t* vertices, MsqVector< 3 >* derivs ) 00300 { 00301 vertices[0] = 0; 00302 derivs[0][0] = -1.0; 00303 derivs[0][1] = -1.0; 00304 derivs[0][2] = -1.0; 00305 00306 vertices[1] = 1; 00307 derivs[1][0] = 1.0; 00308 derivs[1][1] = 0.0; 00309 derivs[1][2] = 0.0; 00310 00311 vertices[2] = 2; 00312 derivs[2][0] = 0.0; 00313 derivs[2][1] = 1.0; 00314 derivs[2][2] = 0.0; 00315 00316 vertices[3] = 3; 00317 derivs[3][0] = 0.0; 00318 derivs[3][1] = 0.0; 00319 derivs[3][2] = 1.0; 00320 } 00321 00322 static const unsigned edges[][2] = { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 3 }, { 2, 3 } }; 00323 00324 static void derivatives_at_corner( unsigned corner, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs, 00325 size_t& num_vtx ) 00326 { 00327 // begin with derivatives for linear tetrahedron 00328 num_vtx = 4; 00329 get_linear_derivatives( vertices, derivs ); 00330 00331 // adjust for the presence of mid-edge nodes 00332 switch( corner ) 00333 { 00334 case 0: 00335 if( nodeset.mid_edge_node( 0 ) ) 00336 { 00337 vertices[num_vtx] = 4; 00338 derivs[num_vtx][0] = 4.0; 00339 derivs[num_vtx][1] = 0.0; 00340 derivs[num_vtx][2] = 0.0; 00341 derivs[0][0] -= 2.0; 00342 derivs[1][0] -= 2.0; 00343 ++num_vtx; 00344 } 00345 if( nodeset.mid_edge_node( 2 ) ) 00346 { 00347 vertices[num_vtx] = 6; 00348 derivs[num_vtx][0] = 0.0; 00349 derivs[num_vtx][1] = 4.0; 00350 derivs[num_vtx][2] = 0.0; 00351 derivs[0][1] -= 2.0; 00352 derivs[2][1] -= 2.0; 00353 ++num_vtx; 00354 } 00355 if( nodeset.mid_edge_node( 3 ) ) 00356 { 00357 vertices[num_vtx] = 7; 00358 derivs[num_vtx][0] = 0.0; 00359 derivs[num_vtx][1] = 0.0; 00360 derivs[num_vtx][2] = 4.0; 00361 derivs[0][2] -= 2.0; 00362 derivs[3][2] -= 2.0; 00363 ++num_vtx; 00364 } 00365 break; 00366 00367 case 1: 00368 if( nodeset.mid_edge_node( 0 ) ) 00369 { 00370 vertices[num_vtx] = 4; 00371 derivs[num_vtx][0] = -4.0; 00372 derivs[num_vtx][1] = -4.0; 00373 derivs[num_vtx][2] = -4.0; 00374 derivs[0][0] += 2.0; 00375 derivs[0][1] += 2.0; 00376 derivs[0][2] += 2.0; 00377 derivs[1][0] += 2.0; 00378 derivs[1][1] += 2.0; 00379 derivs[1][2] += 2.0; 00380 ++num_vtx; 00381 } 00382 if( nodeset.mid_edge_node( 1 ) ) 00383 { 00384 vertices[num_vtx] = 5; 00385 derivs[num_vtx][0] = 0.0; 00386 derivs[num_vtx][1] = 4.0; 00387 derivs[num_vtx][2] = 0.0; 00388 derivs[1][1] -= 2.0; 00389 derivs[2][1] -= 2.0; 00390 ++num_vtx; 00391 } 00392 if( nodeset.mid_edge_node( 4 ) ) 00393 { 00394 vertices[num_vtx] = 8; 00395 derivs[num_vtx][0] = 0.0; 00396 derivs[num_vtx][1] = 0.0; 00397 derivs[num_vtx][2] = 4.0; 00398 derivs[1][2] -= 2.0; 00399 derivs[3][2] -= 2.0; 00400 ++num_vtx; 00401 } 00402 break; 00403 00404 case 2: 00405 if( nodeset.mid_edge_node( 1 ) ) 00406 { 00407 vertices[num_vtx] = 5; 00408 derivs[num_vtx][0] = 4.0; 00409 derivs[num_vtx][1] = 0.0; 00410 derivs[num_vtx][2] = 0.0; 00411 derivs[1][0] -= 2.0; 00412 derivs[2][0] -= 2.0; 00413 ++num_vtx; 00414 } 00415 if( nodeset.mid_edge_node( 2 ) ) 00416 { 00417 vertices[num_vtx] = 6; 00418 derivs[num_vtx][0] = -4.0; 00419 derivs[num_vtx][1] = -4.0; 00420 derivs[num_vtx][2] = -4.0; 00421 derivs[0][0] += 2.0; 00422 derivs[0][1] += 2.0; 00423 derivs[0][2] += 2.0; 00424 derivs[2][0] += 2.0; 00425 derivs[2][1] += 2.0; 00426 derivs[2][2] += 2.0; 00427 ++num_vtx; 00428 } 00429 if( nodeset.mid_edge_node( 5 ) ) 00430 { 00431 vertices[num_vtx] = 9; 00432 derivs[num_vtx][0] = 0.0; 00433 derivs[num_vtx][1] = 0.0; 00434 derivs[num_vtx][2] = 4.0; 00435 derivs[2][2] -= 2.0; 00436 derivs[3][2] -= 2.0; 00437 ++num_vtx; 00438 } 00439 break; 00440 00441 case 3: 00442 if( nodeset.mid_edge_node( 3 ) ) 00443 { 00444 vertices[num_vtx] = 7; 00445 derivs[num_vtx][0] = -4.0; 00446 derivs[num_vtx][1] = -4.0; 00447 derivs[num_vtx][2] = -4.0; 00448 derivs[0][0] += 2.0; 00449 derivs[0][1] += 2.0; 00450 derivs[0][2] += 2.0; 00451 derivs[3][0] += 2.0; 00452 derivs[3][1] += 2.0; 00453 derivs[3][2] += 2.0; 00454 ++num_vtx; 00455 } 00456 if( nodeset.mid_edge_node( 4 ) ) 00457 { 00458 vertices[num_vtx] = 8; 00459 derivs[num_vtx][0] = 4.0; 00460 derivs[num_vtx][1] = 0.0; 00461 derivs[num_vtx][2] = 0.0; 00462 derivs[1][0] -= 2.0; 00463 derivs[3][0] -= 2.0; 00464 ++num_vtx; 00465 } 00466 00467 if( nodeset.mid_edge_node( 5 ) ) 00468 { 00469 vertices[num_vtx] = 9; 00470 derivs[num_vtx][0] = 0.0; 00471 derivs[num_vtx][1] = 4.0; 00472 derivs[num_vtx][2] = 0.0; 00473 derivs[2][1] -= 2.0; 00474 derivs[3][1] -= 2.0; 00475 ++num_vtx; 00476 } 00477 break; 00478 } 00479 } 00480 00481 static void derivatives_at_mid_edge( unsigned edge, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs, 00482 size_t& num_vtx ) 00483 { 00484 int sign; 00485 num_vtx = 2; 00486 switch( edge ) 00487 { 00488 case 0: 00489 vertices[0] = 0; 00490 derivs[0][0] = -1.0; 00491 derivs[0][1] = -1.0; 00492 derivs[0][2] = -1.0; 00493 00494 vertices[1] = 1; 00495 derivs[1][0] = 1.0; 00496 derivs[1][1] = 0.0; 00497 derivs[1][2] = 0.0; 00498 00499 if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 2 ) ) 00500 { 00501 vertices[num_vtx] = 2; 00502 sign = 1 - 2 * nodeset.mid_edge_node( 1 ); 00503 derivs[num_vtx][0] = 0.0; 00504 derivs[num_vtx][1] = sign; 00505 derivs[num_vtx][2] = 0.0; 00506 ++num_vtx; 00507 } 00508 00509 if( nodeset.mid_edge_node( 3 ) == nodeset.mid_edge_node( 4 ) ) 00510 { 00511 vertices[num_vtx] = 3; 00512 sign = 1 - 2 * nodeset.mid_edge_node( 3 ); 00513 derivs[num_vtx][0] = 0.0; 00514 derivs[num_vtx][1] = 0.0; 00515 derivs[num_vtx][2] = sign; 00516 ++num_vtx; 00517 } 00518 00519 if( nodeset.mid_edge_node( 0 ) ) 00520 { 00521 vertices[num_vtx] = 4; 00522 derivs[num_vtx][0] = 0.0; 00523 derivs[num_vtx][1] = -2.0; 00524 derivs[num_vtx][2] = -2.0; 00525 derivs[0][1] += 1.0; 00526 derivs[0][2] += 1.0; 00527 derivs[1][1] += 1.0; 00528 derivs[1][2] += 1.0; 00529 ++num_vtx; 00530 } 00531 00532 if( nodeset.mid_edge_node( 1 ) ) 00533 { 00534 vertices[num_vtx] = 5; 00535 derivs[num_vtx][0] = 0.0; 00536 derivs[num_vtx][1] = 2.0; 00537 derivs[num_vtx][2] = 0.0; 00538 derivs[1][1] -= 1.0; 00539 ++num_vtx; 00540 } 00541 00542 if( nodeset.mid_edge_node( 2 ) ) 00543 { 00544 vertices[num_vtx] = 6; 00545 derivs[num_vtx][0] = 0.0; 00546 derivs[num_vtx][1] = 2.0; 00547 derivs[num_vtx][2] = 0.0; 00548 derivs[0][1] -= 1.0; 00549 ++num_vtx; 00550 } 00551 00552 if( nodeset.mid_edge_node( 3 ) ) 00553 { 00554 vertices[num_vtx] = 7; 00555 derivs[num_vtx][0] = 0.0; 00556 derivs[num_vtx][1] = 0.0; 00557 derivs[num_vtx][2] = 2.0; 00558 derivs[0][2] -= 1.0; 00559 ++num_vtx; 00560 } 00561 00562 if( nodeset.mid_edge_node( 4 ) ) 00563 { 00564 vertices[num_vtx] = 8; 00565 derivs[num_vtx][0] = 0.0; 00566 derivs[num_vtx][1] = 0.0; 00567 derivs[num_vtx][2] = 2.0; 00568 derivs[1][2] -= 1.0; 00569 ++num_vtx; 00570 } 00571 break; 00572 00573 case 1: 00574 vertices[0] = 1; 00575 derivs[0][0] = 1.0; 00576 derivs[0][1] = 0.0; 00577 derivs[0][2] = 0.0; 00578 00579 vertices[1] = 2; 00580 derivs[1][0] = 0.0; 00581 derivs[1][1] = 1.0; 00582 derivs[1][2] = 0.0; 00583 00584 if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 2 ) ) 00585 { 00586 vertices[num_vtx] = 0; 00587 sign = 2 * nodeset.mid_edge_node( 0 ) - 1; 00588 derivs[num_vtx][0] = sign; 00589 derivs[num_vtx][1] = sign; 00590 derivs[num_vtx][2] = sign; 00591 ++num_vtx; 00592 } 00593 00594 if( nodeset.mid_edge_node( 4 ) == nodeset.mid_edge_node( 5 ) ) 00595 { 00596 vertices[num_vtx] = 3; 00597 sign = 1 - 2 * nodeset.mid_edge_node( 4 ); 00598 derivs[num_vtx][0] = 0.0; 00599 derivs[num_vtx][1] = 0.0; 00600 derivs[num_vtx][2] = sign; 00601 ++num_vtx; 00602 } 00603 00604 if( nodeset.mid_edge_node( 0 ) ) 00605 { 00606 vertices[num_vtx] = 4; 00607 derivs[num_vtx][0] = -2.0; 00608 derivs[num_vtx][1] = -2.0; 00609 derivs[num_vtx][2] = -2.0; 00610 ++num_vtx; 00611 derivs[0][0] += 1.0; 00612 derivs[0][1] += 1.0; 00613 derivs[0][2] += 1.0; 00614 } 00615 00616 if( nodeset.mid_edge_node( 1 ) ) 00617 { 00618 vertices[num_vtx] = 5; 00619 derivs[num_vtx][0] = 2.0; 00620 derivs[num_vtx][1] = 2.0; 00621 derivs[num_vtx][2] = 0.0; 00622 ++num_vtx; 00623 derivs[0][0] -= 1.0; 00624 derivs[0][1] -= 1.0; 00625 derivs[1][0] -= 1.0; 00626 derivs[1][1] -= 1.0; 00627 } 00628 00629 if( nodeset.mid_edge_node( 2 ) ) 00630 { 00631 vertices[num_vtx] = 6; 00632 derivs[num_vtx][0] = -2.0; 00633 derivs[num_vtx][1] = -2.0; 00634 derivs[num_vtx][2] = -2.0; 00635 ++num_vtx; 00636 derivs[1][0] += 1.0; 00637 derivs[1][1] += 1.0; 00638 derivs[1][2] += 1.0; 00639 } 00640 00641 if( nodeset.mid_edge_node( 4 ) ) 00642 { 00643 vertices[num_vtx] = 8; 00644 derivs[num_vtx][0] = 0.0; 00645 derivs[num_vtx][1] = 0.0; 00646 derivs[num_vtx][2] = 2.0; 00647 ++num_vtx; 00648 derivs[0][2] -= 1.0; 00649 } 00650 00651 if( nodeset.mid_edge_node( 5 ) ) 00652 { 00653 vertices[num_vtx] = 9; 00654 derivs[num_vtx][0] = 0.0; 00655 derivs[num_vtx][1] = 0.0; 00656 derivs[num_vtx][2] = 2.0; 00657 ++num_vtx; 00658 derivs[1][2] -= 1.0; 00659 } 00660 break; 00661 00662 case 2: 00663 vertices[0] = 0; 00664 derivs[0][0] = -1.0; 00665 derivs[0][1] = -1.0; 00666 derivs[0][2] = -1.0; 00667 00668 vertices[1] = 2; 00669 derivs[1][0] = 0.0; 00670 derivs[1][1] = 1.0; 00671 derivs[1][2] = 0.0; 00672 00673 if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 1 ) ) 00674 { 00675 vertices[num_vtx] = 1; 00676 sign = 1 - 2 * nodeset.mid_edge_node( 0 ); 00677 derivs[num_vtx][0] = sign; 00678 derivs[num_vtx][1] = 0.0; 00679 derivs[num_vtx][2] = 0.0; 00680 ++num_vtx; 00681 } 00682 00683 if( nodeset.mid_edge_node( 3 ) == nodeset.mid_edge_node( 5 ) ) 00684 { 00685 vertices[num_vtx] = 3; 00686 sign = 1 - 2 * nodeset.mid_edge_node( 3 ); 00687 derivs[num_vtx][0] = 0.0; 00688 derivs[num_vtx][1] = 0.0; 00689 derivs[num_vtx][2] = sign; 00690 ++num_vtx; 00691 } 00692 00693 if( nodeset.mid_edge_node( 0 ) ) 00694 { 00695 vertices[num_vtx] = 4; 00696 derivs[num_vtx][0] = 2.0; 00697 derivs[num_vtx][1] = 0.0; 00698 derivs[num_vtx][2] = 0.0; 00699 ++num_vtx; 00700 derivs[0][0] -= 1.0; 00701 } 00702 00703 if( nodeset.mid_edge_node( 1 ) ) 00704 { 00705 vertices[num_vtx] = 5; 00706 derivs[num_vtx][0] = 2.0; 00707 derivs[num_vtx][1] = 0.0; 00708 derivs[num_vtx][2] = 0.0; 00709 ++num_vtx; 00710 derivs[1][0] -= 1.0; 00711 } 00712 00713 if( nodeset.mid_edge_node( 2 ) ) 00714 { 00715 vertices[num_vtx] = 6; 00716 derivs[num_vtx][0] = -2.0; 00717 derivs[num_vtx][1] = 0.0; 00718 derivs[num_vtx][2] = -2.0; 00719 ++num_vtx; 00720 derivs[0][0] += 1.0; 00721 derivs[0][2] += 1.0; 00722 derivs[1][0] += 1.0; 00723 derivs[1][2] += 1.0; 00724 } 00725 00726 if( nodeset.mid_edge_node( 3 ) ) 00727 { 00728 vertices[num_vtx] = 7; 00729 derivs[num_vtx][0] = 0.0; 00730 derivs[num_vtx][1] = 0.0; 00731 derivs[num_vtx][2] = 2.0; 00732 ++num_vtx; 00733 derivs[0][2] -= 1.0; 00734 } 00735 00736 if( nodeset.mid_edge_node( 5 ) ) 00737 { 00738 vertices[num_vtx] = 9; 00739 derivs[num_vtx][0] = 0.0; 00740 derivs[num_vtx][1] = 0.0; 00741 derivs[num_vtx][2] = 2.0; 00742 ++num_vtx; 00743 derivs[1][2] -= 1.0; 00744 } 00745 break; 00746 00747 case 3: 00748 vertices[0] = 0; 00749 derivs[0][0] = -1.0; 00750 derivs[0][1] = -1.0; 00751 derivs[0][2] = -1.0; 00752 00753 vertices[1] = 3; 00754 derivs[1][0] = 0.0; 00755 derivs[1][1] = 0.0; 00756 derivs[1][2] = 1.0; 00757 00758 if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 4 ) ) 00759 { 00760 vertices[num_vtx] = 1; 00761 sign = 1 - 2 * nodeset.mid_edge_node( 0 ); 00762 derivs[num_vtx][0] = sign; 00763 derivs[num_vtx][1] = 0.0; 00764 derivs[num_vtx][2] = 0.0; 00765 ++num_vtx; 00766 } 00767 00768 if( nodeset.mid_edge_node( 2 ) == nodeset.mid_edge_node( 5 ) ) 00769 { 00770 vertices[num_vtx] = 2; 00771 sign = 1 - 2 * nodeset.mid_edge_node( 2 ); 00772 derivs[num_vtx][0] = 0.0; 00773 derivs[num_vtx][1] = sign; 00774 derivs[num_vtx][2] = 0.0; 00775 ++num_vtx; 00776 } 00777 00778 if( nodeset.mid_edge_node( 0 ) ) 00779 { 00780 vertices[num_vtx] = 4; 00781 derivs[num_vtx][0] = 2.0; 00782 derivs[num_vtx][1] = 0.0; 00783 derivs[num_vtx][2] = 0.0; 00784 ++num_vtx; 00785 derivs[0][0] -= 1.0; 00786 } 00787 00788 if( nodeset.mid_edge_node( 2 ) ) 00789 { 00790 vertices[num_vtx] = 6; 00791 derivs[num_vtx][0] = 0.0; 00792 derivs[num_vtx][1] = 2.0; 00793 derivs[num_vtx][2] = 0.0; 00794 ++num_vtx; 00795 derivs[0][1] -= 1.0; 00796 } 00797 00798 if( nodeset.mid_edge_node( 3 ) ) 00799 { 00800 vertices[num_vtx] = 7; 00801 derivs[num_vtx][0] = -2.0; 00802 derivs[num_vtx][1] = -2.0; 00803 derivs[num_vtx][2] = 0.0; 00804 ++num_vtx; 00805 derivs[0][0] += 1.0; 00806 derivs[0][1] += 1.0; 00807 derivs[1][0] += 1.0; 00808 derivs[1][1] += 1.0; 00809 } 00810 00811 if( nodeset.mid_edge_node( 4 ) ) 00812 { 00813 vertices[num_vtx] = 8; 00814 derivs[num_vtx][0] = 2.0; 00815 derivs[num_vtx][1] = 0.0; 00816 derivs[num_vtx][2] = 0.0; 00817 ++num_vtx; 00818 derivs[1][0] -= 1.0; 00819 } 00820 00821 if( nodeset.mid_edge_node( 5 ) ) 00822 { 00823 vertices[num_vtx] = 9; 00824 derivs[num_vtx][0] = 0.0; 00825 derivs[num_vtx][1] = 2.0; 00826 derivs[num_vtx][2] = 0.0; 00827 ++num_vtx; 00828 derivs[1][1] -= 1.0; 00829 } 00830 break; 00831 00832 case 4: 00833 vertices[0] = 1; 00834 derivs[0][0] = 1.0; 00835 derivs[0][1] = 0.0; 00836 derivs[0][2] = 0.0; 00837 00838 vertices[1] = 3; 00839 derivs[1][0] = 0.0; 00840 derivs[1][1] = 0.0; 00841 derivs[1][2] = 1.0; 00842 00843 if( nodeset.mid_edge_node( 0 ) == nodeset.mid_edge_node( 3 ) ) 00844 { 00845 vertices[num_vtx] = 0; 00846 sign = 2 * nodeset.mid_edge_node( 0 ) - 1; 00847 derivs[num_vtx][0] = sign; 00848 derivs[num_vtx][1] = sign; 00849 derivs[num_vtx][2] = sign; 00850 ++num_vtx; 00851 } 00852 00853 if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 5 ) ) 00854 { 00855 vertices[num_vtx] = 2; 00856 sign = 1 - 2 * nodeset.mid_edge_node( 1 ); 00857 derivs[num_vtx][0] = 0.0; 00858 derivs[num_vtx][1] = sign; 00859 derivs[num_vtx][2] = 0.0; 00860 ++num_vtx; 00861 } 00862 00863 if( nodeset.mid_edge_node( 0 ) ) 00864 { 00865 vertices[num_vtx] = 4; 00866 derivs[num_vtx][0] = -2.0; 00867 derivs[num_vtx][1] = -2.0; 00868 derivs[num_vtx][2] = -2.0; 00869 ++num_vtx; 00870 derivs[0][0] += 1.0; 00871 derivs[0][1] += 1.0; 00872 derivs[0][2] += 1.0; 00873 } 00874 00875 if( nodeset.mid_edge_node( 1 ) ) 00876 { 00877 vertices[num_vtx] = 5; 00878 derivs[num_vtx][0] = 0.0; 00879 derivs[num_vtx][1] = 2.0; 00880 derivs[num_vtx][2] = 0.0; 00881 ++num_vtx; 00882 derivs[0][1] -= 1.0; 00883 } 00884 00885 if( nodeset.mid_edge_node( 3 ) ) 00886 { 00887 vertices[num_vtx] = 7; 00888 derivs[num_vtx][0] = -2.0; 00889 derivs[num_vtx][1] = -2.0; 00890 derivs[num_vtx][2] = -2.0; 00891 ++num_vtx; 00892 derivs[1][0] += 1.0; 00893 derivs[1][1] += 1.0; 00894 derivs[1][2] += 1.0; 00895 } 00896 00897 if( nodeset.mid_edge_node( 4 ) ) 00898 { 00899 vertices[num_vtx] = 8; 00900 derivs[num_vtx][0] = 2.0; 00901 derivs[num_vtx][1] = 0.0; 00902 derivs[num_vtx][2] = 2.0; 00903 ++num_vtx; 00904 derivs[0][0] -= 1.0; 00905 derivs[0][2] -= 1.0; 00906 derivs[1][0] -= 1.0; 00907 derivs[1][2] -= 1.0; 00908 } 00909 00910 if( nodeset.mid_edge_node( 5 ) ) 00911 { 00912 vertices[num_vtx] = 9; 00913 derivs[num_vtx][0] = 0.0; 00914 derivs[num_vtx][1] = 2.0; 00915 derivs[num_vtx][2] = 0.0; 00916 ++num_vtx; 00917 derivs[1][1] -= 1.0; 00918 } 00919 break; 00920 00921 case 5: 00922 vertices[0] = 2; 00923 derivs[0][0] = 0.0; 00924 derivs[0][1] = 1.0; 00925 derivs[0][2] = 0.0; 00926 00927 vertices[1] = 3; 00928 derivs[1][0] = 0.0; 00929 derivs[1][1] = 0.0; 00930 derivs[1][2] = 1.0; 00931 00932 if( nodeset.mid_edge_node( 2 ) == nodeset.mid_edge_node( 3 ) ) 00933 { 00934 vertices[num_vtx] = 0; 00935 sign = 2 * nodeset.mid_edge_node( 2 ) - 1; 00936 derivs[num_vtx][0] = sign; 00937 derivs[num_vtx][1] = sign; 00938 derivs[num_vtx][2] = sign; 00939 ++num_vtx; 00940 } 00941 00942 if( nodeset.mid_edge_node( 1 ) == nodeset.mid_edge_node( 4 ) ) 00943 { 00944 vertices[num_vtx] = 1; 00945 sign = 1 - 2 * nodeset.mid_edge_node( 1 ); 00946 derivs[num_vtx][0] = sign; 00947 derivs[num_vtx][1] = 0.0; 00948 derivs[num_vtx][2] = 0.0; 00949 ++num_vtx; 00950 } 00951 00952 if( nodeset.mid_edge_node( 1 ) ) 00953 { 00954 vertices[num_vtx] = 5; 00955 derivs[num_vtx][0] = 2.0; 00956 derivs[num_vtx][1] = 0.0; 00957 derivs[num_vtx][2] = 0.0; 00958 ++num_vtx; 00959 derivs[0][0] -= 1.0; 00960 } 00961 00962 if( nodeset.mid_edge_node( 2 ) ) 00963 { 00964 vertices[num_vtx] = 6; 00965 derivs[num_vtx][0] = -2.0; 00966 derivs[num_vtx][1] = -2.0; 00967 derivs[num_vtx][2] = -2.0; 00968 ++num_vtx; 00969 derivs[0][0] += 1.0; 00970 derivs[0][1] += 1.0; 00971 derivs[0][2] += 1.0; 00972 } 00973 00974 if( nodeset.mid_edge_node( 3 ) ) 00975 { 00976 vertices[num_vtx] = 7; 00977 derivs[num_vtx][0] = -2.0; 00978 derivs[num_vtx][1] = -2.0; 00979 derivs[num_vtx][2] = -2.0; 00980 ++num_vtx; 00981 derivs[1][0] += 1.0; 00982 derivs[1][1] += 1.0; 00983 derivs[1][2] += 1.0; 00984 } 00985 00986 if( nodeset.mid_edge_node( 4 ) ) 00987 { 00988 vertices[num_vtx] = 8; 00989 derivs[num_vtx][0] = 2.0; 00990 derivs[num_vtx][1] = 0.0; 00991 derivs[num_vtx][2] = 0.0; 00992 ++num_vtx; 00993 derivs[1][0] -= 1.0; 00994 } 00995 00996 if( nodeset.mid_edge_node( 5 ) ) 00997 { 00998 vertices[num_vtx] = 9; 00999 derivs[num_vtx][0] = 0.0; 01000 derivs[num_vtx][1] = 2.0; 01001 derivs[num_vtx][2] = 2.0; 01002 ++num_vtx; 01003 derivs[0][1] -= 1.0; 01004 derivs[0][2] -= 1.0; 01005 derivs[1][1] -= 1.0; 01006 derivs[1][2] -= 1.0; 01007 } 01008 break; 01009 } 01010 } 01011 01012 // Derivatives of coefficients for mid-edge nodes 01013 01014 const double ft = 4.0 / 3.0; 01015 01016 const double ho_dr[6][4] = { { 0., -ft, ft, 0. }, { 0., ft, ft, ft }, { 0., -ft, -ft, -ft }, 01017 { -ft, -ft, -ft, 0. }, { ft, ft, ft, 0. }, { 0., 0., 0., 0. } }; 01018 01019 const double ho_ds[6][4] = { { -ft, -ft, 0., -ft }, { ft, ft, 0., ft }, { ft, -ft, 0., 0. }, 01020 { -ft, -ft, -ft, 0. }, { 0., 0., 0., 0. }, { ft, ft, ft, 0. } }; 01021 01022 const double ho_dt[6][4] = { { -ft, -ft, 0., -ft }, { 0., 0., 0., 0. }, { 0., -ft, -ft, -ft }, 01023 { 0., -ft, 0., ft }, { ft, ft, 0., ft }, { 0., ft, ft, ft } }; 01024 01025 static void derivatives_at_mid_face( unsigned face, NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs, 01026 size_t& num_vtx ) 01027 { 01028 // begin with derivatives for linear tetrahedron 01029 num_vtx = 4; 01030 get_linear_derivatives( vertices, derivs ); 01031 01032 for( unsigned i = 0; i < 6; ++i ) 01033 if( nodeset.mid_edge_node( i ) ) 01034 { 01035 vertices[num_vtx] = i + 4; 01036 derivs[num_vtx][0] = ho_dr[i][face]; 01037 derivs[num_vtx][1] = ho_ds[i][face]; 01038 derivs[num_vtx][2] = ho_dt[i][face]; 01039 ++num_vtx; 01040 int j = edges[i][0]; 01041 derivs[j][0] -= 0.5 * ho_dr[i][face]; 01042 derivs[j][1] -= 0.5 * ho_ds[i][face]; 01043 derivs[j][2] -= 0.5 * ho_dt[i][face]; 01044 j = edges[i][1]; 01045 derivs[j][0] -= 0.5 * ho_dr[i][face]; 01046 derivs[j][1] -= 0.5 * ho_ds[i][face]; 01047 derivs[j][2] -= 0.5 * ho_dt[i][face]; 01048 } 01049 } 01050 01051 // position (0->r, 1->s, 2->t) of zero-valued term for mid-edge node 01052 static const int zeros[6] = { 0, 2, 1, 2, 1, 0 }; 01053 // value of mid-edge terms 01054 static const int signs[6] = { -1, 1, -1, -1, 1, 1 }; 01055 01056 static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 3 >* derivs, size_t& num_vtx ) 01057 { 01058 01059 bool corners[4] = { false, false, false, false }; 01060 double corner_vals[4][3] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; 01061 01062 num_vtx = 0; 01063 for( unsigned i = 4; i < 10; ++i ) 01064 { 01065 int sign = signs[i - 4]; 01066 int zero = zeros[i - 4]; 01067 01068 if( nodeset.mid_edge_node( i - 4 ) ) 01069 { 01070 vertices[num_vtx] = i; 01071 derivs[num_vtx][0] = (double)sign; 01072 derivs[num_vtx][1] = (double)sign; 01073 derivs[num_vtx][2] = (double)sign; 01074 derivs[num_vtx][zero] = 0.0; 01075 ++num_vtx; 01076 } 01077 else 01078 { 01079 for( unsigned j = 0; j < 2; ++j ) 01080 { 01081 int corner = edges[i - 4][j]; 01082 int v1 = ( zero + 1 ) % 3; 01083 int v2 = ( zero + 2 ) % 3; 01084 corners[corner] = true; 01085 corner_vals[corner][v1] += 0.5 * sign; 01086 corner_vals[corner][v2] += 0.5 * sign; 01087 } 01088 } 01089 } 01090 01091 for( unsigned i = 0; i < 4; ++i ) 01092 if( corners[i] ) 01093 { 01094 vertices[num_vtx] = i; 01095 derivs[num_vtx][0] = corner_vals[i][0]; 01096 derivs[num_vtx][1] = corner_vals[i][1]; 01097 derivs[num_vtx][2] = corner_vals[i][2]; 01098 ++num_vtx; 01099 } 01100 } 01101 01102 void TetLagrangeShape::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out, 01103 MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const 01104 { 01105 if( !nodeset.have_any_mid_node() ) 01106 { 01107 num_vtx = 4; 01108 get_linear_derivatives( vertex_indices_out, d_coeff_d_xi_out ); 01109 return; 01110 } 01111 01112 if( nodeset.have_any_mid_face_node() | nodeset.have_any_mid_region_node() ) 01113 { 01114 MSQ_SETERR( err ) 01115 ( "TetLagrangeShape does not support mid-face/mid-element nodes", MsqError::UNSUPPORTED_ELEMENT ); 01116 return; 01117 } 01118 01119 switch( loc.dimension ) 01120 { 01121 case 0: 01122 derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 01123 break; 01124 case 1: 01125 derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 01126 break; 01127 case 2: 01128 derivatives_at_mid_face( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 01129 break; 01130 case 3: 01131 derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 01132 break; 01133 default: 01134 MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG ); 01135 } 01136 } 01137 01138 void TetLagrangeShape::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const 01139 { 01140 const double a = 1.122462048309373; // 6th root of 2 01141 const double b = 1.7320508075688772; // sqrt(3) 01142 const double c = 1.5874010519681994; // 2 to the 2/3 01143 J( 0, 0 ) = a; 01144 J( 0, 1 ) = 0.5 * a; 01145 J( 0, 2 ) = 0.5 * a; 01146 J( 1, 0 ) = 0.0; 01147 J( 1, 1 ) = 0.5 * a * b; 01148 J( 1, 2 ) = 0.5 * a / b; 01149 J( 2, 0 ) = 0.0; 01150 J( 2, 1 ) = 0.0; 01151 J( 2, 2 ) = c / b; 01152 } 01153 01154 } // namespace MBMesquite