MOAB: Mesh Oriented datABase  (version 5.4.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     [email protected], [email protected], [email protected],
00024     [email protected], [email protected], [email protected]
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 )
00224             {
00225                 tem *= -1;
00226             }
00227 
00228             return tem;
00229 
00230         case QUADRILATERAL:
00231             cross_vec = ( ( verts[vertexIndices[1]] - verts[vertexIndices[0]] ) *
00232                           ( verts[vertexIndices[3]] - verts[vertexIndices[0]] ) );
00233             pd.get_domain_normal_at_element( element_index, surface_normal, err );
00234             MSQ_ERRZERO( err );
00235             tem = ( cross_vec.length() / 2.0 );
00236             // if normals do not point in same general direction, negate area
00237             if( cross_vec % surface_normal < 0 )
00238             {
00239                 tem *= -1;
00240             }
00241             cross_vec = ( ( verts[vertexIndices[3]] - verts[vertexIndices[2]] ) *
00242                           ( verts[vertexIndices[1]] - verts[vertexIndices[2]] ) );
00243             tem2      = ( cross_vec.length() / 2.0 );
00244             // if normals do not point in same general direction, negate area
00245             if( cross_vec % surface_normal < 0 )
00246             {
00247                 tem2 *= -1;
00248                 // test to make sure surface normal existed
00249                 // if(surface_normal.length_squared()<.5){
00250                 // err.set_msg("compute_signed_area called without surface_normal available.");
00251                 //}
00252             }
00253             return ( tem + tem2 );
00254 
00255         default:
00256             MSQ_SETERR( err )
00257             ( "Invalid type of element passed to compute unsigned area.", MsqError::UNSUPPORTED_ELEMENT );
00258             return 0;
00259     };
00260     return 0.0;
00261 }
00262 
00263 /*!Appends the indices (in the vertex array) of the vertices to connected
00264   to vertex_array[vertex_index] to the end of the vector vert_indices.
00265   The connected vertices are right-hand ordered as defined by the
00266   entity.
00267 
00268 */
00269 void MsqMeshEntity::get_connected_vertices( std::size_t vertex_index,
00270                                             std::vector< std::size_t >& vert_indices,
00271                                             MsqError& err )
00272 {
00273     // index is set to the index in the vertexIndices corresponding
00274     // to vertex_index
00275     int index;
00276     for( index = vertex_count() - 1; index >= 0; --index )
00277         if( vertexIndices[index] == vertex_index ) break;
00278     if( index < 0 )
00279     {
00280         MSQ_SETERR( err )( "Invalid vertex index.", MsqError::INVALID_ARG );
00281         return;
00282     }
00283 
00284     unsigned n;
00285     const unsigned* indices = TopologyInfo::adjacent_vertices( mType, index, n );
00286     if( !indices ) MSQ_SETERR( err )( "Element type not available", MsqError::INVALID_ARG );
00287     for( unsigned i = 0; i < n; ++i )
00288         vert_indices.push_back( vertexIndices[indices[i]] );
00289 }
00290 
00291 /*! Gives the normal at the surface point corner_pt ... but if not available,
00292     gives the normalized cross product of corner_vec1 and corner_vec2.
00293   */
00294 /*
00295 void MsqMeshEntity::compute_corner_normal(size_t corner,
00296                                           Vector3D &normal,
00297                                           PatchData &pd,
00298                                           MsqError &err)
00299 {
00300   if ( get_element_type()==TRIANGLE || get_element_type()==QUADRILATERAL )
00301   {
00302       // There are two cases where we cannot get a normal from the
00303       // geometry that are not errors:
00304       // 1) There is no domain set
00305       // 2) The vertex is at a degenerate point on the geometry (e.g.
00306       //     tip of a cone.)
00307 
00308     bool have_normal = false;
00309 
00310       // Get normal from domain
00311     if (pd.domain_set())
00312     {
00313       size_t index = pd.get_element_index(this);
00314       pd.get_domain_normal_at_corner( index, corner, normal, err );MSQ_ERRRTN(err);
00315 
00316       double length = normal.length();
00317       if (length > DBL_EPSILON)
00318       {
00319         have_normal = true;
00320         normal /= length;
00321       }
00322     }
00323 
00324       // If failed to get normal from domain, calculate
00325       // from adjacent vertices.
00326     if ( !have_normal )
00327     {
00328       const int num_corner = this->vertex_count();
00329       const int prev_corner = (corner + num_corner - 1) % num_corner;
00330       const int next_corner = (corner + 1 ) % num_corner;
00331       const size_t this_idx = vertexIndices[corner];
00332       const size_t prev_idx = vertexIndices[prev_corner];
00333       const size_t next_idx = vertexIndices[next_corner];
00334       normal = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
00335              * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
00336       normal.normalize();
00337     }
00338   }
00339   else
00340     MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
00341                     MsqError::INVALID_ARG);
00342 }
00343 */
00344 
00345 void MsqMeshEntity::compute_corner_normals( Vector3D normals[], PatchData& pd, MsqError& err )
00346 {
00347     EntityTopology type = get_element_type();
00348     if( type != TRIANGLE && type != QUADRILATERAL && type != POLYGON )
00349     {
00350         MSQ_SETERR( err )
00351         ( "Should only be used for faces (tri, quads, ...).", MsqError::INVALID_ARG );
00352         return;
00353     }
00354 
00355     // There are two cases where we cannot get a normal from the
00356     // geometry that are not errors:
00357     // 1) There is no domain set
00358     // 2) The vertex is at a degenerate point on the geometry (e.g.
00359     //     tip of a cone.)
00360 
00361     // Get normal from domain
00362     if( pd.domain_set() )
00363     {
00364         size_t index = pd.get_element_index( this );
00365         pd.get_domain_normals_at_corners( index, normals, err );MSQ_ERRRTN( err );
00366     }
00367 
00368     // Check if normals are valid (none are valid if !pd.domain_set())
00369     const unsigned count = vertex_count();
00370     for( unsigned i = 0; i < count; ++i )
00371     {
00372         // If got valid normal from domain,
00373         // make it a unit vector and continue.
00374         if( pd.domain_set() )
00375         {
00376             double length = normals[i].length();
00377             if( length > DBL_EPSILON )
00378             {
00379                 normals[i] /= length;
00380                 continue;
00381             }
00382         }
00383 
00384         const size_t prev_idx = vertexIndices[( i + count - 1 ) % count];
00385         const size_t this_idx = vertexIndices[i];
00386         const size_t next_idx = vertexIndices[( i + 1 ) % count];
00387 
00388         // Calculate normal using edges adjacent to corner
00389         normals[i] = ( pd.vertex_by_index( next_idx ) - pd.vertex_by_index( this_idx ) ) *
00390                      ( pd.vertex_by_index( prev_idx ) - pd.vertex_by_index( this_idx ) );
00391         normals[i].normalize();
00392     }
00393 }
00394 
00395 ostream& operator<<( ostream& stream, const MsqMeshEntity& entity )
00396 {
00397     size_t num_vert = entity.vertex_count();
00398     stream << "MsqMeshEntity " << &entity << " with vertices ";
00399     for( size_t i = 0; i < num_vert; ++i )
00400         stream << entity.vertexIndices[i] << " ";
00401     stream << endl;
00402     return stream;
00403 }
00404 
00405 /*! Get a array of indices that specifies for the given vertex
00406   the correct matrix map.  This is used by the I_DFT point
00407   relaxation methods in the laplacian smoothers.
00408 
00409 */
00410 size_t MsqMeshEntity::get_local_matrix_map_about_vertex( PatchData& pd,
00411                                                          MsqVertex* vert,
00412                                                          size_t local_map_size,
00413                                                          int* local_map,
00414                                                          MsqError& err ) const
00415 {
00416     // i iterates through elem's vertices
00417     int i = 0;
00418     // index is set to the index in the vertexIndices corresponding
00419     // to vertex_index
00420     int index                     = -1;
00421     int return_val                = 0;
00422     const MsqVertex* vertex_array = pd.get_vertex_array( err );
00423     if( err ) return return_val;
00424 
00425     switch( mType )
00426     {
00427         case TRIANGLE:
00428             MSQ_SETERR( err )
00429             ( "Requested function not yet supported for Triangles.", MsqError::NOT_IMPLEMENTED );
00430 
00431             break;
00432 
00433         case QUADRILATERAL:
00434             MSQ_SETERR( err )
00435             ( "Requested function not yet supported for Quadrilaterals.", MsqError::NOT_IMPLEMENTED );
00436 
00437             break;
00438 
00439         case TETRAHEDRON:
00440             if( local_map_size < 4 )
00441             {
00442                 MSQ_SETERR( err )
00443                 ( "Array of incorrect length sent to function.", MsqError::INVALID_ARG );
00444                 return return_val;
00445             }
00446             return_val = 4;
00447             while( i < 4 )
00448             {
00449                 if( &vertex_array[vertexIndices[i]] == vert )
00450                 {
00451                     index = i;
00452                     break;
00453                 }
00454                 ++i;
00455             }
00456             switch( index )
00457             {
00458                 case( 0 ):
00459                     local_map[0] = 0;
00460                     local_map[1] = 1;
00461                     local_map[2] = 2;
00462                     local_map[3] = 3;
00463                     break;
00464                 case( 1 ):
00465                     local_map[0] = 1;
00466                     local_map[1] = 0;
00467                     local_map[2] = 3;
00468                     local_map[3] = 2;
00469                     break;
00470                 case( 2 ):
00471                     local_map[0] = 2;
00472                     local_map[1] = 3;
00473                     local_map[2] = 0;
00474                     local_map[3] = 1;
00475                     break;
00476                 case( 3 ):
00477                     local_map[0] = 3;
00478                     local_map[1] = 2;
00479                     local_map[2] = 1;
00480                     local_map[3] = 0;
00481                     break;
00482                 default:
00483                     local_map[0] = -1;
00484                     local_map[1] = -1;
00485                     local_map[2] = -1;
00486                     local_map[3] = -1;
00487             };
00488 
00489             break;
00490 
00491         case HEXAHEDRON:
00492             MSQ_SETERR( err )
00493             ( "Requested function not yet supported for Hexahedrons.", MsqError::NOT_IMPLEMENTED );
00494 
00495             break;
00496         default:
00497             MSQ_SETERR( err )( "Element type not available", MsqError::UNSUPPORTED_ELEMENT );
00498             break;
00499     }
00500     return return_val;
00501 }
00502 
00503 void MsqMeshEntity::check_element_orientation( PatchData& pd, int& inverted, int& total, MsqError& err )
00504 {
00505     NodeSet all = all_nodes( err );MSQ_ERRRTN( err );
00506     unsigned i;
00507 
00508     if( TopologyInfo::dimension( mType ) == 2 )
00509     {
00510         if( !pd.domain_set() )
00511         {
00512             total    = 0;
00513             inverted = 0;
00514             return;
00515         }
00516 
00517         const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
00518         if( !mf )
00519         {
00520             check_element_orientation_corners( pd, inverted, total, err );
00521             return;
00522         }
00523 
00524         NodeSet sample = mf->sample_points( all );
00525         total          = sample.num_nodes();
00526         inverted       = 0;
00527 
00528         if( sample.have_any_corner_node() )
00529         {
00530             for( i = 0; i < TopologyInfo::corners( mType ); ++i )
00531                 if( sample.corner_node( i ) ) inverted += inverted_jacobian_2d( pd, all, Sample( 0, i ), err );
00532         }
00533         if( sample.have_any_mid_edge_node() )
00534         {
00535             for( i = 0; i < TopologyInfo::edges( mType ); ++i )
00536                 if( sample.mid_edge_node( i ) ) inverted += inverted_jacobian_2d( pd, all, Sample( 1, i ), err );
00537         }
00538         if( sample.have_any_mid_face_node() ) inverted += inverted_jacobian_2d( pd, all, Sample( 2, 0 ), err );
00539     }
00540     else
00541     {
00542         const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
00543         if( !mf )
00544         {
00545             check_element_orientation_corners( pd, inverted, total, err );
00546             return;
00547         }
00548 
00549         NodeSet sample = mf->sample_points( all );
00550         total          = sample.num_nodes();
00551         inverted       = 0;
00552 
00553         if( sample.have_any_corner_node() )
00554         {
00555             for( i = 0; i < TopologyInfo::corners( mType ); ++i )
00556                 if( sample.corner_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 0, i ), err );
00557         }
00558         if( sample.have_any_mid_edge_node() )
00559         {
00560             for( i = 0; i < TopologyInfo::edges( mType ); ++i )
00561                 if( sample.mid_edge_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 1, i ), err );
00562         }
00563         if( sample.have_any_mid_face_node() )
00564         {
00565             for( i = 0; i < TopologyInfo::faces( mType ); ++i )
00566                 if( sample.mid_face_node( i ) ) inverted += inverted_jacobian_3d( pd, all, Sample( 2, i ), err );
00567         }
00568         if( sample.have_any_mid_region_node() )
00569         {
00570             inverted += inverted_jacobian_3d( pd, all, Sample( 3, 0 ), err );
00571         }
00572     }
00573 }
00574 
00575 bool MsqMeshEntity::inverted_jacobian_3d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
00576 {
00577     MsqMatrix< 3, 3 > J;
00578     MsqVector< 3 > junk[27];
00579     size_t junk2[27], junk3;
00580     assert( node_count() <= 27 );
00581 
00582     const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
00583     mf->jacobian( pd, pd.get_element_index( this ), nodes, sample, junk2, junk, junk3, J, err );
00584     MSQ_ERRZERO( err );
00585     // const double size_eps_sqr = sqr_Frobenius( J ) * DBL_EPSILON;
00586     const double d = det( J );
00587     double l1      = J.column( 0 ) % J.column( 0 );
00588     double l2      = J.column( 1 ) % J.column( 1 );
00589     double l3      = J.column( 2 ) % J.column( 2 );
00590     return d < 0 || d * d < DBL_EPSILON * DBL_EPSILON * l1 * l2 * l3;
00591 }
00592 
00593 bool MsqMeshEntity::inverted_jacobian_2d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
00594 {
00595     MsqMatrix< 3, 2 > J;
00596     MsqVector< 2 > junk[9];
00597     size_t junk2[9], junk3;
00598     assert( node_count() <= 9 );
00599 
00600     const int idx               = pd.get_element_index( this );
00601     const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
00602     mf->jacobian( pd, idx, nodes, sample, junk2, junk, junk3, J, err );
00603     MSQ_ERRZERO( err );
00604     const MsqVector< 3 > cross = J.column( 0 ) * J.column( 1 );
00605 
00606     if( pd.domain_set() )
00607     {
00608         Vector3D norm;
00609         pd.get_domain_normal_at_sample( pd.get_element_index( this ), sample, norm, err );
00610         MSQ_ERRZERO( err );
00611 
00612         const MsqVector< 3 > N( &norm[0] );
00613         if( cross % N < 0.0 ) return true;
00614     }
00615 
00616     const double l1 = J.column( 0 ) % J.column( 0 );
00617     const double l2 = J.column( 1 ) % J.column( 1 );
00618     return cross % cross < DBL_EPSILON * DBL_EPSILON * l1 * l2;
00619 }
00620 
00621 NodeSet MsqMeshEntity::all_nodes( MsqError& err ) const
00622 {
00623     bool mid_edge, mid_face, mid_vol;
00624     TopologyInfo::higher_order( mType, node_count(), mid_edge, mid_face, mid_vol, err );
00625     NodeSet result;
00626     result.set_all_corner_nodes( mType );
00627     if( mid_edge ) result.set_all_mid_edge_nodes( mType );
00628     if( mid_face ) result.set_all_mid_face_nodes( mType );
00629     if( mid_vol ) result.set_all_mid_region_nodes( mType );
00630     return result;
00631 }
00632 
00633 void MsqMeshEntity::check_element_orientation_corners( PatchData& pd, int& inverted, int& total, MsqError& err )
00634 {
00635     int num_nodes = node_count();
00636     total = inverted = 0;
00637 
00638     if( node_count() > vertex_count() )
00639     {
00640         MSQ_SETERR( err )
00641         ( "Cannot perform inversion test for higher-order element"
00642           " without mapping function.",
00643           MsqError::UNSUPPORTED_ELEMENT );
00644         return;
00645     }
00646 
00647     const MsqVertex* vertices = pd.get_vertex_array( err );MSQ_ERRRTN( err );
00648 
00649     const Vector3D d_con( 1.0, 1.0, 1.0 );
00650 
00651     int i;
00652     Vector3D coord_vectors[3];
00653     Vector3D center_vector;
00654 
00655     switch( mType )
00656     {
00657         case TRIANGLE:
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             center_vector    = vertices[vertexIndices[0]];
00664             coord_vectors[0] = vertices[vertexIndices[1]] - center_vector;
00665             coord_vectors[1] = vertices[vertexIndices[2]] - center_vector;
00666             total            = 1;
00667             inverted         = ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00668             break;
00669 
00670         case QUADRILATERAL:
00671 
00672             if( !pd.domain_set() ) return;
00673 
00674             pd.get_domain_normal_at_element( this, coord_vectors[2], err );MSQ_ERRRTN( err );
00675             coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();  // Need unit normal
00676 
00677             for( i = 0; i < 4; ++i )
00678             {
00679                 center_vector    = vertices[vertexIndices[i]];
00680                 coord_vectors[0] = vertices[vertexIndices[( i + 1 ) % 4]] - center_vector;
00681                 coord_vectors[1] = vertices[vertexIndices[( i + 3 ) % 4]] - center_vector;
00682                 ++total;
00683                 inverted += ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00684             }
00685             break;
00686 
00687         case TETRAHEDRON:
00688             center_vector    = vertices[vertexIndices[0]];
00689             coord_vectors[0] = vertices[vertexIndices[1]] - center_vector;
00690             coord_vectors[1] = vertices[vertexIndices[2]] - center_vector;
00691             coord_vectors[2] = vertices[vertexIndices[3]] - center_vector;
00692             total            = 1;
00693             inverted         = ( coord_vectors[0] % ( coord_vectors[1] * coord_vectors[2] ) <= 0.0 );
00694             break;
00695 
00696         case POLYGON:
00697 
00698             if( !pd.domain_set() ) return;
00699 
00700             pd.get_domain_normal_at_element( this, coord_vectors[2], err );MSQ_ERRRTN( err );
00701             coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();  // Need unit normal
00702 
00703             for( i = 0; i < num_nodes; ++i )
00704             {
00705                 center_vector    = vertices[vertexIndices[i]];
00706                 coord_vectors[0] = vertices[vertexIndices[( i + 1 ) % num_nodes]] - center_vector;
00707                 coord_vectors[1] = vertices[vertexIndices[( i + num_nodes - 1 ) % num_nodes]] - center_vector;
00708                 ++total;
00709                 inverted += ( coord_vectors[2] % ( coord_vectors[0] * coord_vectors[1] ) <= 0.0 );
00710             }
00711             break;
00712 
00713         default:  // generic code for 3D elements
00714         {
00715             size_t num_corners = corner_count();
00716             unsigned num_adj;
00717             const unsigned* adj_idx;
00718             for( unsigned j = 0; j < num_corners; ++j )
00719             {
00720                 adj_idx = TopologyInfo::adjacent_vertices( mType, j, num_adj );
00721                 if( 3 != num_adj )
00722                 {
00723                     MSQ_SETERR( err )( "Unsupported element type.", MsqError::INVALID_ARG );
00724                     return;
00725                 }
00726 
00727                 center_vector    = vertices[vertexIndices[j]];
00728                 coord_vectors[0] = vertices[vertexIndices[adj_idx[0]]] - center_vector;
00729                 coord_vectors[1] = vertices[vertexIndices[adj_idx[1]]] - center_vector;
00730                 coord_vectors[2] = vertices[vertexIndices[adj_idx[2]]] - center_vector;
00731                 ++total;
00732                 inverted += ( coord_vectors[0] % ( coord_vectors[1] * coord_vectors[2] ) <= 0.0 );
00733             }
00734             break;
00735         }
00736     }  // end switch over element type
00737 }
00738 
00739 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines