MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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