MOAB: Mesh Oriented datABase  (version 5.3.1)
MsqMeshEntity.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in 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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 //
00028 // ORIG-DATE: 16-May-02 at 10:26:21
00029 //  LAST-MOD: 18-Jun-04 at 11:36:07 by Thomas Leurent
00030 //
00031 /*! \file MsqMeshEntity.cpp
00032 
00033 \brief All elements in Mesquite are of type MsqMeshEntity. Their associated
00034 functionality is implemented in this file.
00035 
00036     \author Thomas Leurent
00037     \author Michael Brewer
00038     \author Darryl Melander
00039     \date 2002-05-16
00040  */
00041 
00042 #include "Mesquite.hpp"
00043 #include "MsqMeshEntity.hpp"
00044 #include "MsqVertex.hpp"
00045 #include "PatchData.hpp"
00046 
00047 #include <vector>
00048 #include <ostream>
00049 using std::endl;
00050 using std::ostream;
00051 using std::vector;
00052 
00053 namespace MBMesquite
00054 {
00055 
00056 //! Gets the indices of the vertices of this element.
00057 //! The indices are only valid in the PatchData from which
00058 //! this element was retrieved.
00059 //! The order of the vertices is the canonical order for this
00060 //! element's type.
00061 void MsqMeshEntity::get_vertex_indices( std::vector< std::size_t >& vertices ) const
00062 {
00063     vertices.resize( vertex_count() );
00064     std::copy( vertexIndices, vertexIndices + vertex_count(), vertices.begin() );
00065 }
00066 
00067 //! Gets the indices of the vertices of this element.
00068 //! The indices are only valid in the PatchData from which
00069 //! this element was retrieved.
00070 //! The order of the vertices is the canonical order for this
00071 //! element's type.
00072 //! The indices are placed appended to the end of the list.
00073 //! The list is not cleared before appending this entity's vertices.
00074 void MsqMeshEntity::append_vertex_indices( std::vector< std::size_t >& vertex_list ) const
00075 {
00076     vertex_list.insert( vertex_list.end(), vertexIndices, vertexIndices + vertex_count() );
00077 }
00078 
00079 void MsqMeshEntity::get_node_indices( std::vector< std::size_t >& indices ) const
00080 {
00081     indices.resize( node_count() );
00082     std::copy( vertexIndices, vertexIndices + node_count(), indices.begin() );
00083 }
00084 
00085 void MsqMeshEntity::append_node_indices( std::vector< std::size_t >& indices ) const
00086 {
00087     indices.insert( indices.end(), vertexIndices, vertexIndices + node_count() );
00088 }
00089 
00090 /*! The centroid of an element containing n vertices with equal masses is located at
00091   \f[ \b{x} = \frac{ \sum_{i=1}^{n} \b{x}_i }{ n }  \f]
00092   where \f$ \b{x}_i  ,\, i=1,...,n\f$ are the vertices coordinates.
00093 */
00094 void MsqMeshEntity::get_centroid( Vector3D& centroid, const PatchData& pd, MsqError& err ) const
00095 {
00096     centroid               = 0.0;
00097     const MsqVertex* vtces = pd.get_vertex_array( err );MSQ_ERRRTN( err );
00098     size_t nve = vertex_count();
00099     for( size_t i = 0; i < nve; ++i )
00100         centroid += vtces[vertexIndices[i]];
00101     centroid /= nve;
00102 }
00103 
00104 static inline double corner_volume( const Vector3D& v0, const Vector3D& v1, const Vector3D& v2, const Vector3D& v3 )
00105 {
00106     return ( v1 - v0 ) * ( v2 - v0 ) % ( v3 - v0 );
00107 }
00108 
00109 /*!
00110   \brief Computes the area of the given element.  Returned value is
00111   always non-negative.  If the entity passed is not a two-dimensional
00112   element, an error is set.*/
00113 double MsqMeshEntity::compute_unsigned_area( PatchData& pd, MsqError& err )
00114 {
00115     const MsqVertex* verts = pd.get_vertex_array( err );
00116     MSQ_ERRZERO( err );
00117     double tem = 0.0;
00118     switch( mType )
00119     {
00120 
00121         case TRIANGLE:
00122             tem = ( ( verts[vertexIndices[1]] - verts[vertexIndices[0]] ) *
00123                     ( verts[vertexIndices[2]] - verts[vertexIndices[0]] ) )
00124                       .length();
00125             return 0.5 * tem;
00126 
00127         case QUADRILATERAL:
00128             tem = ( ( verts[vertexIndices[1]] - verts[vertexIndices[0]] ) *
00129                     ( verts[vertexIndices[3]] - verts[vertexIndices[0]] ) )
00130                       .length();
00131             tem += ( ( verts[vertexIndices[3]] - verts[vertexIndices[2]] ) *
00132                      ( verts[vertexIndices[1]] - verts[vertexIndices[2]] ) )
00133                        .length();
00134             return ( tem / 2.0 );
00135 
00136         case POLYGON:
00137             // assume convex
00138             for( unsigned i = 1; i < numVertexIndices - 1; ++i )
00139                 tem += ( ( verts[vertexIndices[i]] - verts[vertexIndices[0]] ) *
00140                          ( verts[vertexIndices[i + 1]] - verts[vertexIndices[0]] ) )
00141                            .length();
00142             return 0.5 * tem;
00143 
00144         case TETRAHEDRON:
00145             return 1.0 / 6.0 *
00146                    fabs( corner_volume( verts[vertexIndices[0]], verts[vertexIndices[1]], verts[vertexIndices[2]],
00147                                         verts[vertexIndices[3]] ) );
00148 
00149         case PYRAMID: {
00150             Vector3D m =
00151                 verts[vertexIndices[0]] + verts[vertexIndices[1]] + verts[vertexIndices[2]] + verts[vertexIndices[3]];
00152             Vector3D t1 = verts[vertexIndices[0]] - verts[vertexIndices[2]];
00153             Vector3D t2 = verts[vertexIndices[1]] - verts[vertexIndices[3]];
00154             tem         = ( ( t1 + t2 ) * ( t1 - t2 ) ) % ( verts[vertexIndices[4]] - 0.25 * m );
00155             return ( 1.0 / 12.0 ) * fabs( tem );
00156         }
00157 
00158         case PRISM: {
00159             tem = corner_volume( verts[vertexIndices[0]], verts[vertexIndices[1]], verts[vertexIndices[2]],
00160                                  verts[vertexIndices[3]] );
00161 
00162             tem += corner_volume( verts[vertexIndices[1]], verts[vertexIndices[2]], verts[vertexIndices[3]],
00163                                   verts[vertexIndices[4]] );
00164 
00165             tem += corner_volume( verts[vertexIndices[2]], verts[vertexIndices[3]], verts[vertexIndices[4]],
00166                                   verts[vertexIndices[5]] );
00167 
00168             return 1.0 / 6.0 * fabs( tem );
00169         }
00170 
00171         case HEXAHEDRON: {
00172 
00173             tem = corner_volume( verts[vertexIndices[1]], verts[vertexIndices[2]], verts[vertexIndices[0]],
00174                                  verts[vertexIndices[5]] );
00175 
00176             tem += corner_volume( verts[vertexIndices[3]], verts[vertexIndices[0]], verts[vertexIndices[2]],
00177                                   verts[vertexIndices[7]] );
00178 
00179             tem += corner_volume( verts[vertexIndices[4]], verts[vertexIndices[7]], verts[vertexIndices[5]],
00180                                   verts[vertexIndices[0]] );
00181 
00182             tem += corner_volume( verts[vertexIndices[6]], verts[vertexIndices[5]], verts[vertexIndices[7]],
00183                                   verts[vertexIndices[2]] );
00184 
00185             tem += corner_volume( verts[vertexIndices[5]], verts[vertexIndices[2]], verts[vertexIndices[0]],
00186                                   verts[vertexIndices[7]] );
00187 
00188             return ( 1.0 / 6.0 ) * fabs( tem );
00189         }
00190 
00191         default:
00192             MSQ_SETERR( err )
00193             ( "Invalid type of element passed to compute unsigned area.", MsqError::UNSUPPORTED_ELEMENT );
00194             return 0;
00195     }
00196     return 0;
00197 }
00198 
00199 /*!
00200   \brief Computes the area of the given element.  Returned value can be
00201   negative.  If the entity passed is not a two-dimensional element, an
00202   error is set.*/
00203 double MsqMeshEntity::compute_signed_area( PatchData& pd, MsqError& err )
00204 {
00205     const MsqVertex* verts = pd.get_vertex_array( err );
00206     MSQ_ERRZERO( err );
00207     double tem  = 0.0;
00208     double tem2 = 0.0;
00209     Vector3D surface_normal;
00210     Vector3D cross_vec;
00211     size_t element_index = pd.get_element_index( this );
00212 
00213     switch( mType )
00214     {
00215 
00216         case TRIANGLE:
00217             cross_vec = ( ( verts[vertexIndices[1]] - verts[vertexIndices[0]] ) *
00218                           ( verts[vertexIndices[2]] - verts[vertexIndices[0]] ) );
00219             pd.get_domain_normal_at_element( element_index, surface_normal, err );
00220             MSQ_ERRZERO( err );
00221             tem = ( cross_vec.length() / 2.0 );
00222             // if normals do not point in same general direction, negate area
00223             if( cross_vec % surface_normal < 0 ) { tem *= -1; }
00224 
00225             return tem;
00226 
00227         case QUADRILATERAL:
00228             cross_vec = ( ( verts[vertexIndices[1]] - verts[vertexIndices[0]] ) *
00229                           ( verts[vertexIndices[3]] - verts[vertexIndices[0]] ) );
00230             pd.get_domain_normal_at_element( element_index, surface_normal, err );
00231             MSQ_ERRZERO( err );
00232             tem = ( cross_vec.length() / 2.0 );
00233             // if normals do not point in same general direction, negate area
00234             if( cross_vec % surface_normal < 0 ) { tem *= -1; }
00235             cross_vec = ( ( verts[vertexIndices[3]] - verts[vertexIndices[2]] ) *
00236                           ( verts[vertexIndices[1]] - verts[vertexIndices[2]] ) );
00237             tem2      = ( cross_vec.length() / 2.0 );
00238             // if normals do not point in same general direction, negate area
00239             if( cross_vec % surface_normal < 0 )
00240             {
00241                 tem2 *= -1;
00242                 // test to make sure surface normal existed
00243                 // if(surface_normal.length_squared()<.5){
00244                 // err.set_msg("compute_signed_area called without surface_normal available.");
00245                 //}
00246             }
00247             return ( tem + tem2 );
00248 
00249         default:
00250             MSQ_SETERR( err )
00251             ( "Invalid type of element passed to compute unsigned area.", MsqError::UNSUPPORTED_ELEMENT );
00252             return 0;
00253     };
00254     return 0.0;
00255 }
00256 
00257 /*!Appends the indices (in the vertex array) of the vertices to connected
00258   to vertex_array[vertex_index] to the end of the vector vert_indices.
00259   The connected vertices are right-hand ordered as defined by the
00260   entity.
00261 
00262 */
00263 void MsqMeshEntity::get_connected_vertices( std::size_t vertex_index, std::vector< std::size_t >& vert_indices,
00264                                             MsqError& err )
00265 {
00266     // index is set to the index in the vertexIndices corresponding
00267     // to vertex_index
00268     int index;
00269     for( index = vertex_count() - 1; index >= 0; --index )
00270         if( vertexIndices[index] == vertex_index ) break;
00271     if( index < 0 )
00272     {
00273         MSQ_SETERR( err )( "Invalid vertex index.", MsqError::INVALID_ARG );
00274         return;
00275     }
00276 
00277     unsigned n;
00278     const unsigned* indices = TopologyInfo::adjacent_vertices( mType, index, n );
00279     if( !indices ) MSQ_SETERR( err )( "Element type not available", MsqError::INVALID_ARG );
00280     for( unsigned i = 0; i < n; ++i )
00281         vert_indices.push_back( vertexIndices[indices[i]] );
00282 }
00283 
00284 /*! Gives the normal at the surface point corner_pt ... but if not available,
00285     gives the normalized cross product of corner_vec1 and corner_vec2.
00286   */
00287 /*
00288 void MsqMeshEntity::compute_corner_normal(size_t corner,
00289                                           Vector3D &normal,
00290                                           PatchData &pd,
00291                                           MsqError &err)
00292 {
00293   if ( get_element_type()==TRIANGLE || get_element_type()==QUADRILATERAL )
00294   {
00295       // There are two cases where we cannot get a normal from the
00296       // geometry that are not errors:
00297       // 1) There is no domain set
00298       // 2) The vertex is at a degenerate point on the geometry (e.g.
00299       //     tip of a cone.)
00300 
00301     bool have_normal = false;
00302 
00303       // Get normal from domain
00304     if (pd.domain_set())
00305     {
00306       size_t index = pd.get_element_index(this);
00307       pd.get_domain_normal_at_corner( index, corner, normal, err );MSQ_ERRRTN(err);
00308 
00309       double length = normal.length();
00310       if (length > DBL_EPSILON)
00311       {
00312         have_normal = true;
00313         normal /= length;
00314       }
00315     }
00316 
00317       // If failed to get normal from domain, calculate
00318       // from adjacent vertices.
00319     if ( !have_normal )
00320     {
00321       const int num_corner = this->vertex_count();
00322       const int prev_corner = (corner + num_corner - 1) % num_corner;
00323       const int next_corner = (corner + 1 ) % num_corner;
00324       const size_t this_idx = vertexIndices[corner];
00325       const size_t prev_idx = vertexIndices[prev_corner];
00326       const size_t next_idx = vertexIndices[next_corner];
00327       normal = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
00328              * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
00329       normal.normalize();
00330     }
00331   }
00332   else
00333     MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
00334                     MsqError::INVALID_ARG);
00335 }
00336 */
00337 
00338 void MsqMeshEntity::compute_corner_normals( Vector3D normals[], PatchData& pd, MsqError& err )
00339 {
00340     EntityTopology type = get_element_type();
00341     if( type != TRIANGLE && type != QUADRILATERAL && type != POLYGON )
00342     {
00343         MSQ_SETERR( err )
00344         ( "Should only be used for faces (tri, quads, ...).", MsqError::INVALID_ARG );
00345         return;
00346     }
00347 
00348     // There are two cases where we cannot get a normal from the
00349     // geometry that are not errors:
00350     // 1) There is no domain set
00351     // 2) The vertex is at a degenerate point on the geometry (e.g.
00352     //     tip of a cone.)
00353 
00354     // Get normal from domain
00355     if( pd.domain_set() )
00356     {
00357         size_t index = pd.get_element_index( this );
00358         pd.get_domain_normals_at_corners( index, normals, err );MSQ_ERRRTN( err );
00359     }
00360 
00361     // Check if normals are valid (none are valid if !pd.domain_set())
00362     const unsigned count = vertex_count();
00363     for( unsigned i = 0; i < count; ++i )
00364     {
00365         // If got valid normal from domain,
00366         // make it a unit vector and continue.
00367         if( pd.domain_set() )
00368         {
00369             double length = normals[i].length();
00370             if( length > DBL_EPSILON )
00371             {
00372                 normals[i] /= length;
00373                 continue;
00374             }
00375         }
00376 
00377         const size_t prev_idx = vertexIndices[( i + count - 1 ) % count];
00378         const size_t this_idx = vertexIndices[i];
00379         const size_t next_idx = vertexIndices[( i + 1 ) % count];
00380 
00381         // Calculate normal using edges adjacent to corner
00382         normals[i] = ( pd.vertex_by_index( next_idx ) - pd.vertex_by_index( this_idx ) ) *
00383                      ( pd.vertex_by_index( prev_idx ) - pd.vertex_by_index( this_idx ) );
00384         normals[i].normalize();
00385     }
00386 }
00387 
00388 ostream& operator<<( ostream& stream, const MsqMeshEntity& entity )
00389 {
00390     size_t num_vert = entity.vertex_count();
00391     stream << "MsqMeshEntity " << &entity << " with vertices ";
00392     for( size_t i = 0; i < num_vert; ++i )
00393         stream << entity.vertexIndices[i] << " ";
00394     stream << endl;
00395     return stream;
00396 }
00397 
00398 /*! Get a array of indices that specifies for the given vertex
00399   the correct matrix map.  This is used by the I_DFT point
00400   relaxation methods in the laplacian smoothers.
00401 
00402 */
00403 size_t MsqMeshEntity::get_local_matrix_map_about_vertex( PatchData& pd, MsqVertex* vert, size_t local_map_size,
00404                                                          int* local_map, MsqError& err ) const
00405 {
00406     // i iterates through elem's vertices
00407     int i = 0;
00408     // index is set to the index in the vertexIndices corresponding
00409     // to vertex_index
00410     int index                     = -1;
00411     int return_val                = 0;
00412     const MsqVertex* vertex_array = pd.get_vertex_array( err );
00413     if( err ) return return_val;
00414 
00415     switch( mType )
00416     {
00417         case TRIANGLE:
00418             MSQ_SETERR( err )
00419             ( "Requested function not yet supported for Triangles.", MsqError::NOT_IMPLEMENTED );
00420 
00421             break;
00422 
00423         case QUADRILATERAL:
00424             MSQ_SETERR( err )
00425             ( "Requested function not yet supported for Quadrilaterals.", MsqError::NOT_IMPLEMENTED );
00426 
00427             break;
00428 
00429         case TETRAHEDRON:
00430             if( local_map_size < 4 )
00431             {
00432                 MSQ_SETERR( err )
00433                 ( "Array of incorrect length sent to function.", MsqError::INVALID_ARG );
00434                 return return_val;
00435             }
00436             return_val = 4;
00437             while( i < 4 )
00438             {
00439                 if( &vertex_array[vertexIndices[i]] == vert )
00440                 {
00441                     index = i;
00442                     break;
00443                 }
00444                 ++i;
00445             }
00446             switch( index )
00447             {
00448                 case( 0 ):
00449                     local_map[0] = 0;
00450                     local_map[1] = 1;
00451                     local_map[2] = 2;
00452                     local_map[3] = 3;
00453                     break;
00454                 case( 1 ):
00455                     local_map[0] = 1;
00456                     local_map[1] = 0;
00457                     local_map[2] = 3;
00458                     local_map[3] = 2;
00459                     break;
00460                 case( 2 ):
00461                     local_map[0] = 2;
00462                     local_map[1] = 3;
00463                     local_map[2] = 0;
00464                     local_map[3] = 1;
00465                     break;
00466                 case( 3 ):
00467                     local_map[0] = 3;
00468                     local_map[1] = 2;
00469                     local_map[2] = 1;
00470                     local_map[3] = 0;
00471                     break;
00472                 default:
00473                     local_map[0] = -1;
00474                     local_map[1] = -1;
00475                     local_map[2] = -1;
00476                     local_map[3] = -1;
00477             };
00478 
00479             break;
00480 
00481         case HEXAHEDRON:
00482             MSQ_SETERR( err )
00483             ( "Requested function not yet supported for Hexahedrons.", MsqError::NOT_IMPLEMENTED );
00484 
00485             break;
00486         default:
00487             MSQ_SETERR( err )( "Element type not available", MsqError::UNSUPPORTED_ELEMENT );
00488             break;
00489     }
00490     return return_val;
00491 }
00492 
00493 void MsqMeshEntity::check_element_orientation( PatchData& pd, int& inverted, int& total, MsqError& err )
00494 {
00495     NodeSet all = all_nodes( err );MSQ_ERRRTN( err );
00496     unsigned i;
00497 
00498     if( TopologyInfo::dimension( mType ) == 2 )
00499     {
00500         if( !pd.domain_set() )
00501         {
00502             total    = 0;
00503             inverted = 0;
00504             return;
00505         }
00506 
00507         const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
00508         if( !mf )
00509         {
00510             check_element_orientation_corners( pd, inverted, total, err );
00511             return;
00512         }
00513 
00514         NodeSet sample = mf->sample_points( all );
00515         total          = sample.num_nodes();
00516         inverted       = 0;
00517 
00518         if( sample.have_any_corner_node() )
00519         {
00520             for( i = 0; i < TopologyInfo::corners( mType ); ++i )
00521                 if( sample.corner_node( i ) ) inverted += inverted_jacobian_2d( pd, all, Sample( 0, i ), err );
00522         }
00523         if( sample.have_any_mid_edge_node() )
00524         {
00525             for( i = 0; i < TopologyInfo::edges( mType ); ++i )
00526                 if( sample.mid_edge_node( i ) ) inverted += inverted_jacobian_2d( pd, all, Sample( 1, i ), err );
00527         }
00528         if( sample.have_any_mid_face_node() ) inverted += inverted_jacobian_2d( pd, all, Sample( 2, 0 ), err );
00529     }
00530     else
00531     {
00532         const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
00533         if( !mf )
00534         {
00535             check_element_orientation_corners( pd, inverted, total, err );
00536             return;
00537         }
00538 
00539         NodeSet sample = mf->sample_points( all );
00540         total          = sample.num_nodes();
00541         inverted       = 0;
00542 
00543         if( sample.have_any_corner_node() )
00544         {
00545             for( i = 0; i < TopologyInfo::corners( mType ); ++i )
00546                 if( sample.corner_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 0, i ), err );
00547         }
00548         if( sample.have_any_mid_edge_node() )
00549         {
00550             for( i = 0; i < TopologyInfo::edges( mType ); ++i )
00551                 if( sample.mid_edge_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 1, i ), err );
00552         }
00553         if( sample.have_any_mid_face_node() )
00554         {
00555             for( i = 0; i < TopologyInfo::faces( mType ); ++i )
00556                 if( sample.mid_face_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 2, i ), err );
00557         }
00558         if( sample.have_any_mid_region_node() ) { inverted += inverted_jacobian_3d( pd, all, Sample( 3, 0 ), err ); }
00559     }
00560 }
00561 
00562 bool MsqMeshEntity::inverted_jacobian_3d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
00563 {
00564     MsqMatrix< 3, 3 > J;
00565     MsqVector< 3 > junk[27];
00566     size_t junk2[27], junk3;
00567     assert( node_count() <= 27 );
00568 
00569     const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
00570     mf->jacobian( pd, pd.get_element_index( this ), nodes, sample, junk2, junk, junk3, J, err );
00571     MSQ_ERRZERO( err );
00572     // const double size_eps_sqr = sqr_Frobenius( J ) * DBL_EPSILON;
00573     const double d = det( J );
00574     double l1      = J.column( 0 ) % J.column( 0 );
00575     double l2      = J.column( 1 ) % J.column( 1 );
00576     double l3      = J.column( 2 ) % J.column( 2 );
00577     return d < 0 || d * d < DBL_EPSILON * DBL_EPSILON * l1 * l2 * l3;
00578 }
00579 
00580 bool MsqMeshEntity::inverted_jacobian_2d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
00581 {
00582     MsqMatrix< 3, 2 > J;
00583     MsqVector< 2 > junk[9];
00584     size_t junk2[9], junk3;
00585     assert( node_count() <= 9 );
00586 
00587     const int idx               = pd.get_element_index( this );
00588     const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
00589     mf->jacobian( pd, idx, nodes, sample, junk2, junk, junk3, J, err );
00590     MSQ_ERRZERO( err );
00591     const MsqVector< 3 > cross = J.column( 0 ) * J.column( 1 );
00592 
00593     if( pd.domain_set() )
00594     {
00595         Vector3D norm;
00596         pd.get_domain_normal_at_sample( pd.get_element_index( this ), sample, norm, err );
00597         MSQ_ERRZERO( err );
00598 
00599         const MsqVector< 3 > N( &norm[0] );
00600         if( cross % N < 0.0 ) return true;
00601     }
00602 
00603     const double l1 = J.column( 0 ) % J.column( 0 );
00604     const double l2 = J.column( 1 ) % J.column( 1 );
00605     return cross % cross < DBL_EPSILON * DBL_EPSILON * l1 * l2;
00606 }
00607 
00608 NodeSet MsqMeshEntity::all_nodes( MsqError& err ) const
00609 {
00610     bool mid_edge, mid_face, mid_vol;
00611     TopologyInfo::higher_order( mType, node_count(), mid_edge, mid_face, mid_vol, err );
00612     NodeSet result;
00613     result.set_all_corner_nodes( mType );
00614     if( mid_edge ) result.set_all_mid_edge_nodes( mType );
00615     if( mid_face ) result.set_all_mid_face_nodes( mType );
00616     if( mid_vol ) result.set_all_mid_region_nodes( mType );
00617     return result;
00618 }
00619 
00620 void MsqMeshEntity::check_element_orientation_corners( PatchData& pd, int& inverted, int& total, MsqError& err )
00621 {
00622     int num_nodes = node_count();
00623     total = inverted = 0;
00624 
00625     if( node_count() > vertex_count() )
00626     {
00627         MSQ_SETERR( err )
00628         ( "Cannot perform inversion test for higher-order element"
00629           " without mapping function.",
00630           MsqError::UNSUPPORTED_ELEMENT );
00631         return;
00632     }
00633 
00634     const MsqVertex* vertices = pd.get_vertex_array( err );MSQ_ERRRTN( err );
00635 
00636     const Vector3D d_con( 1.0, 1.0, 1.0 );
00637 
00638     int i;
00639     Vector3D coord_vectors[3];
00640     Vector3D center_vector;
00641 
00642     switch( mType )
00643     {
00644         case TRIANGLE:
00645 
00646             if( !pd.domain_set() ) return;
00647 
00648             pd.get_domain_normal_at_element( this, coord_vectors[2], err );MSQ_ERRRTN( err );
00649             coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();  // Need unit normal
00650             center_vector    = vertices[vertexIndices[0]];
00651             coord_vectors[0] = vertices[vertexIndices[1]] - center_vector;
00652             coord_vectors[1] = vertices[vertexIndices[2]] - center_vector;
00653             total            = 1;
00654             inverted         = ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00655             break;
00656 
00657         case QUADRILATERAL:
00658 
00659             if( !pd.domain_set() ) return;
00660 
00661             pd.get_domain_normal_at_element( this, coord_vectors[2], err );MSQ_ERRRTN( err );
00662             coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();  // Need unit normal
00663 
00664             for( i = 0; i < 4; ++i )
00665             {
00666                 center_vector    = vertices[vertexIndices[i]];
00667                 coord_vectors[0] = vertices[vertexIndices[( i + 1 ) % 4]] - center_vector;
00668                 coord_vectors[1] = vertices[vertexIndices[( i + 3 ) % 4]] - center_vector;
00669                 ++total;
00670                 inverted += ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00671             }
00672             break;
00673 
00674         case TETRAHEDRON:
00675             center_vector    = vertices[vertexIndices[0]];
00676             coord_vectors[0] = vertices[vertexIndices[1]] - center_vector;
00677             coord_vectors[1] = vertices[vertexIndices[2]] - center_vector;
00678             coord_vectors[2] = vertices[vertexIndices[3]] - center_vector;
00679             total            = 1;
00680             inverted         = ( coord_vectors[0] % ( coord_vectors[1] * coord_vectors[2] ) <= 0.0 );
00681             break;
00682 
00683         case POLYGON:
00684 
00685             if( !pd.domain_set() ) return;
00686 
00687             pd.get_domain_normal_at_element( this, coord_vectors[2], err );MSQ_ERRRTN( err );
00688             coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();  // Need unit normal
00689 
00690             for( i = 0; i < num_nodes; ++i )
00691             {
00692                 center_vector    = vertices[vertexIndices[i]];
00693                 coord_vectors[0] = vertices[vertexIndices[( i + 1 ) % num_nodes]] - center_vector;
00694                 coord_vectors[1] = vertices[vertexIndices[( i + num_nodes - 1 ) % num_nodes]] - center_vector;
00695                 ++total;
00696                 inverted += ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00697             }
00698             break;
00699 
00700         default:  // generic code for 3D elements
00701         {
00702             size_t num_corners = corner_count();
00703             unsigned num_adj;
00704             const unsigned* adj_idx;
00705             for( unsigned j = 0; j < num_corners; ++j )
00706             {
00707                 adj_idx = TopologyInfo::adjacent_vertices( mType, j, num_adj );
00708                 if( 3 != num_adj )
00709                 {
00710                     MSQ_SETERR( err )( "Unsupported element type.", MsqError::INVALID_ARG );
00711                     return;
00712                 }
00713 
00714                 center_vector    = vertices[vertexIndices[j]];
00715                 coord_vectors[0] = vertices[vertexIndices[adj_idx[0]]] - center_vector;
00716                 coord_vectors[1] = vertices[vertexIndices[adj_idx[1]]] - center_vector;
00717                 coord_vectors[2] = vertices[vertexIndices[adj_idx[2]]] - center_vector;
00718                 ++total;
00719                 inverted += ( coord_vectors[0] % ( coord_vectors[1] * coord_vectors[2] ) <= 0.0 );
00720             }
00721             break;
00722         }
00723     }  // end switch over element type
00724 }
00725 
00726 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines