MOAB: Mesh Oriented datABase  (version 5.2.1)
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 <assert.h>
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines