MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "SphereDecomp.hpp" 00002 #include "moab/MeshTopoUtil.hpp" 00003 #include "moab/Range.hpp" 00004 #include "moab/CN.hpp" 00005 #include <math.h> 00006 #include <assert.h> 00007 #include <iostream> 00008 00009 #define RR \ 00010 if( MB_SUCCESS != result ) return result 00011 00012 const char* SUBDIV_VERTICES_TAG_NAME = "subdiv_vertices"; 00013 00014 using namespace moab; 00015 00016 SphereDecomp::SphereDecomp( Interface* impl ) 00017 { 00018 mbImpl = impl; 00019 } 00020 00021 ErrorCode SphereDecomp::build_sphere_mesh( const char* sphere_radii_tag_name, EntityHandle* hex_set ) 00022 { 00023 ErrorCode result = mbImpl->tag_get_handle( sphere_radii_tag_name, 1, MB_TYPE_DOUBLE, sphereRadiiTag );RR; 00024 00025 // need to make sure all interior edges and faces are created 00026 Range all_verts; 00027 result = mbImpl->get_entities_by_type( 0, MBVERTEX, all_verts );RR; 00028 MeshTopoUtil mtu( mbImpl ); 00029 result = mtu.construct_aentities( all_verts );RR; 00030 00031 // create tag to hold vertices 00032 result = mbImpl->tag_get_handle( SUBDIV_VERTICES_TAG_NAME, 9, MB_TYPE_HANDLE, subdivVerticesTag, 00033 MB_TAG_DENSE | MB_TAG_EXCL );RR; 00034 00035 // compute nodal positions for each dimension element 00036 result = compute_nodes( 1 );RR; 00037 result = compute_nodes( 2 );RR; 00038 result = compute_nodes( 3 );RR; 00039 00040 // build hex elements 00041 std::vector< EntityHandle > sphere_hexes, interstic_hexes; 00042 result = build_hexes( sphere_hexes, interstic_hexes );RR; 00043 00044 result = mbImpl->tag_delete( subdivVerticesTag );RR; 00045 00046 if( NULL != hex_set ) 00047 { 00048 if( 0 == *hex_set ) 00049 { 00050 EntityHandle this_set; 00051 // make a new set 00052 result = mbImpl->create_meshset( MESHSET_SET, this_set );RR; 00053 *hex_set = this_set; 00054 } 00055 00056 // save all the hexes to this set 00057 result = mbImpl->add_entities( *hex_set, &sphere_hexes[0], sphere_hexes.size() );RR; 00058 result = mbImpl->add_entities( *hex_set, &interstic_hexes[0], interstic_hexes.size() );RR; 00059 } 00060 00061 return result; 00062 } 00063 00064 ErrorCode SphereDecomp::compute_nodes( const int dim ) 00065 { 00066 // get facets of that dimension 00067 Range these_ents; 00068 const EntityType the_types[4] = { MBVERTEX, MBEDGE, MBTRI, MBTET }; 00069 00070 ErrorCode result = mbImpl->get_entities_by_dimension( 0, dim, these_ents );RR; 00071 assert( mbImpl->type_from_handle( *these_ents.begin() ) == the_types[dim] && 00072 mbImpl->type_from_handle( *these_ents.rbegin() ) == the_types[dim] ); 00073 00074 EntityHandle subdiv_vertices[9]; 00075 MeshTopoUtil mtu( mbImpl ); 00076 double avg_pos[3], vert_pos[12], new_vert_pos[12], new_new_vert_pos[3]; 00077 double radii[4], unitv[3]; 00078 int num_verts = CN::VerticesPerEntity( the_types[dim] ); 00079 00080 for( Range::iterator rit = these_ents.begin(); rit != these_ents.end(); ++rit ) 00081 { 00082 00083 // get vertices 00084 const EntityHandle* connect; 00085 int num_connect; 00086 result = mbImpl->get_connectivity( *rit, connect, num_connect );RR; 00087 00088 // compute center 00089 result = mtu.get_average_position( connect, num_connect, avg_pos );RR; 00090 00091 // create center vertex 00092 result = mbImpl->create_vertex( avg_pos, subdiv_vertices[num_verts] );RR; 00093 00094 // get coords of other vertices 00095 result = mbImpl->get_coords( connect, num_connect, vert_pos );RR; 00096 00097 // get radii associated with each vertex 00098 result = mbImpl->tag_get_data( sphereRadiiTag, connect, num_connect, radii );RR; 00099 00100 // compute subdiv vertex position for each vertex 00101 for( int i = 0; i < num_verts; i++ ) 00102 { 00103 for( int j = 0; j < 3; j++ ) 00104 unitv[j] = avg_pos[j] - vert_pos[3 * i + j]; 00105 double vlength = sqrt( unitv[0] * unitv[0] + unitv[1] * unitv[1] + unitv[2] * unitv[2] ); 00106 if( vlength < radii[i] ) 00107 { 00108 std::cout << "Radius too large at vertex " << i << std::endl; 00109 result = MB_FAILURE; 00110 continue; 00111 } 00112 00113 for( int j = 0; j < 3; j++ ) 00114 unitv[j] /= vlength; 00115 00116 for( int j = 0; j < 3; j++ ) 00117 new_vert_pos[3 * i + j] = vert_pos[3 * i + j] + radii[i] * unitv[j]; 00118 00119 // create vertex at this position 00120 ErrorCode tmp_result = mbImpl->create_vertex( &new_vert_pos[3 * i], subdiv_vertices[i] ); 00121 if( MB_SUCCESS != tmp_result ) result = tmp_result; 00122 } 00123 00124 if( MB_SUCCESS != result ) return result; 00125 00126 // compute subdiv vertex positions for vertices inside spheres; just mid-pt between 00127 // previous subdiv vertex and corner vertex 00128 for( int i = 0; i < num_verts; i++ ) 00129 { 00130 for( int j = 0; j < 3; j++ ) 00131 new_new_vert_pos[j] = .5 * ( vert_pos[3 * i + j] + new_vert_pos[3 * i + j] ); 00132 00133 result = mbImpl->create_vertex( new_new_vert_pos, subdiv_vertices[num_verts + 1 + i] ); 00134 } 00135 00136 // set the tag 00137 result = mbImpl->tag_set_data( subdivVerticesTag, &( *rit ), 1, subdiv_vertices );RR; 00138 } 00139 00140 return result; 00141 } 00142 00143 ErrorCode SphereDecomp::build_hexes( std::vector< EntityHandle >& sphere_hexes, 00144 std::vector< EntityHandle >& interstic_hexes ) 00145 { 00146 // build hexes inside each tet element separately 00147 Range tets; 00148 ErrorCode result = mbImpl->get_entities_by_type( 0, MBTET, tets );RR; 00149 00150 for( Range::iterator vit = tets.begin(); vit != tets.end(); ++vit ) 00151 { 00152 result = subdivide_tet( *vit, sphere_hexes, interstic_hexes );RR; 00153 } 00154 00155 return MB_SUCCESS; 00156 } 00157 00158 ErrorCode SphereDecomp::subdivide_tet( EntityHandle tet, std::vector< EntityHandle >& sphere_hexes, 00159 std::vector< EntityHandle >& interstic_hexes ) 00160 { 00161 // 99: (#subdiv_verts/entity=9) * (#edges=6 + #faces=4 + 1=tet) 00162 EntityHandle subdiv_verts[99]; 00163 00164 // get tet connectivity 00165 std::vector< EntityHandle > tet_conn; 00166 ErrorCode result = mbImpl->get_connectivity( &tet, 1, tet_conn );RR; 00167 00168 for( int dim = 1; dim <= 3; dim++ ) 00169 { 00170 // get entities of this dimension 00171 std::vector< EntityHandle > ents; 00172 if( dim != 3 ) 00173 { 00174 result = mbImpl->get_adjacencies( &tet, 1, dim, false, ents );RR; 00175 } 00176 else 00177 ents.push_back( tet ); 00178 00179 // for each, get subdiv verts & put into vector 00180 for( std::vector< EntityHandle >::iterator vit = ents.begin(); vit != ents.end(); ++vit ) 00181 { 00182 result = retrieve_subdiv_verts( tet, *vit, &tet_conn[0], dim, subdiv_verts );RR; 00183 } 00184 } 00185 00186 // ok, subdiv_verts are in canonical order; now create the hexes, using pre-computed templates 00187 00188 // Templates are specified in terms of the vertices making up each hex; vertices are specified 00189 // by specifying the facet index and type they resolve, and the index of that vertex in that 00190 // facet's subdivision vertices list. 00191 00192 // Each facet is subdivided into: 00193 // - a mid vertex 00194 // - one vertex for each corner vertex on the facet (located on a line between the mid vertex 00195 // and 00196 // the corresponding corner vertex, a distance equal to the sphere radius away from the corner 00197 // vertex) 00198 // - one vertex midway between each corner vertex and the corresponding "sphere surface" vertex 00199 // For edges, tris and tets this gives 5, 7 and 9 subdivision vertices, respectively. 00200 // Subdivision vertices appear in the list in the order: sphere surface vertices, mid vertex, 00201 // sphere interior vertices. In each of those sub lists, vertices are listed in the canonical 00202 // order of the corresponding corner vertices for that facet. 00203 00204 // Subdivision vertices for facetes are indexed by listing the facet type they resolve (EDGE, 00205 // FACE, TET), the index of that facet (integer = 0..5, 0..3, 0 for edges, tris, tet, resp), and 00206 // subdivision index (AINDEX..EINDEX for edges, AINDEX..GINDEX for tris, AINDEX..IINDEX for 00207 // tets). 00208 00209 // Subdivision vertices for all facets of a tet are stored in one subdivision vertex vector, in 00210 // order of increasing facet dimension and index (index varies fastest). The ESV, FSV, and TSV 00211 // macros are used to compute the indices into that vector for various parameters. The CV macro 00212 // is used to index into the tet connectivity vector. 00213 00214 // Subdivision templates for splitting the tet into 28 hexes were derived by hand, and are 00215 // listed below (using the indexing scheme described above). 00216 00217 #define EDGE 0 00218 #define FACE 1 00219 #define TET 2 00220 #define AINDEX 0 00221 #define BINDEX 1 00222 #define CINDEX 2 00223 #define DINDEX 3 00224 #define EINDEX 4 00225 #define FINDEX 5 00226 #define GINDEX 6 00227 #define HINDEX 7 00228 #define IINDEX 8 00229 #define V0INDEX 0 00230 #define V1INDEX 1 00231 #define V2INDEX 2 00232 #define V3INDEX 3 00233 #define CV( a ) tet_conn[a] 00234 #define ESV( a, b ) subdiv_verts[a * 9 + b] 00235 #define FSV( a, b ) subdiv_verts[54 + a * 9 + b] 00236 #define TSV( a, b ) subdiv_verts[90 + a * 9 + b] 00237 00238 EntityHandle this_connect[8], this_hex; 00239 00240 // first, interstices hexes, three per vertex/spherical surface 00241 // V0: 00242 int i = 0; 00243 this_connect[i++] = ESV( 0, AINDEX ); 00244 this_connect[i++] = ESV( 0, CINDEX ); 00245 this_connect[i++] = FSV( 3, DINDEX ); 00246 this_connect[i++] = FSV( 3, AINDEX ); 00247 this_connect[i++] = FSV( 0, AINDEX ); 00248 this_connect[i++] = FSV( 0, DINDEX ); 00249 this_connect[i++] = TSV( 0, EINDEX ); 00250 this_connect[i++] = TSV( 0, AINDEX ); 00251 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00252 interstic_hexes.push_back( this_hex ); 00253 00254 i = 0; 00255 this_connect[i++] = FSV( 0, AINDEX ); 00256 this_connect[i++] = FSV( 0, DINDEX ); 00257 this_connect[i++] = TSV( 0, EINDEX ); 00258 this_connect[i++] = TSV( 0, AINDEX ); 00259 this_connect[i++] = ESV( 3, AINDEX ); 00260 this_connect[i++] = ESV( 3, CINDEX ); 00261 this_connect[i++] = FSV( 2, DINDEX ); 00262 this_connect[i++] = FSV( 2, AINDEX ); 00263 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00264 interstic_hexes.push_back( this_hex ); 00265 00266 i = 0; 00267 this_connect[i++] = FSV( 3, AINDEX ); 00268 this_connect[i++] = FSV( 3, DINDEX ); 00269 this_connect[i++] = ESV( 2, CINDEX ); 00270 this_connect[i++] = ESV( 2, BINDEX ); 00271 this_connect[i++] = TSV( 0, AINDEX ); 00272 this_connect[i++] = TSV( 0, EINDEX ); 00273 this_connect[i++] = FSV( 2, DINDEX ); 00274 this_connect[i++] = FSV( 2, AINDEX ); 00275 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00276 interstic_hexes.push_back( this_hex ); 00277 00278 // V1: 00279 i = 0; 00280 this_connect[i++] = ESV( 0, CINDEX ); 00281 this_connect[i++] = ESV( 0, BINDEX ); 00282 this_connect[i++] = FSV( 3, CINDEX ); 00283 this_connect[i++] = FSV( 3, DINDEX ); 00284 this_connect[i++] = FSV( 0, DINDEX ); 00285 this_connect[i++] = FSV( 0, BINDEX ); 00286 this_connect[i++] = TSV( 0, BINDEX ); 00287 this_connect[i++] = TSV( 0, EINDEX ); 00288 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00289 interstic_hexes.push_back( this_hex ); 00290 00291 i = 0; 00292 this_connect[i++] = FSV( 0, DINDEX ); 00293 this_connect[i++] = FSV( 0, BINDEX ); 00294 this_connect[i++] = TSV( 0, BINDEX ); 00295 this_connect[i++] = TSV( 0, EINDEX ); 00296 this_connect[i++] = ESV( 4, CINDEX ); 00297 this_connect[i++] = ESV( 4, AINDEX ); 00298 this_connect[i++] = FSV( 1, AINDEX ); 00299 this_connect[i++] = FSV( 1, DINDEX ); 00300 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00301 interstic_hexes.push_back( this_hex ); 00302 00303 i = 0; 00304 this_connect[i++] = FSV( 1, DINDEX ); 00305 this_connect[i++] = FSV( 1, AINDEX ); 00306 this_connect[i++] = TSV( 0, BINDEX ); 00307 this_connect[i++] = TSV( 0, EINDEX ); 00308 this_connect[i++] = ESV( 1, CINDEX ); 00309 this_connect[i++] = ESV( 1, AINDEX ); 00310 this_connect[i++] = FSV( 3, CINDEX ); 00311 this_connect[i++] = FSV( 3, DINDEX ); 00312 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00313 interstic_hexes.push_back( this_hex ); 00314 00315 // V2: 00316 i = 0; 00317 this_connect[i++] = FSV( 3, DINDEX ); 00318 this_connect[i++] = ESV( 1, CINDEX ); 00319 this_connect[i++] = ESV( 1, BINDEX ); 00320 this_connect[i++] = FSV( 3, BINDEX ); 00321 this_connect[i++] = TSV( 0, EINDEX ); 00322 this_connect[i++] = FSV( 1, DINDEX ); 00323 this_connect[i++] = FSV( 1, BINDEX ); 00324 this_connect[i++] = TSV( 0, CINDEX ); 00325 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00326 interstic_hexes.push_back( this_hex ); 00327 00328 i = 0; 00329 this_connect[i++] = TSV( 0, EINDEX ); 00330 this_connect[i++] = FSV( 1, DINDEX ); 00331 this_connect[i++] = FSV( 1, BINDEX ); 00332 this_connect[i++] = TSV( 0, CINDEX ); 00333 this_connect[i++] = FSV( 2, DINDEX ); 00334 this_connect[i++] = ESV( 5, CINDEX ); 00335 this_connect[i++] = ESV( 5, AINDEX ); 00336 this_connect[i++] = FSV( 2, CINDEX ); 00337 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00338 interstic_hexes.push_back( this_hex ); 00339 00340 i = 0; 00341 this_connect[i++] = TSV( 0, CINDEX ); 00342 this_connect[i++] = FSV( 2, CINDEX ); 00343 this_connect[i++] = ESV( 2, AINDEX ); 00344 this_connect[i++] = FSV( 3, BINDEX ); 00345 this_connect[i++] = TSV( 0, EINDEX ); 00346 this_connect[i++] = FSV( 2, DINDEX ); 00347 this_connect[i++] = ESV( 2, CINDEX ); 00348 this_connect[i++] = FSV( 3, DINDEX ); 00349 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00350 interstic_hexes.push_back( this_hex ); 00351 00352 // V3: 00353 i = 0; 00354 this_connect[i++] = TSV( 0, EINDEX ); 00355 this_connect[i++] = FSV( 1, DINDEX ); 00356 this_connect[i++] = ESV( 5, CINDEX ); 00357 this_connect[i++] = FSV( 2, DINDEX ); 00358 this_connect[i++] = TSV( 0, DINDEX ); 00359 this_connect[i++] = FSV( 1, CINDEX ); 00360 this_connect[i++] = ESV( 5, BINDEX ); 00361 this_connect[i++] = FSV( 2, BINDEX ); 00362 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00363 interstic_hexes.push_back( this_hex ); 00364 00365 i = 0; 00366 this_connect[i++] = FSV( 0, DINDEX ); 00367 this_connect[i++] = ESV( 4, CINDEX ); 00368 this_connect[i++] = FSV( 1, DINDEX ); 00369 this_connect[i++] = TSV( 0, EINDEX ); 00370 this_connect[i++] = FSV( 0, CINDEX ); 00371 this_connect[i++] = ESV( 4, BINDEX ); 00372 this_connect[i++] = FSV( 1, CINDEX ); 00373 this_connect[i++] = TSV( 0, DINDEX ); 00374 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00375 interstic_hexes.push_back( this_hex ); 00376 00377 i = 0; 00378 this_connect[i++] = ESV( 3, CINDEX ); 00379 this_connect[i++] = FSV( 0, DINDEX ); 00380 this_connect[i++] = TSV( 0, EINDEX ); 00381 this_connect[i++] = FSV( 2, DINDEX ); 00382 this_connect[i++] = ESV( 3, BINDEX ); 00383 this_connect[i++] = FSV( 0, CINDEX ); 00384 this_connect[i++] = TSV( 0, DINDEX ); 00385 this_connect[i++] = FSV( 2, BINDEX ); 00386 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00387 interstic_hexes.push_back( this_hex ); 00388 00389 // now, the sphere interiors, four hexes per vertex sphere 00390 00391 // V0: 00392 i = 0; 00393 this_connect[i++] = CV( V0INDEX ); 00394 this_connect[i++] = ESV( 0, DINDEX ); 00395 this_connect[i++] = FSV( 3, EINDEX ); 00396 this_connect[i++] = ESV( 2, EINDEX ); 00397 this_connect[i++] = ESV( 3, DINDEX ); 00398 this_connect[i++] = FSV( 0, EINDEX ); 00399 this_connect[i++] = TSV( 0, FINDEX ); 00400 this_connect[i++] = FSV( 2, EINDEX ); 00401 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00402 sphere_hexes.push_back( this_hex ); 00403 00404 i = 0; 00405 this_connect[i++] = ESV( 0, DINDEX ); 00406 this_connect[i++] = ESV( 0, AINDEX ); 00407 this_connect[i++] = FSV( 3, AINDEX ); 00408 this_connect[i++] = FSV( 3, EINDEX ); 00409 this_connect[i++] = FSV( 0, EINDEX ); 00410 this_connect[i++] = FSV( 0, AINDEX ); 00411 this_connect[i++] = TSV( 0, AINDEX ); 00412 this_connect[i++] = TSV( 0, FINDEX ); 00413 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00414 sphere_hexes.push_back( this_hex ); 00415 00416 i = 0; 00417 this_connect[i++] = FSV( 3, EINDEX ); 00418 this_connect[i++] = FSV( 3, AINDEX ); 00419 this_connect[i++] = ESV( 2, BINDEX ); 00420 this_connect[i++] = ESV( 2, EINDEX ); 00421 this_connect[i++] = TSV( 0, FINDEX ); 00422 this_connect[i++] = TSV( 0, AINDEX ); 00423 this_connect[i++] = FSV( 2, AINDEX ); 00424 this_connect[i++] = FSV( 2, EINDEX ); 00425 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00426 sphere_hexes.push_back( this_hex ); 00427 00428 i = 0; 00429 this_connect[i++] = TSV( 0, FINDEX ); 00430 this_connect[i++] = TSV( 0, AINDEX ); 00431 this_connect[i++] = FSV( 2, AINDEX ); 00432 this_connect[i++] = FSV( 2, EINDEX ); 00433 this_connect[i++] = FSV( 0, EINDEX ); 00434 this_connect[i++] = FSV( 0, AINDEX ); 00435 this_connect[i++] = ESV( 3, AINDEX ); 00436 this_connect[i++] = ESV( 3, DINDEX ); 00437 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00438 sphere_hexes.push_back( this_hex ); 00439 00440 // V1: 00441 i = 0; 00442 this_connect[i++] = CV( V1INDEX ); 00443 this_connect[i++] = ESV( 1, DINDEX ); 00444 this_connect[i++] = FSV( 3, GINDEX ); 00445 this_connect[i++] = ESV( 0, EINDEX ); 00446 this_connect[i++] = ESV( 4, DINDEX ); 00447 this_connect[i++] = FSV( 1, EINDEX ); 00448 this_connect[i++] = TSV( 0, GINDEX ); 00449 this_connect[i++] = FSV( 0, FINDEX ); 00450 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00451 sphere_hexes.push_back( this_hex ); 00452 00453 i = 0; 00454 this_connect[i++] = FSV( 3, GINDEX ); 00455 this_connect[i++] = ESV( 1, DINDEX ); 00456 this_connect[i++] = ESV( 1, AINDEX ); 00457 this_connect[i++] = FSV( 3, CINDEX ); 00458 this_connect[i++] = TSV( 0, GINDEX ); 00459 this_connect[i++] = FSV( 1, EINDEX ); 00460 this_connect[i++] = FSV( 1, AINDEX ); 00461 this_connect[i++] = TSV( 0, BINDEX ); 00462 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00463 sphere_hexes.push_back( this_hex ); 00464 00465 i = 0; 00466 this_connect[i++] = TSV( 0, GINDEX ); 00467 this_connect[i++] = FSV( 1, EINDEX ); 00468 this_connect[i++] = FSV( 1, AINDEX ); 00469 this_connect[i++] = TSV( 0, BINDEX ); 00470 this_connect[i++] = FSV( 0, FINDEX ); 00471 this_connect[i++] = ESV( 4, DINDEX ); 00472 this_connect[i++] = ESV( 4, AINDEX ); 00473 this_connect[i++] = FSV( 0, BINDEX ); 00474 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00475 sphere_hexes.push_back( this_hex ); 00476 00477 i = 0; 00478 this_connect[i++] = ESV( 0, BINDEX ); 00479 this_connect[i++] = ESV( 0, EINDEX ); 00480 this_connect[i++] = FSV( 3, GINDEX ); 00481 this_connect[i++] = FSV( 3, CINDEX ); 00482 this_connect[i++] = FSV( 0, BINDEX ); 00483 this_connect[i++] = FSV( 0, FINDEX ); 00484 this_connect[i++] = TSV( 0, GINDEX ); 00485 this_connect[i++] = TSV( 0, BINDEX ); 00486 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00487 sphere_hexes.push_back( this_hex ); 00488 00489 // V2: 00490 i = 0; 00491 this_connect[i++] = ESV( 1, BINDEX ); 00492 this_connect[i++] = ESV( 1, EINDEX ); 00493 this_connect[i++] = FSV( 3, FINDEX ); 00494 this_connect[i++] = FSV( 3, BINDEX ); 00495 this_connect[i++] = FSV( 1, BINDEX ); 00496 this_connect[i++] = FSV( 1, FINDEX ); 00497 this_connect[i++] = TSV( 0, HINDEX ); 00498 this_connect[i++] = TSV( 0, CINDEX ); 00499 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00500 sphere_hexes.push_back( this_hex ); 00501 00502 i = 0; 00503 this_connect[i++] = FSV( 3, FINDEX ); 00504 this_connect[i++] = ESV( 1, EINDEX ); 00505 this_connect[i++] = CV( V2INDEX ); 00506 this_connect[i++] = ESV( 2, DINDEX ); 00507 this_connect[i++] = TSV( 0, HINDEX ); 00508 this_connect[i++] = FSV( 1, FINDEX ); 00509 this_connect[i++] = ESV( 5, DINDEX ); 00510 this_connect[i++] = FSV( 2, GINDEX ); 00511 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00512 sphere_hexes.push_back( this_hex ); 00513 00514 i = 0; 00515 this_connect[i++] = TSV( 0, HINDEX ); 00516 this_connect[i++] = FSV( 1, FINDEX ); 00517 this_connect[i++] = ESV( 5, DINDEX ); 00518 this_connect[i++] = FSV( 2, GINDEX ); 00519 this_connect[i++] = TSV( 0, CINDEX ); 00520 this_connect[i++] = FSV( 1, BINDEX ); 00521 this_connect[i++] = ESV( 5, AINDEX ); 00522 this_connect[i++] = FSV( 2, CINDEX ); 00523 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00524 sphere_hexes.push_back( this_hex ); 00525 00526 i = 0; 00527 this_connect[i++] = FSV( 3, BINDEX ); 00528 this_connect[i++] = FSV( 3, FINDEX ); 00529 this_connect[i++] = ESV( 2, DINDEX ); 00530 this_connect[i++] = ESV( 2, AINDEX ); 00531 this_connect[i++] = TSV( 0, CINDEX ); 00532 this_connect[i++] = TSV( 0, HINDEX ); 00533 this_connect[i++] = FSV( 2, GINDEX ); 00534 this_connect[i++] = FSV( 2, CINDEX ); 00535 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00536 sphere_hexes.push_back( this_hex ); 00537 00538 // V3: 00539 i = 0; 00540 this_connect[i++] = FSV( 0, CINDEX ); 00541 this_connect[i++] = ESV( 4, BINDEX ); 00542 this_connect[i++] = FSV( 1, CINDEX ); 00543 this_connect[i++] = TSV( 0, DINDEX ); 00544 this_connect[i++] = FSV( 0, GINDEX ); 00545 this_connect[i++] = ESV( 4, EINDEX ); 00546 this_connect[i++] = FSV( 1, GINDEX ); 00547 this_connect[i++] = TSV( 0, IINDEX ); 00548 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00549 sphere_hexes.push_back( this_hex ); 00550 00551 i = 0; 00552 this_connect[i++] = ESV( 3, BINDEX ); 00553 this_connect[i++] = FSV( 0, CINDEX ); 00554 this_connect[i++] = TSV( 0, DINDEX ); 00555 this_connect[i++] = FSV( 2, BINDEX ); 00556 this_connect[i++] = ESV( 3, EINDEX ); 00557 this_connect[i++] = FSV( 0, GINDEX ); 00558 this_connect[i++] = TSV( 0, IINDEX ); 00559 this_connect[i++] = FSV( 2, FINDEX ); 00560 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00561 sphere_hexes.push_back( this_hex ); 00562 00563 i = 0; 00564 this_connect[i++] = TSV( 0, DINDEX ); 00565 this_connect[i++] = FSV( 1, CINDEX ); 00566 this_connect[i++] = ESV( 5, BINDEX ); 00567 this_connect[i++] = FSV( 2, BINDEX ); 00568 this_connect[i++] = TSV( 0, IINDEX ); 00569 this_connect[i++] = FSV( 1, GINDEX ); 00570 this_connect[i++] = ESV( 5, EINDEX ); 00571 this_connect[i++] = FSV( 2, FINDEX ); 00572 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00573 sphere_hexes.push_back( this_hex ); 00574 00575 i = 0; 00576 this_connect[i++] = FSV( 0, GINDEX ); 00577 this_connect[i++] = ESV( 4, EINDEX ); 00578 this_connect[i++] = FSV( 1, GINDEX ); 00579 this_connect[i++] = TSV( 0, IINDEX ); 00580 this_connect[i++] = ESV( 3, EINDEX ); 00581 this_connect[i++] = CV( V3INDEX ); 00582 this_connect[i++] = ESV( 5, EINDEX ); 00583 this_connect[i++] = FSV( 2, FINDEX ); 00584 result = mbImpl->create_element( MBHEX, this_connect, 8, this_hex );RR; 00585 sphere_hexes.push_back( this_hex ); 00586 00587 return result; 00588 } 00589 00590 ErrorCode SphereDecomp::retrieve_subdiv_verts( EntityHandle tet, EntityHandle this_ent, const EntityHandle* tet_conn, 00591 const int dim, EntityHandle* subdiv_verts ) 00592 { 00593 // get the subdiv verts for this entity 00594 ErrorCode result; 00595 00596 // if it's a tet, just put them on the end & return 00597 if( tet == this_ent ) 00598 { 00599 result = mbImpl->tag_get_data( subdivVerticesTag, &this_ent, 1, &subdiv_verts[90] ); 00600 return MB_SUCCESS; 00601 } 00602 00603 // if it's a sub-entity, need to find index, relative orientation, and offset 00604 // get connectivity of sub-entity 00605 std::vector< EntityHandle > this_conn; 00606 result = mbImpl->get_connectivity( &this_ent, 1, this_conn );RR; 00607 00608 // get relative orientation 00609 std::vector< int > conn_tet_indices( this_conn.size() ); 00610 for( size_t i = 0; i < this_conn.size(); ++i ) 00611 conn_tet_indices[i] = std::find( tet_conn, tet_conn + 4, this_conn[i] ) - tet_conn; 00612 int sense, side_no, offset; 00613 int success = CN::SideNumber( MBTET, &conn_tet_indices[0], this_conn.size(), dim, side_no, sense, offset ); 00614 if( -1 == success ) return MB_FAILURE; 00615 00616 // start of this entity's subdiv_verts; edges go first, then preceding sides, then this one; 00617 // this assumes 6 edges/tet 00618 EntityHandle* subdiv_start = &subdiv_verts[( ( dim - 1 ) * 6 + side_no ) * 9]; 00619 00620 // get subdiv_verts and put them into proper place 00621 result = mbImpl->tag_get_data( subdivVerticesTag, &this_ent, 1, subdiv_start ); 00622 00623 // could probably do this more elegantly, but isn't worth it 00624 #define SWITCH( a, b ) \ 00625 { \ 00626 EntityHandle tmp_handle = a; \ 00627 a = b; \ 00628 b = tmp_handle; \ 00629 } 00630 switch( dim ) 00631 { 00632 case 1: 00633 if( offset != 0 || sense == -1 ) 00634 { 00635 SWITCH( subdiv_start[0], subdiv_start[1] ); 00636 SWITCH( subdiv_start[3], subdiv_start[4] ); 00637 } 00638 break; 00639 case 2: 00640 // rotate first 00641 if( 0 != offset ) 00642 { 00643 std::rotate( subdiv_start, subdiv_start + offset, subdiv_start + 3 ); 00644 std::rotate( subdiv_start + 4, subdiv_start + 4 + offset, subdiv_start + 7 ); 00645 } 00646 // now flip, if necessary 00647 if( -1 == sense ) 00648 { 00649 SWITCH( subdiv_start[1], subdiv_start[2] ); 00650 SWITCH( subdiv_start[5], subdiv_start[6] ); 00651 } 00652 break; 00653 default: 00654 return MB_FAILURE; 00655 } 00656 00657 // ok, we're done 00658 return MB_SUCCESS; 00659 }