MOAB: Mesh Oriented datABase
(version 5.4.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2004 Lawrence Livermore National Laboratory. Under 00005 the terms of Contract B545069 with the University of Wisconsin -- 00006 Madison, Lawrence Livermore National Laboratory 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] 00024 00025 ***************************************************************** */ 00026 00027 #include "MsqError.hpp" 00028 #include "TopologyInfo.hpp" 00029 00030 #include <cstring> 00031 #include <cassert> 00032 00033 namespace MBMesquite 00034 { 00035 00036 TopologyInfo TopologyInfo::instance; 00037 00038 const char long_polygon_name[] = "Polygon"; 00039 const char long_triangle_name[] = "Triangle"; 00040 const char long_quadrilateral_name[] = "Quadrilateral"; 00041 const char long_polyhedron_name[] = "Polyhedron"; 00042 const char long_tetrahedron_name[] = "Tetrahedron"; 00043 const char long_hexahedron_name[] = "Hexahedron"; 00044 const char long_prism_name[] = "Prism"; 00045 const char long_pyramid_name[] = "Pyramd"; 00046 const char long_septahedron_name[] = "Septahedron"; 00047 const char short_polygon_name[] = "Polygon"; 00048 const char short_triangle_name[] = "Tri"; 00049 const char short_quadrilateral_name[] = "Quad"; 00050 const char short_polyhedron_name[] = "Polyhedron"; 00051 const char short_tetrahedron_name[] = "Tet"; 00052 const char short_hexahedron_name[] = "Hex"; 00053 const char short_prism_name[] = "Pri"; 00054 const char short_pyramid_name[] = "Pyr"; 00055 const char short_septahedron_name[] = "Sept"; 00056 00057 TopologyInfo::TopologyInfo() 00058 { 00059 memset( dimMap, 0, sizeof( dimMap ) ); 00060 memset( adjMap, 0, sizeof( adjMap ) ); 00061 memset( edgeMap, 0, sizeof( edgeMap ) ); 00062 memset( faceMap, 0, sizeof( faceMap ) ); 00063 memset( vertAdjMap, 0, sizeof( vertAdjMap ) ); 00064 memset( shortNames, 0, sizeof( shortNames ) ); 00065 memset( longNames, 0, sizeof( longNames ) ); 00066 00067 longNames[POLYGON] = long_polygon_name; 00068 longNames[TRIANGLE] = long_triangle_name; 00069 longNames[QUADRILATERAL] = long_quadrilateral_name; 00070 longNames[POLYHEDRON] = long_polyhedron_name; 00071 longNames[TETRAHEDRON] = long_tetrahedron_name; 00072 longNames[HEXAHEDRON] = long_hexahedron_name; 00073 longNames[PRISM] = long_prism_name; 00074 longNames[PYRAMID] = long_pyramid_name; 00075 longNames[SEPTAHEDRON] = long_septahedron_name; 00076 00077 shortNames[POLYGON] = short_polygon_name; 00078 shortNames[TRIANGLE] = short_triangle_name; 00079 shortNames[QUADRILATERAL] = short_quadrilateral_name; 00080 shortNames[POLYHEDRON] = short_polyhedron_name; 00081 shortNames[TETRAHEDRON] = short_tetrahedron_name; 00082 shortNames[HEXAHEDRON] = short_hexahedron_name; 00083 shortNames[PRISM] = short_prism_name; 00084 shortNames[PYRAMID] = short_pyramid_name; 00085 shortNames[SEPTAHEDRON] = short_septahedron_name; 00086 00087 dimMap[POLYGON] = 2; 00088 dimMap[TRIANGLE] = 2; 00089 dimMap[QUADRILATERAL] = 2; 00090 dimMap[POLYHEDRON] = 3; 00091 dimMap[TETRAHEDRON] = 3; 00092 dimMap[HEXAHEDRON] = 3; 00093 dimMap[PRISM] = 3; 00094 dimMap[PYRAMID] = 3; 00095 dimMap[SEPTAHEDRON] = 3; 00096 00097 adjMap[TRIANGLE][0] = 3; 00098 adjMap[TRIANGLE][1] = 3; 00099 adjMap[TRIANGLE][2] = 1; 00100 adjMap[TRIANGLE][3] = 0; 00101 00102 adjMap[QUADRILATERAL][0] = 4; 00103 adjMap[QUADRILATERAL][1] = 4; 00104 adjMap[QUADRILATERAL][2] = 1; 00105 adjMap[QUADRILATERAL][3] = 0; 00106 00107 adjMap[TETRAHEDRON][0] = 4; 00108 adjMap[TETRAHEDRON][1] = 6; 00109 adjMap[TETRAHEDRON][2] = 4; 00110 adjMap[TETRAHEDRON][3] = 1; 00111 00112 adjMap[HEXAHEDRON][0] = 8; 00113 adjMap[HEXAHEDRON][1] = 12; 00114 adjMap[HEXAHEDRON][2] = 6; 00115 adjMap[HEXAHEDRON][3] = 1; 00116 00117 adjMap[PRISM][0] = 6; 00118 adjMap[PRISM][1] = 9; 00119 adjMap[PRISM][2] = 5; 00120 adjMap[PRISM][3] = 1; 00121 00122 adjMap[PYRAMID][0] = 5; 00123 adjMap[PYRAMID][1] = 8; 00124 adjMap[PYRAMID][2] = 5; 00125 adjMap[PYRAMID][3] = 1; 00126 00127 adjMap[SEPTAHEDRON][0] = 7; 00128 adjMap[SEPTAHEDRON][1] = 11; 00129 adjMap[SEPTAHEDRON][2] = 6; /* See description in TSTT mesh interface doc */ 00130 adjMap[SEPTAHEDRON][3] = 1; 00131 00132 int side; 00133 for( side = 0; side < 3; ++side ) 00134 { 00135 edgeMap[TRIANGLE - FIRST_FACE][side][0] = side; 00136 edgeMap[TRIANGLE - FIRST_FACE][side][1] = ( side + 1 ) % 3; 00137 } 00138 for( side = 0; side < 4; ++side ) 00139 { 00140 edgeMap[QUADRILATERAL - FIRST_FACE][side][0] = side; 00141 edgeMap[QUADRILATERAL - FIRST_FACE][side][1] = ( side + 1 ) % 4; 00142 } 00143 for( side = 0; side < 3; ++side ) 00144 { 00145 edgeMap[TETRAHEDRON - FIRST_FACE][side][0] = side; 00146 edgeMap[TETRAHEDRON - FIRST_FACE][side][1] = ( side + 1 ) % 3; 00147 } 00148 for( side = 3; side < 6; ++side ) 00149 { 00150 edgeMap[TETRAHEDRON - FIRST_FACE][side][0] = side - 3; 00151 edgeMap[TETRAHEDRON - FIRST_FACE][side][1] = 3; 00152 } 00153 for( side = 0; side < 4; ++side ) 00154 { 00155 edgeMap[HEXAHEDRON - FIRST_FACE][side][0] = side; 00156 edgeMap[HEXAHEDRON - FIRST_FACE][side][1] = ( side + 1 ) % 4; 00157 } 00158 for( side = 4; side < 8; ++side ) 00159 { 00160 edgeMap[HEXAHEDRON - FIRST_FACE][side][0] = side - 4; 00161 edgeMap[HEXAHEDRON - FIRST_FACE][side][1] = side; 00162 } 00163 for( side = 8; side < 12; ++side ) 00164 { 00165 edgeMap[HEXAHEDRON - FIRST_FACE][side][0] = side - 4; 00166 edgeMap[HEXAHEDRON - FIRST_FACE][side][1] = 4 + ( side + 1 ) % 4; 00167 } 00168 for( side = 0; side < 3; ++side ) 00169 { 00170 edgeMap[PRISM - FIRST_FACE][side][0] = side; 00171 edgeMap[PRISM - FIRST_FACE][side][1] = ( side + 1 ) % 3; 00172 } 00173 for( side = 3; side < 6; ++side ) 00174 { 00175 edgeMap[PRISM - FIRST_FACE][side][0] = side - 3; 00176 edgeMap[PRISM - FIRST_FACE][side][1] = side; 00177 } 00178 for( side = 6; side < 9; ++side ) 00179 { 00180 edgeMap[PRISM - FIRST_FACE][side][0] = side - 3; 00181 edgeMap[PRISM - FIRST_FACE][side][1] = 3 + ( side + 1 ) % 3; 00182 } 00183 for( side = 0; side < 4; ++side ) 00184 { 00185 edgeMap[PYRAMID - FIRST_FACE][side][0] = side; 00186 edgeMap[PYRAMID - FIRST_FACE][side][1] = ( side + 1 ) % 4; 00187 } 00188 for( side = 4; side < 8; ++side ) 00189 { 00190 edgeMap[PYRAMID - FIRST_FACE][side][0] = side - 4; 00191 edgeMap[PYRAMID - FIRST_FACE][side][1] = 4; 00192 } 00193 00194 for( side = 0; side < 3; ++side ) 00195 { 00196 faceMap[TETRAHEDRON - FIRST_VOL][side][0] = 3; 00197 faceMap[TETRAHEDRON - FIRST_VOL][side][1] = side; 00198 faceMap[TETRAHEDRON - FIRST_VOL][side][2] = ( side + 1 ) % 3; 00199 faceMap[TETRAHEDRON - FIRST_VOL][side][3] = 3; 00200 } 00201 faceMap[TETRAHEDRON - FIRST_VOL][3][0] = 3; 00202 faceMap[TETRAHEDRON - FIRST_VOL][3][1] = 2; 00203 faceMap[TETRAHEDRON - FIRST_VOL][3][2] = 1; 00204 faceMap[TETRAHEDRON - FIRST_VOL][3][3] = 0; 00205 00206 for( side = 0; side < 4; ++side ) 00207 { 00208 faceMap[HEXAHEDRON - FIRST_VOL][side][0] = 4; 00209 faceMap[HEXAHEDRON - FIRST_VOL][side][1] = side; 00210 faceMap[HEXAHEDRON - FIRST_VOL][side][2] = ( side + 1 ) % 4; 00211 faceMap[HEXAHEDRON - FIRST_VOL][side][3] = 4 + ( side + 1 ) % 4; 00212 faceMap[HEXAHEDRON - FIRST_VOL][side][4] = side + 4; 00213 } 00214 faceMap[HEXAHEDRON - FIRST_VOL][4][0] = 4; 00215 faceMap[HEXAHEDRON - FIRST_VOL][4][1] = 3; 00216 faceMap[HEXAHEDRON - FIRST_VOL][4][2] = 2; 00217 faceMap[HEXAHEDRON - FIRST_VOL][4][3] = 1; 00218 faceMap[HEXAHEDRON - FIRST_VOL][4][4] = 0; 00219 faceMap[HEXAHEDRON - FIRST_VOL][5][0] = 4; 00220 faceMap[HEXAHEDRON - FIRST_VOL][5][1] = 4; 00221 faceMap[HEXAHEDRON - FIRST_VOL][5][2] = 5; 00222 faceMap[HEXAHEDRON - FIRST_VOL][5][3] = 6; 00223 faceMap[HEXAHEDRON - FIRST_VOL][5][4] = 7; 00224 00225 for( side = 0; side < 4; ++side ) 00226 { 00227 faceMap[PYRAMID - FIRST_VOL][side][0] = 3; 00228 faceMap[PYRAMID - FIRST_VOL][side][1] = side; 00229 faceMap[PYRAMID - FIRST_VOL][side][2] = ( side + 1 ) % 4; 00230 faceMap[PYRAMID - FIRST_VOL][side][3] = 4; 00231 } 00232 faceMap[PYRAMID - FIRST_VOL][4][0] = 4; 00233 faceMap[PYRAMID - FIRST_VOL][4][1] = 3; 00234 faceMap[PYRAMID - FIRST_VOL][4][2] = 2; 00235 faceMap[PYRAMID - FIRST_VOL][4][3] = 1; 00236 faceMap[PYRAMID - FIRST_VOL][4][4] = 0; 00237 00238 for( side = 0; side < 3; ++side ) 00239 { 00240 faceMap[PRISM - FIRST_VOL][side][0] = 4; 00241 faceMap[PRISM - FIRST_VOL][side][1] = side; 00242 faceMap[PRISM - FIRST_VOL][side][2] = ( side + 1 ) % 3; 00243 faceMap[PRISM - FIRST_VOL][side][3] = 3 + ( side + 1 ) % 3; 00244 faceMap[PRISM - FIRST_VOL][side][4] = side + 3; 00245 } 00246 faceMap[PRISM - FIRST_VOL][3][0] = 3; 00247 faceMap[PRISM - FIRST_VOL][3][1] = 2; 00248 faceMap[PRISM - FIRST_VOL][3][2] = 1; 00249 faceMap[PRISM - FIRST_VOL][3][3] = 0; 00250 faceMap[PRISM - FIRST_VOL][4][0] = 3; 00251 faceMap[PRISM - FIRST_VOL][4][1] = 3; 00252 faceMap[PRISM - FIRST_VOL][4][2] = 4; 00253 faceMap[PRISM - FIRST_VOL][4][3] = 5; 00254 00255 int i; 00256 for( i = 0; i < 3; ++i ) 00257 { 00258 vertAdjMap[TRIANGLE - FIRST_FACE][i][0] = 2; 00259 vertAdjMap[TRIANGLE - FIRST_FACE][i][1] = ( i + 1 ) % 3; 00260 vertAdjMap[TRIANGLE - FIRST_FACE][i][2] = ( i + 2 ) % 3; 00261 } 00262 00263 for( i = 0; i < 4; ++i ) 00264 { 00265 vertAdjMap[QUADRILATERAL - FIRST_FACE][i][0] = 2; 00266 vertAdjMap[QUADRILATERAL - FIRST_FACE][i][1] = ( i + 1 ) % 4; 00267 vertAdjMap[QUADRILATERAL - FIRST_FACE][i][2] = ( i + 3 ) % 4; 00268 } 00269 00270 unsigned tet_corner_data[] = { 1, 2, 3, 0, 3, 2, 3, 0, 1, 2, 1, 0 }; 00271 for( i = 0; i < 4; ++i ) 00272 { 00273 vertAdjMap[TETRAHEDRON - FIRST_FACE][i][0] = 3; 00274 for( unsigned j = 0; j < 3; ++j ) 00275 vertAdjMap[TETRAHEDRON - FIRST_FACE][i][j + 1] = tet_corner_data[3 * i + j]; 00276 } 00277 00278 for( i = 0; i < 4; ++i ) 00279 { 00280 vertAdjMap[PYRAMID - FIRST_FACE][i][0] = 3; 00281 vertAdjMap[PYRAMID - FIRST_FACE][i][1] = ( i + 1 ) % 4; 00282 vertAdjMap[PYRAMID - FIRST_FACE][i][2] = ( i + 3 ) % 4; 00283 vertAdjMap[PYRAMID - FIRST_FACE][i][3] = 4; 00284 } 00285 vertAdjMap[PYRAMID - FIRST_FACE][4][0] = 4; 00286 for( i = 0; i < 4; i++ ) 00287 vertAdjMap[PYRAMID - FIRST_FACE][4][i + 1] = 3 - i; 00288 00289 for( i = 0; i < 4; ++i ) 00290 { 00291 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][0] = 3; 00292 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][1] = ( i + 1 ) % 4; 00293 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][2] = ( i + 3 ) % 4; 00294 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][3] = i + 4; 00295 } 00296 for( i = 4; i < 8; ++i ) 00297 { 00298 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][0] = 3; 00299 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][1] = ( i + 3 ) % 4 + 4; 00300 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][2] = ( i + 1 ) % 4 + 4; 00301 vertAdjMap[HEXAHEDRON - FIRST_FACE][i][3] = i - 4; 00302 } 00303 00304 for( i = 0; i < 3; ++i ) 00305 { 00306 vertAdjMap[PRISM - FIRST_FACE][i][0] = 3; 00307 vertAdjMap[PRISM - FIRST_FACE][i][1] = ( i + 1 ) % 3; 00308 vertAdjMap[PRISM - FIRST_FACE][i][2] = ( i + 2 ) % 3; 00309 vertAdjMap[PRISM - FIRST_FACE][i][3] = i + 3; 00310 } 00311 for( i = 3; i < 6; ++i ) 00312 { 00313 vertAdjMap[PRISM - FIRST_FACE][i][0] = 3; 00314 vertAdjMap[PRISM - FIRST_FACE][i][1] = ( i + 2 ) % 3 + 3; 00315 vertAdjMap[PRISM - FIRST_FACE][i][2] = ( i + 1 ) % 3 + 3; 00316 vertAdjMap[PRISM - FIRST_FACE][i][3] = i - 3; 00317 } 00318 00319 // Build reverse vertex-vertex adjacency index map 00320 const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON }; 00321 const int num_types = sizeof( types ) / sizeof( types[0] ); 00322 for( i = 0; i < num_types; ++i ) 00323 { 00324 const unsigned num_vert = corners( types[i] ); 00325 for( unsigned v = 0; v < num_vert; ++v ) 00326 { 00327 unsigned num_v_adj; 00328 const unsigned* v_adj = adjacent_vertices( types[i], v, num_v_adj ); 00329 unsigned* reverse = revVertAdjIdx[types[i] - FIRST_FACE][v]; 00330 reverse[0] = num_v_adj; 00331 00332 for( unsigned j = 0; j < num_v_adj; ++j ) 00333 { 00334 unsigned num_j_adj, k; 00335 const unsigned* j_adj = adjacent_vertices( types[i], v_adj[j], num_j_adj ); 00336 for( k = 0; k < num_j_adj && j_adj[k] != v; ++k ) 00337 ; 00338 assert( k < num_j_adj ); // If this fails, vertAdjMap is corrupt! 00339 reverse[j + 1] = k; 00340 } 00341 } 00342 } 00343 } 00344 00345 void TopologyInfo::higher_order( EntityTopology topo, 00346 unsigned num_nodes, 00347 bool& midedge, 00348 bool& midface, 00349 bool& midvol, 00350 MsqError& err ) 00351 { 00352 int ho = higher_order( topo, num_nodes, err ); 00353 midedge = (bool)( ( ho & ( 1 << 1 ) ) >> 1 ); 00354 midface = (bool)( ( ho & ( 1 << 2 ) ) >> 2 ); 00355 midvol = (bool)( ( ho & ( 1 << 3 ) ) >> 3 ); 00356 } 00357 00358 int TopologyInfo::higher_order( EntityTopology topo, unsigned num_nodes, MsqError& err ) 00359 { 00360 int result = 0; 00361 if( topo == POLYGON ) // polygons currently do not have higher order elements 00362 return 0; 00363 00364 if( topo >= MIXED || num_nodes < instance.adjMap[topo][0] ) 00365 { 00366 MSQ_SETERR( err )( "Invalid element topology", MsqError::INVALID_ARG ); 00367 return 0; 00368 } 00369 00370 unsigned dim = instance.dimMap[topo]; 00371 assert( num_nodes >= instance.adjMap[topo][0] ); 00372 unsigned nodes = num_nodes - instance.adjMap[topo][0]; 00373 unsigned edges = instance.adjMap[topo][1]; 00374 unsigned faces = instance.adjMap[topo][2]; 00375 if( edges && nodes >= edges ) 00376 { 00377 nodes -= edges; 00378 result |= 1 << 1; 00379 } 00380 if( faces && nodes >= faces ) 00381 { 00382 nodes -= faces; 00383 result |= 1 << 2; 00384 } 00385 if( 1 == nodes ) 00386 { 00387 if( 2 == dim ) 00388 { 00389 nodes -= 1; 00390 result |= 1 << 2; 00391 } 00392 else if( 3 == dim ) 00393 { 00394 nodes -= 1; 00395 result |= 1 << 3; 00396 } 00397 } 00398 00399 if( nodes ) 00400 { 00401 MSQ_SETERR( err )( "Invalid element topology", MsqError::INVALID_STATE ); 00402 } 00403 00404 return result; 00405 } 00406 00407 int TopologyInfo::higher_order_from_side( EntityTopology topo, 00408 unsigned num_nodes, 00409 unsigned side_dimension, 00410 unsigned side_number, 00411 MsqError& err ) 00412 { 00413 bool mids[4] = { true }; 00414 higher_order( topo, num_nodes, mids[1], mids[2], mids[3], err ); 00415 MSQ_ERRZERO( err ); 00416 00417 if( side_dimension > dimension( topo ) || side_number > adjacent( topo, side_dimension ) ) 00418 { 00419 MSQ_SETERR( err )( MsqError::INVALID_ARG, "Invalid side number: %u\n", side_number ); 00420 return 0; 00421 } 00422 00423 if( !mids[side_dimension] ) return -1; 00424 00425 int result = side_number; 00426 switch( side_dimension ) 00427 { 00428 case 3: 00429 if( mids[2] ) result += faces( topo ); 00430 case 2: 00431 if( mids[1] ) result += edges( topo ); 00432 case 1: 00433 result += corners( topo ); 00434 case 0: 00435 break; 00436 default: 00437 MSQ_SETERR( err )( MsqError::INVALID_ARG, "Invalid dimension: %u\n", side_dimension ); 00438 return 0; 00439 } 00440 return result; 00441 } 00442 00443 void TopologyInfo::side_from_higher_order( EntityTopology topo, 00444 unsigned num_nodes, 00445 unsigned node_number, 00446 unsigned& side_dim_out, 00447 unsigned& side_num_out, 00448 MsqError& err ) 00449 { 00450 bool midedge, midface, midvol; 00451 higher_order( topo, num_nodes, midedge, midface, midvol, err );MSQ_ERRRTN( err ); 00452 side_num_out = node_number; 00453 00454 if( side_num_out < corners( topo ) ) 00455 { 00456 side_dim_out = 0; 00457 return; 00458 } 00459 side_num_out -= corners( topo ); 00460 00461 if( midedge ) 00462 { 00463 if( side_num_out < edges( topo ) ) 00464 { 00465 side_dim_out = 1; 00466 return; 00467 } 00468 side_num_out -= edges( topo ); 00469 } 00470 00471 if( midface ) 00472 { 00473 if( side_num_out < faces( topo ) ) 00474 { 00475 side_dim_out = 2; 00476 return; 00477 } 00478 side_num_out -= faces( topo ); 00479 } 00480 00481 if( midvol && side_num_out == 0 ) 00482 { 00483 side_dim_out = 3; 00484 return; 00485 } 00486 00487 MSQ_SETERR( err )( MsqError::INVALID_ARG, "Invalid node index\n" ); 00488 } 00489 00490 const unsigned* TopologyInfo::edge_vertices( EntityTopology topo, unsigned edge, MsqError& err ) 00491 { 00492 if( topo < (EntityTopology)FIRST_FACE || topo > (EntityTopology)LAST_VOL || edge >= edges( topo ) ) 00493 { 00494 MSQ_SETERR( err )( MsqError::INVALID_ARG ); 00495 topo = (EntityTopology)FIRST_FACE; 00496 edge = 0; 00497 } 00498 00499 return instance.edgeMap[topo - FIRST_FACE][edge]; 00500 } 00501 00502 const unsigned* TopologyInfo::edge_vertices( EntityTopology topo, unsigned edge ) 00503 { 00504 if( topo < (EntityTopology)FIRST_FACE || topo > (EntityTopology)LAST_VOL || edge >= edges( topo ) ) 00505 { 00506 return 0; 00507 } 00508 return instance.edgeMap[topo - FIRST_FACE][edge]; 00509 } 00510 00511 const unsigned* TopologyInfo::face_vertices( EntityTopology topo, unsigned face, unsigned& length, MsqError& err ) 00512 { 00513 if( topo < (EntityTopology)FIRST_VOL || topo > (EntityTopology)LAST_VOL || face >= faces( topo ) ) 00514 { 00515 MSQ_SETERR( err )( MsqError::INVALID_ARG ); 00516 topo = (EntityTopology)FIRST_VOL; 00517 face = 0; 00518 } 00519 00520 length = instance.faceMap[topo - FIRST_VOL][face][0]; 00521 return instance.faceMap[topo - FIRST_VOL][face] + 1; 00522 } 00523 const unsigned* TopologyInfo::face_vertices( EntityTopology topo, unsigned face, unsigned& length ) 00524 { 00525 if( topo < (EntityTopology)FIRST_VOL || topo > (EntityTopology)LAST_VOL || face >= faces( topo ) ) 00526 { 00527 return 0; 00528 } 00529 00530 length = instance.faceMap[topo - FIRST_VOL][face][0]; 00531 return instance.faceMap[topo - FIRST_VOL][face] + 1; 00532 } 00533 00534 const unsigned* TopologyInfo::side_vertices( EntityTopology topo, 00535 unsigned dim, 00536 unsigned side, 00537 unsigned& count_out, 00538 MsqError& err ) 00539 { 00540 static const unsigned all[] = { 0, 1, 2, 3, 4, 5, 6, 7 }; 00541 const unsigned* result; 00542 00543 if( dim != 0 && dim == dimension( topo ) ) 00544 { 00545 count_out = corners( topo ); 00546 result = all; 00547 } 00548 else if( dim == 1 ) 00549 { 00550 count_out = 2; 00551 result = edge_vertices( topo, side, err ); 00552 } 00553 else if( dim == 2 ) 00554 { 00555 result = face_vertices( topo, side, count_out, err ); 00556 } 00557 else 00558 { 00559 MSQ_SETERR( err )( MsqError::INVALID_ARG ); 00560 count_out = 0; 00561 result = 0; 00562 } 00563 return result; 00564 } 00565 const unsigned* TopologyInfo::side_vertices( EntityTopology topo, unsigned dim, unsigned side, unsigned& count_out ) 00566 { 00567 static const unsigned all[] = { 0, 1, 2, 3, 4, 5, 6, 7 }; 00568 const unsigned* result; 00569 00570 if( dim != 0 && dim == dimension( topo ) ) 00571 { 00572 count_out = corners( topo ); 00573 result = all; 00574 } 00575 else if( dim == 1 ) 00576 { 00577 count_out = 2; 00578 result = edge_vertices( topo, side ); 00579 } 00580 else if( dim == 2 ) 00581 { 00582 result = face_vertices( topo, side, count_out ); 00583 } 00584 else 00585 { 00586 result = 0; 00587 } 00588 return result; 00589 } 00590 00591 void TopologyInfo::side_number( EntityTopology topo, 00592 unsigned num_nodes, 00593 unsigned node_index, 00594 unsigned& side_dim_out, 00595 unsigned& side_num_out, 00596 MsqError& err ) 00597 { 00598 if( topo >= (EntityTopology)MIXED || num_nodes < instance.adjMap[topo][0] ) 00599 { 00600 MSQ_SETERR( err )( "Invalid element topology", MsqError::INVALID_ARG ); 00601 return; 00602 } 00603 00604 unsigned nodes = instance.adjMap[topo][0]; 00605 unsigned edges = instance.adjMap[topo][1]; 00606 unsigned faces = instance.adjMap[topo][2]; 00607 side_num_out = node_index; 00608 00609 if( side_num_out < nodes ) 00610 { 00611 side_dim_out = 0; 00612 return; 00613 } 00614 num_nodes -= nodes; 00615 side_num_out -= nodes; 00616 00617 if( edges && num_nodes >= edges ) 00618 { 00619 if( side_num_out < edges ) 00620 { 00621 side_dim_out = 1; 00622 return; 00623 } 00624 num_nodes -= edges; 00625 side_num_out -= edges; 00626 } 00627 if( faces && num_nodes >= faces ) 00628 { 00629 if( side_num_out < faces ) 00630 { 00631 side_dim_out = 2; 00632 return; 00633 } 00634 num_nodes -= faces; 00635 side_num_out -= faces; 00636 } 00637 if( side_num_out == 0 ) 00638 { 00639 side_dim_out = instance.dimMap[topo]; 00640 side_num_out = 0; 00641 return; 00642 } 00643 00644 MSQ_SETERR( err )( MsqError::INVALID_ARG ); 00645 } 00646 00647 const unsigned* TopologyInfo::adjacent_vertices( EntityTopology topo, unsigned index, unsigned& num_adj_out ) 00648 { 00649 const unsigned count = corners( topo ); 00650 if( !count || index >= count ) 00651 { 00652 num_adj_out = 0; 00653 return 0; 00654 } 00655 00656 const unsigned* vect = instance.vertAdjMap[topo - FIRST_FACE][index]; 00657 num_adj_out = vect[0]; 00658 return vect + 1; 00659 } 00660 00661 const unsigned* TopologyInfo::reverse_vertex_adjacency_offsets( EntityTopology topo, 00662 unsigned index, 00663 unsigned& num_adj_out ) 00664 { 00665 const unsigned count = corners( topo ); 00666 if( !count || index >= count ) 00667 { 00668 num_adj_out = 0; 00669 return 0; 00670 } 00671 00672 const unsigned* vect = instance.revVertAdjIdx[topo - FIRST_FACE][index]; 00673 num_adj_out = vect[0]; 00674 return vect + 1; 00675 } 00676 00677 bool TopologyInfo::compare_sides( const size_t* verts1, 00678 EntityTopology type1, 00679 unsigned side1, 00680 const size_t* verts2, 00681 EntityTopology type2, 00682 unsigned side2, 00683 unsigned side_dim, 00684 MsqError& err ) 00685 { 00686 const unsigned *conn1, *conn2; 00687 unsigned len1, len2; 00688 00689 conn1 = side_vertices( type1, side_dim, side1, len1, err ); 00690 MSQ_ERRZERO( err ); 00691 conn2 = side_vertices( type2, side_dim, side2, len2, err ); 00692 MSQ_ERRZERO( err ); 00693 00694 // obviously not the same if different number of vertices 00695 // (triangular face cannot match quadrilateral face) 00696 if( len1 != len2 ) return false; 00697 00698 // Find location (i) in vertices of side of second element 00699 // that matches the first vertex in the side of the first 00700 // element. 00701 unsigned i, j; 00702 for( i = 0; i < len2; ++i ) 00703 if( verts1[conn1[0]] == verts2[conn2[i]] ) break; 00704 // If not found, then no match 00705 if( i == len2 ) return false; 00706 00707 // Try comparing side connectivity in forward order 00708 for( j = 1; j < len1; ++j ) 00709 if( verts1[conn1[j]] != verts2[conn2[( i + j ) % len2]] ) break; 00710 // If they match, we're done 00711 if( j == len1 ) return true; 00712 00713 // Try comparing in reverse order 00714 for( j = 1; j < len1; ++j ) 00715 if( verts1[conn1[j]] != verts2[conn2[( i + len2 - j ) % len2]] ) return false; 00716 // If here, matched in reverse order 00717 return true; 00718 } 00719 00720 unsigned TopologyInfo::find_edge( EntityTopology topo, 00721 const unsigned* side_vertices, 00722 bool& reversed_out, 00723 MsqError& err ) 00724 { 00725 if( dimension( topo ) <= 1 ) 00726 { 00727 MSQ_SETERR( err )( MsqError::INVALID_ARG, "Invalid element dimension" ); 00728 return (unsigned)-1; 00729 } 00730 00731 for( unsigned i = 0; i < edges( topo ); ++i ) 00732 { 00733 const unsigned* edge = edge_vertices( topo, i, err ); 00734 MSQ_ERRZERO( err ); 00735 00736 if( edge[0] == side_vertices[0] && edge[1] == side_vertices[1] ) 00737 { 00738 reversed_out = false; 00739 return i; 00740 } 00741 00742 if( edge[0] == side_vertices[1] && edge[1] == side_vertices[0] ) 00743 { 00744 reversed_out = true; 00745 return i; 00746 } 00747 } 00748 00749 MSQ_SETERR( err )( MsqError::INVALID_ARG, "No such edge" ); 00750 return (unsigned)-1; 00751 } 00752 00753 unsigned TopologyInfo::find_face( EntityTopology topo, 00754 const unsigned* side_vertices, 00755 unsigned num_vertices, 00756 bool& reversed_out, 00757 MsqError& err ) 00758 { 00759 if( dimension( topo ) <= 2 ) 00760 { 00761 MSQ_SETERR( err )( MsqError::INVALID_ARG, "Invalid element dimension" ); 00762 return (unsigned)-1; 00763 } 00764 00765 for( unsigned i = 0; i < faces( topo ); ++i ) 00766 { 00767 unsigned j, n, offset; 00768 const unsigned* face = face_vertices( topo, i, n, err ); 00769 MSQ_ERRZERO( err ); 00770 if( n != num_vertices ) continue; 00771 00772 for( offset = 0; offset < num_vertices; ++offset ) 00773 if( face[offset] == side_vertices[0] ) break; 00774 if( offset == num_vertices ) continue; 00775 00776 for( j = 1; j < num_vertices; ++j ) 00777 if( side_vertices[j] != face[( offset + j ) % num_vertices] ) break; 00778 if( j == num_vertices ) 00779 { 00780 reversed_out = false; 00781 return i; 00782 } 00783 00784 for( j = 1; j < num_vertices; ++j ) 00785 if( side_vertices[j] != face[( offset + num_vertices - j ) % num_vertices] ) break; 00786 if( j == num_vertices ) 00787 { 00788 reversed_out = true; 00789 return i; 00790 } 00791 } 00792 00793 MSQ_SETERR( err )( MsqError::INVALID_ARG, "No such face" ); 00794 return (unsigned)-1; 00795 } 00796 00797 void TopologyInfo::find_side( EntityTopology topo, 00798 const unsigned* side_vertices, 00799 unsigned num_vertices, 00800 unsigned& dimension_out, 00801 unsigned& number_out, 00802 bool& reversed_out, 00803 MsqError& err ) 00804 { 00805 switch( num_vertices ) 00806 { 00807 case 1: 00808 dimension_out = 0; 00809 number_out = *side_vertices; 00810 reversed_out = false; 00811 if( *side_vertices >= corners( topo ) ) MSQ_SETERR( err ) 00812 ( MsqError::INVALID_ARG, "Invalid corner number: %u\n", *side_vertices ); 00813 break; 00814 case 2: 00815 dimension_out = 1; 00816 number_out = find_edge( topo, side_vertices, reversed_out, err );MSQ_CHKERR( err ); 00817 break; 00818 case 3: 00819 case 4: 00820 dimension_out = 2; 00821 number_out = find_face( topo, side_vertices, num_vertices, reversed_out, err );MSQ_CHKERR( err ); 00822 break; 00823 default: 00824 MSQ_SETERR( err ) 00825 ( MsqError::UNSUPPORTED_ELEMENT, "Invalid number of side vertices: %u\n", num_vertices ); 00826 break; 00827 } 00828 } 00829 00830 } // namespace MBMesquite