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