Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
skin.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <set>
00003 #include <limits>
00004 #include <ctime>
00005 #include <vector>
00006 #include <cstdlib>
00007 #include <cstdio>
00008 #include <cassert>
00009 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00010 #include <unistd.h>
00011 #include <sys/types.h>
00012 #include <sys/stat.h>
00013 #endif
00014 #include <fcntl.h>
00015 #include "moab/Interface.hpp"
00016 #include "MBTagConventions.hpp"
00017 #include "moab/Core.hpp"
00018 #include "moab/Range.hpp"
00019 #include "moab/Skinner.hpp"
00020 #include "moab/AdaptiveKDTree.hpp"
00021 #include "moab/CN.hpp"
00022 
00023 using namespace moab;
00024 
00025 static void get_time_mem( double& tot_time, double& tot_mem );
00026 
00027 // Different platforms follow different conventions for usage
00028 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00029 #include <sys/resource.h>
00030 #endif
00031 #ifdef SOLARIS
00032 extern "C" int getrusage( int, struct rusage* );
00033 #ifndef RUSAGE_SELF
00034 #include </usr/ucbinclude/sys/rusage.h>
00035 #endif
00036 #endif
00037 
00038 const char DEFAULT_FIXED_TAG[] = "fixed";
00039 const int MIN_EDGE_LEN_DENOM   = 4;
00040 
00041 #define CHKERROR( A )                                                        \
00042     do                                                                       \
00043     {                                                                        \
00044         if( MB_SUCCESS != ( A ) )                                            \
00045         {                                                                    \
00046             std::cerr << "Internal error at line " << __LINE__ << std::endl; \
00047             return 3;                                                        \
00048         }                                                                    \
00049     } while( false )
00050 
00051 static ErrorCode merge_duplicate_vertices( Interface&, double epsilon );
00052 static ErrorCode min_edge_length( Interface&, double& result );
00053 
00054 static void usage( const char* argv0, bool help = false )
00055 {
00056     std::ostream& str = help ? std::cout : std::cerr;
00057 
00058     str << "Usage: " << argv0
00059         << " [-b <block_num> [-b ...] ] [-l] [-m] [-M <n>] [-p] [-s <sideset_num>] [-S] [-t|-T "
00060            "<name>] [-w] [-v|-V <n>]"
00061         << " <input_file> [<output_file>]" << std::endl;
00062     str << "Help : " << argv0 << " -h" << std::endl;
00063     if( !help ) exit( 1 );
00064 
00065     str << "Options: " << std::endl;
00066     str << "-a : Compute skin using vert-elem adjacencies (more memory, less time)." << std::endl;
00067     str << "-b <block_num> : Compute skin only for material set/block <block_num>." << std::endl;
00068     str << "-p : Print cpu & memory performance." << std::endl;
00069     str << "-s <sideset_num> : Put skin in neumann set/sideset <sideset_num>." << std::endl;
00070     str << "-S : Look for and use structured mesh information to speed up skinning." << std::endl;
00071     str << "-t : Set '" << DEFAULT_FIXED_TAG << "' tag on skin vertices." << std::endl;
00072     str << "-T <name> : Create tag with specified name and set to 1 on skin vertices." << std::endl;
00073     str << "-w : Write out whole mesh (otherwise just writes skin)." << std::endl;
00074     str << "-m : consolidate duplicate vertices" << std::endl;
00075     str << "-M <n> : consolidate duplicate vertices with specified tolerance. "
00076            "(Default: min_edge_length/"
00077         << MIN_EDGE_LEN_DENOM << ")" << std::endl;
00078     str << "-l : List total numbers of entities and vertices in skin." << std::endl;
00079     exit( 0 );
00080 }
00081 
00082 int main( int argc, char* argv[] )
00083 {
00084     int i = 1;
00085     std::vector< int > matsets;
00086     int neuset_num = -1;
00087     bool write_tag = false, write_whole_mesh = false;
00088     bool print_perf         = false;
00089     bool use_vert_elem_adjs = false;
00090     bool merge_vertices     = false;
00091     double merge_epsilon    = -1;
00092     bool list_skin          = false;
00093     bool use_scd            = false;
00094     const char* fixed_tag   = DEFAULT_FIXED_TAG;
00095     const char *input_file = 0, *output_file = 0;
00096 
00097     bool no_more_flags = false;
00098     char* endptr       = 0;
00099     long block         = 0;  // initialize to eliminate compiler warning
00100     while( i < argc )
00101     {
00102         if( !no_more_flags && argv[i][0] == '-' )
00103         {
00104             const int f = i++;
00105             for( int j = 1; argv[f][j]; ++j )
00106             {
00107                 switch( argv[f][j] )
00108                 {
00109                     case 'a':
00110                         use_vert_elem_adjs = true;
00111                         break;
00112                     case 'p':
00113                         print_perf = true;
00114                         break;
00115                     case 't':
00116                         write_tag = true;
00117                         break;
00118                     case 'w':
00119                         write_whole_mesh = true;
00120                         break;
00121                     case 'm':
00122                         merge_vertices = true;
00123                         break;
00124                     case '-':
00125                         no_more_flags = true;
00126                         break;
00127                     case 'h':
00128                         usage( argv[0], true );
00129                         break;
00130                     case 'l':
00131                         list_skin = true;
00132                         break;
00133                     case 'S':
00134                         use_scd = true;
00135                         break;
00136                     case 'b':
00137                         if( i == argc || 0 >= ( block = strtol( argv[i], &endptr, 0 ) ) || *endptr )
00138                         {
00139                             std::cerr << "Expected positive integer following '-b' flag" << std::endl;
00140                             usage( argv[0] );
00141                         }
00142                         matsets.push_back( (int)block );
00143                         ++i;
00144                         break;
00145                     case 's':
00146                         if( i == argc || 0 >= ( neuset_num = strtol( argv[i], &endptr, 0 ) ) || *endptr )
00147                         {
00148                             std::cerr << "Expected positive integer following '-s' flag" << std::endl;
00149                             usage( argv[0] );
00150                         }
00151                         ++i;
00152                         break;
00153                     case 'T':
00154                         if( i == argc || argv[i][0] == '-' )
00155                         {
00156                             std::cerr << "Expected tag name following '-T' flag" << std::endl;
00157                             usage( argv[0] );
00158                         }
00159                         fixed_tag = argv[i++];
00160                         break;
00161                     case 'M':
00162                         if( i == argc || 0.0 > ( merge_epsilon = strtod( argv[i], &endptr ) ) || *endptr )
00163                         {
00164                             std::cerr << "Expected positive numeric value following '-M' flag" << std::endl;
00165                             usage( argv[0] );
00166                         }
00167                         merge_vertices = true;
00168                         ++i;
00169                         break;
00170                     default:
00171                         std::cerr << "Unrecognized flag: '" << argv[f][j] << "'" << std::endl;
00172                         usage( argv[0] );
00173                         break;
00174                 }
00175             }
00176         }
00177         else if( input_file && output_file )
00178         {
00179             std::cerr << "Extra argument: " << argv[i] << std::endl;
00180             usage( argv[0] );
00181         }
00182         else if( input_file )
00183         {
00184             output_file = argv[i++];
00185         }
00186         else
00187         {
00188             input_file = argv[i++];
00189         }
00190     }
00191 
00192     if( !input_file )
00193     {
00194         std::cerr << "No input file specified" << std::endl;
00195         usage( argv[0] );
00196     }
00197 
00198     ErrorCode result;
00199     Core mbimpl;
00200     Interface* iface = &mbimpl;
00201 
00202     if( print_perf )
00203     {
00204         double tmp_time1, tmp_mem1;
00205         get_time_mem( tmp_time1, tmp_mem1 );
00206         std::cout << "Before reading: cpu time = " << tmp_time1 << ", memory = " << tmp_mem1 / 1.0e6 << "MB."
00207                   << std::endl;
00208     }
00209 
00210     // read input file
00211     result = iface->load_mesh( input_file );
00212     if( MB_SUCCESS != result )
00213     {
00214         std::cerr << "Failed to load \"" << input_file << "\"." << std::endl;
00215         return 2;
00216     }
00217     std::cerr << "Read \"" << input_file << "\"" << std::endl;
00218     if( print_perf )
00219     {
00220         double tmp_time2, tmp_mem2;
00221         get_time_mem( tmp_time2, tmp_mem2 );
00222         std::cout << "After reading: cpu time = " << tmp_time2 << ", memory = " << tmp_mem2 / 1.0e6 << "MB."
00223                   << std::endl;
00224     }
00225 
00226     if( merge_vertices )
00227     {
00228         if( merge_epsilon < 0.0 )
00229         {
00230             if( MB_SUCCESS != min_edge_length( *iface, merge_epsilon ) )
00231             {
00232                 std::cerr << "Error determining minimum edge length" << std::endl;
00233                 return 1;
00234             }
00235             merge_epsilon /= MIN_EDGE_LEN_DENOM;
00236         }
00237         if( MB_SUCCESS != merge_duplicate_vertices( *iface, merge_epsilon ) )
00238         {
00239             std::cerr << "Error merging duplicate vertices" << std::endl;
00240             return 1;
00241         }
00242     }
00243 
00244     // get entities of largest dimension
00245     int dim = 4;
00246     Range entities;
00247     while( entities.empty() && dim > 1 )
00248     {
00249         dim--;
00250         result = iface->get_entities_by_dimension( 0, dim, entities );CHKERROR( result );
00251     }
00252 
00253     Range skin_ents;
00254     Tag matset_tag = 0, neuset_tag = 0;
00255     result = iface->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matset_tag );
00256     if( MB_SUCCESS != result ) return 1;
00257     result = iface->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, neuset_tag );
00258     if( MB_SUCCESS != result ) return 1;
00259 
00260     if( matsets.empty() )
00261         skin_ents = entities;
00262     else
00263     {
00264         // get all entities in the specified blocks
00265         if( 0 == matset_tag )
00266         {
00267             std::cerr << "Couldn't find any material sets in this mesh." << std::endl;
00268             return 1;
00269         }
00270 
00271         for( std::vector< int >::iterator vit = matsets.begin(); vit != matsets.end(); ++vit )
00272         {
00273             int this_matset             = *vit;
00274             const void* this_matset_ptr = &this_matset;
00275             Range this_range, ent_range;
00276             result =
00277                 iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &matset_tag, &this_matset_ptr, 1, this_range );
00278             if( MB_SUCCESS != result )
00279             {
00280                 std::cerr << "Trouble getting material set #" << *vit << std::endl;
00281                 return 1;
00282             }
00283             else if( this_range.empty() )
00284             {
00285                 std::cerr << "Warning: couldn't find material set " << *vit << std::endl;
00286                 continue;
00287             }
00288 
00289             result = iface->get_entities_by_dimension( *this_range.begin(), dim, ent_range, true );
00290             if( MB_SUCCESS != result ) continue;
00291             skin_ents.merge( ent_range );
00292         }
00293     }
00294 
00295     if( skin_ents.empty() )
00296     {
00297         std::cerr << "No entities for which to compute skin; exiting." << std::endl;
00298         return 1;
00299     }
00300 
00301     if( use_vert_elem_adjs )
00302     {
00303         // make a call which we know will generate vert-elem adjs
00304         Range dum_range;
00305         result = iface->get_adjacencies( &( *skin_ents.begin() ), 1, 1, false, dum_range );
00306         if( MB_SUCCESS != result ) return 1;
00307     }
00308 
00309     double tmp_time = 0.0, tmp_mem = 0.0;
00310     if( print_perf )
00311     {
00312         get_time_mem( tmp_time, tmp_mem );
00313         std::cout << "Before skinning: cpu time = " << tmp_time << ", memory = " << tmp_mem / 1.0e6 << "MB."
00314                   << std::endl;
00315     }
00316 
00317     // skin the mesh
00318     Range forward_lower, reverse_lower;
00319     Skinner tool( iface );
00320     if( use_scd )
00321         result = tool.find_skin( 0, skin_ents, false, forward_lower, NULL, false, true, true );
00322     else
00323         result = tool.find_skin( 0, skin_ents, false, forward_lower, &reverse_lower );
00324     Range boundary;
00325     boundary.merge( forward_lower );
00326     boundary.merge( reverse_lower );
00327     if( MB_SUCCESS != result || boundary.empty() )
00328     {
00329         std::cerr << "Mesh skinning failed." << std::endl;
00330         return 3;
00331     }
00332 
00333     if( list_skin )
00334     {
00335         Range skin_verts;
00336         result = iface->get_adjacencies( boundary, 0, true, skin_verts, Interface::UNION );
00337         std::cout << "Skin has ";
00338         if( skin_ents.num_of_dimension( 3 ) )
00339             std::cout << boundary.num_of_dimension( 2 ) << " faces and ";
00340         else if( skin_ents.num_of_dimension( 2 ) )
00341             std::cout << boundary.num_of_dimension( 1 ) << " edges and ";
00342         std::cout << skin_verts.size() << " vertices." << std::endl;
00343     }
00344     if( write_tag )
00345     {
00346         // get tag handle
00347         Tag tag;
00348         int zero = 0;
00349         result   = iface->tag_get_handle( fixed_tag, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );CHKERROR( result );
00350 
00351         // Set tags
00352         std::vector< int > ones;
00353         Range bverts;
00354         result = iface->get_adjacencies( boundary, 0, false, bverts, Interface::UNION );
00355         if( MB_SUCCESS != result )
00356         {
00357             std::cerr << "Trouble getting vertices on boundary." << std::endl;
00358             return 1;
00359         }
00360         ones.resize( bverts.size(), 1 );
00361         result = iface->tag_set_data( tag, bverts, &ones[0] );CHKERROR( result );
00362     }
00363 
00364     if( -1 != neuset_num )
00365     {
00366         // create a neumann set with these entities
00367         if( 0 == neuset_tag )
00368         {
00369             result = iface->tag_get_handle( "NEUMANN_SET_TAG_NAME", 1, MB_TYPE_INTEGER, neuset_tag,
00370                                             MB_TAG_SPARSE | MB_TAG_CREAT );
00371             if( MB_SUCCESS != result || 0 == neuset_tag ) return 1;
00372         }
00373 
00374         // always create a forward neumann set, assuming we have something in the set
00375         EntityHandle forward_neuset = 0;
00376         result                      = iface->create_meshset( MESHSET_SET, forward_neuset );
00377         if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
00378         result = iface->tag_set_data( neuset_tag, &forward_neuset, 1, &neuset_num );
00379         if( MB_SUCCESS != result ) return 1;
00380 
00381         if( !forward_lower.empty() )
00382         {
00383             result = iface->add_entities( forward_neuset, forward_lower );
00384             if( MB_SUCCESS != result ) return 1;
00385         }
00386         if( !reverse_lower.empty() )
00387         {
00388             EntityHandle reverse_neuset = 1;
00389             result                      = iface->create_meshset( MESHSET_SET, reverse_neuset );
00390             if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
00391 
00392             result = iface->add_entities( reverse_neuset, reverse_lower );
00393             if( MB_SUCCESS != result ) return 1;
00394             Tag sense_tag;
00395             int dum_sense = 0;
00396             result        = iface->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
00397                                                    &dum_sense );
00398             if( result != MB_SUCCESS ) return 1;
00399             int sense_val = -1;
00400             result        = iface->tag_set_data( neuset_tag, &reverse_neuset, 1, &sense_val );
00401             if( MB_SUCCESS != result ) return 0;
00402             result = iface->add_entities( forward_neuset, &reverse_neuset, 1 );
00403             if( MB_SUCCESS != result ) return 0;
00404         }
00405     }
00406 
00407     if( NULL != output_file && write_whole_mesh )
00408     {
00409 
00410         // write output file
00411         result = iface->write_mesh( output_file );
00412         if( MB_SUCCESS != result )
00413         {
00414             std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
00415             return 2;
00416         }
00417         std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
00418     }
00419     else if( NULL != output_file )
00420     {
00421         // write only skin; write them as one set
00422         EntityHandle skin_set;
00423         result = iface->create_meshset( MESHSET_SET, skin_set );
00424         if( MB_SUCCESS != result ) return 1;
00425         result = iface->add_entities( skin_set, forward_lower );
00426         if( MB_SUCCESS != result ) return 1;
00427         result = iface->add_entities( skin_set, reverse_lower );
00428         if( MB_SUCCESS != result ) return 1;
00429 
00430         int dum = 10000;
00431         result  = iface->tag_set_data( matset_tag, &skin_set, 1, &dum );
00432         if( MB_SUCCESS != result ) return 1;
00433 
00434         result = iface->write_mesh( output_file, &skin_set, 1 );
00435         if( MB_SUCCESS != result )
00436         {
00437             std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
00438             return 2;
00439         }
00440         std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
00441     }
00442 
00443     if( print_perf )
00444     {
00445         double tot_time, tot_mem;
00446         get_time_mem( tot_time, tot_mem );
00447         std::cout << "Total cpu time = " << tot_time << " seconds." << std::endl;
00448         std::cout << "Total skin cpu time = " << tot_time - tmp_time << " seconds." << std::endl;
00449         std::cout << "Total memory = " << tot_mem / 1024 << " MB." << std::endl;
00450         std::cout << "Total skin memory = " << ( tot_mem - tmp_mem ) / 1024 << " MB." << std::endl;
00451         std::cout << "Entities: " << std::endl;
00452         iface->list_entities( 0, 0 );
00453     }
00454 
00455     return 0;
00456 }
00457 
00458 #if defined( _MSC_VER ) || defined( __MINGW32__ )
00459 void get_time_mem( double& tot_time, double& tot_mem )
00460 {
00461     tot_time = (double)clock() / CLOCKS_PER_SEC;
00462     tot_mem  = 0;
00463 }
00464 #else
00465 void get_time_mem( double& tot_time, double& tot_mem )
00466 {
00467     struct rusage r_usage;
00468     getrusage( RUSAGE_SELF, &r_usage );
00469     double utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
00470     double stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
00471     tot_time     = utime + stime;
00472     tot_mem      = 0;
00473     if( 0 != r_usage.ru_maxrss )
00474     {
00475         tot_mem = (double)r_usage.ru_maxrss;
00476     }
00477     else
00478     {
00479         // this machine doesn't return rss - try going to /proc
00480         // print the file name to open
00481         char file_str[4096], dum_str[4096];
00482         int file_ptr = open( "/proc/self/stat", O_RDONLY );
00483         int file_len = read( file_ptr, file_str, sizeof( file_str ) - 1 );
00484         if( file_len <= 0 )
00485         {
00486             close( file_ptr );
00487             return;
00488         }
00489 
00490         close( file_ptr );
00491         file_str[file_len] = '\0';
00492         // read the preceding fields and the ones we really want...
00493         int dum_int;
00494         unsigned int dum_uint, vm_size, rss;
00495         int num_fields = sscanf( file_str,
00496                                  "%d "                 // pid
00497                                  "%s "                 // comm
00498                                  "%c "                 // state
00499                                  "%d %d %d %d %d "     // ppid, pgrp, session, tty, tpgid
00500                                  "%u %u %u %u %u "     // flags, minflt, cminflt, majflt, cmajflt
00501                                  "%d %d %d %d %d %d "  // utime, stime, cutime, cstime, counter, priority
00502                                  "%u %u "              // timeout, itrealvalue
00503                                  "%d "                 // starttime
00504                                  "%u %u",              // vsize, rss
00505                                  &dum_int, dum_str, dum_str, &dum_int, &dum_int, &dum_int, &dum_int, &dum_int,
00506                                  &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_int, &dum_int, &dum_int,
00507                                  &dum_int, &dum_int, &dum_int, &dum_uint, &dum_uint, &dum_int, &vm_size, &rss );
00508         if( num_fields == 24 ) tot_mem = ( (double)vm_size );
00509     }
00510 }
00511 #endif
00512 
00513 ErrorCode min_edge_length( Interface& moab, double& result )
00514 {
00515     double sqr_result = std::numeric_limits< double >::max();
00516 
00517     ErrorCode rval;
00518     Range entities;
00519     rval = moab.get_entities_by_handle( 0, entities );
00520     if( MB_SUCCESS != rval ) return rval;
00521     Range::iterator i = entities.upper_bound( MBVERTEX );
00522     entities.erase( entities.begin(), i );
00523     i = entities.lower_bound( MBENTITYSET );
00524     entities.erase( i, entities.end() );
00525 
00526     std::vector< EntityHandle > storage;
00527     for( i = entities.begin(); i != entities.end(); ++i )
00528     {
00529         EntityType t = moab.type_from_handle( *i );
00530         const EntityHandle* conn;
00531         int conn_len, indices[2];
00532         rval = moab.get_connectivity( *i, conn, conn_len, true, &storage );
00533         if( MB_SUCCESS != rval ) return rval;
00534 
00535         int num_edges = CN::NumSubEntities( t, 1 );
00536         for( int j = 0; j < num_edges; ++j )
00537         {
00538             CN::SubEntityVertexIndices( t, 1, j, indices );
00539             EntityHandle v[2] = { conn[indices[0]], conn[indices[1]] };
00540             if( v[0] == v[1] ) continue;
00541 
00542             double c[6];
00543             rval = moab.get_coords( v, 2, c );
00544             if( MB_SUCCESS != rval ) return rval;
00545 
00546             c[0] -= c[3];
00547             c[1] -= c[4];
00548             c[2] -= c[5];
00549             double len_sqr = c[0] * c[0] + c[1] * c[1] + c[2] * c[2];
00550             if( len_sqr < sqr_result ) sqr_result = len_sqr;
00551         }
00552     }
00553 
00554     result = sqrt( sqr_result );
00555     return MB_SUCCESS;
00556 }
00557 
00558 ErrorCode merge_duplicate_vertices( Interface& moab, const double epsilon )
00559 {
00560     ErrorCode rval;
00561     Range verts;
00562     rval = moab.get_entities_by_type( 0, MBVERTEX, verts );
00563     if( MB_SUCCESS != rval ) return rval;
00564 
00565     AdaptiveKDTree tree( &moab );
00566     EntityHandle root;
00567     rval = tree.build_tree( verts, &root );
00568     if( MB_SUCCESS != rval )
00569     {
00570         fprintf( stderr, "Failed to build kD-tree.\n" );
00571         return rval;
00572     }
00573 
00574     std::set< EntityHandle > dead_verts;
00575     std::vector< EntityHandle > leaves;
00576     for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
00577     {
00578         double coords[3];
00579         rval = moab.get_coords( &*i, 1, coords );
00580         if( MB_SUCCESS != rval ) return rval;
00581 
00582         leaves.clear();
00583         ;
00584         rval = tree.distance_search( coords, epsilon, leaves, epsilon, epsilon );
00585         if( MB_SUCCESS != rval ) return rval;
00586 
00587         Range near;
00588         for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
00589         {
00590             Range tmp;
00591             rval = moab.get_entities_by_type( *j, MBVERTEX, tmp );
00592             if( MB_SUCCESS != rval ) return rval;
00593             near.merge( tmp.begin(), tmp.end() );
00594         }
00595 
00596         Range::iterator v = near.find( *i );
00597         assert( v != near.end() );
00598         near.erase( v );
00599 
00600         EntityHandle merge = 0;
00601         for( Range::iterator j = near.begin(); j != near.end(); ++j )
00602         {
00603             if( *j < *i && dead_verts.find( *j ) != dead_verts.end() ) continue;
00604 
00605             double coords2[3];
00606             rval = moab.get_coords( &*j, 1, coords2 );
00607             if( MB_SUCCESS != rval ) return rval;
00608 
00609             coords2[0] -= coords[0];
00610             coords2[1] -= coords[1];
00611             coords2[2] -= coords[2];
00612             double dsqr = coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2];
00613             if( dsqr <= epsilon * epsilon )
00614             {
00615                 merge = *j;
00616                 break;
00617             }
00618         }
00619 
00620         if( merge )
00621         {
00622             dead_verts.insert( *i );
00623             rval = moab.merge_entities( merge, *i, false, true );
00624             if( MB_SUCCESS != rval ) return rval;
00625         }
00626     }
00627 
00628     if( dead_verts.empty() )
00629         std::cout << "No duplicate/coincident vertices." << std::endl;
00630     else
00631         std::cout << "Merged and deleted " << dead_verts.size() << " vertices "
00632                   << "coincident within " << epsilon << std::endl;
00633 
00634     return MB_SUCCESS;
00635 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines