skin.cpp

Go to the documentation of this file.

#include <iostream>
#include <set>
#include <limits>
#include <ctime>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <cassert>
#if !defined( _MSC_VER ) && !defined( __MINGW32__ )
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#endif
#include <fcntl.h>
#include "moab/Interface.hpp"
#include "MBTagConventions.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/Skinner.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/CN.hpp"
 
using namespace moab;
 
static void get_time_mem( double& tot_time, double& tot_mem );
 
// Different platforms follow different conventions for usage
#if !defined( _MSC_VER ) && !defined( __MINGW32__ )
#include <sys/resource.h>
#endif
#ifdef SOLARIS
extern "C" int getrusage( int, struct rusage* );
#ifndef RUSAGE_SELF
#include </usr/ucbinclude/sys/rusage.h>
#endif
#endif
 
const char DEFAULT_FIXED_TAG[] = "fixed";
const int MIN_EDGE_LEN_DENOM = 4;
 
#define CHKERROR( A ) \
 do \
 { \
 if( MB_SUCCESS != ( A ) ) \
 { \
 std::cerr << "Internal error at line " << __LINE__ << std::endl; \
 return 3; \
 } \
 } while( false )
 
static ErrorCode merge_duplicate_vertices( Interface&, double epsilon );
static ErrorCode min_edge_length( Interface&, double& result );
 
static void usage( const char* argv0, bool help = false )
{
 std::ostream& str = help ? std::cout : std::cerr;
 
 str << "Usage: " << argv0
 << " [-b <block_num> [-b ...] ] [-l] [-m] [-M <n>] [-p] [-s <sideset_num>] [-S] [-t|-T "
 "<name>] [-w] [-v|-V <n>]"
 << " <input_file> [<output_file>]" << std::endl;
 str << "Help : " << argv0 << " -h" << std::endl;
 if( !help ) exit( 1 );
 
 str << "Options: " << std::endl;
 str << "-a : Compute skin using vert-elem adjacencies (more memory, less time)." << std::endl;
 str << "-b <block_num> : Compute skin only for material set/block <block_num>." << std::endl;
 str << "-p : Print cpu & memory performance." << std::endl;
 str << "-s <sideset_num> : Put skin in neumann set/sideset <sideset_num>." << std::endl;
 str << "-S : Look for and use structured mesh information to speed up skinning." << std::endl;
 str << "-t : Set '" << DEFAULT_FIXED_TAG << "' tag on skin vertices." << std::endl;
 str << "-T <name> : Create tag with specified name and set to 1 on skin vertices." << std::endl;
 str << "-w : Write out whole mesh (otherwise just writes skin)." << std::endl;
 str << "-m : consolidate duplicate vertices" << std::endl;
 str << "-M <n> : consolidate duplicate vertices with specified tolerance. "
 "(Default: min_edge_length/"
 << MIN_EDGE_LEN_DENOM << ")" << std::endl;
 str << "-l : List total numbers of entities and vertices in skin." << std::endl;
 exit( 0 );
}
 
int main( int argc, char* argv[] )
{
 int i = 1;
 std::vector< int > matsets;
 int neuset_num = -1;
 bool write_tag = false, write_whole_mesh = false;
 bool print_perf = false;
 bool use_vert_elem_adjs = false;
 bool merge_vertices = false;
 double merge_epsilon = -1;
 bool list_skin = false;
 bool use_scd = false;
 const char* fixed_tag = DEFAULT_FIXED_TAG;
 const char *input_file = 0, *output_file = 0;
 
 bool no_more_flags = false;
 char* endptr = 0;
 long block = 0; // initialize to eliminate compiler warning
 while( i < argc )
 {
 if( !no_more_flags && argv[i][0] == '-' )
 {
 const int f = i++;
 for( int j = 1; argv[f][j]; ++j )
 {
 switch( argv[f][j] )
 {
 case 'a':
 use_vert_elem_adjs = true;
 break;
 case 'p':
 print_perf = true;
 break;
 case 't':
 write_tag = true;
 break;
 case 'w':
 write_whole_mesh = true;
 break;
 case 'm':
 merge_vertices = true;
 break;
 case '-':
 no_more_flags = true;
 break;
 case 'h':
 usage( argv[0], true );
 break;
 case 'l':
 list_skin = true;
 break;
 case 'S':
 use_scd = true;
 break;
 case 'b':
 if( i == argc || 0 >= ( block = strtol( argv[i], &endptr, 0 ) ) || *endptr )
 {
 std::cerr << "Expected positive integer following '-b' flag" << std::endl;
 usage( argv[0] );
 }
 matsets.push_back( (int)block );
 ++i;
 break;
 case 's':
 if( i == argc || 0 >= ( neuset_num = strtol( argv[i], &endptr, 0 ) ) || *endptr )
 {
 std::cerr << "Expected positive integer following '-s' flag" << std::endl;
 usage( argv[0] );
 }
 ++i;
 break;
 case 'T':
 if( i == argc || argv[i][0] == '-' )
 {
 std::cerr << "Expected tag name following '-T' flag" << std::endl;
 usage( argv[0] );
 }
 fixed_tag = argv[i++];
 break;
 case 'M':
 if( i == argc || 0.0 > ( merge_epsilon = strtod( argv[i], &endptr ) ) || *endptr )
 {
 std::cerr << "Expected positive numeric value following '-M' flag" << std::endl;
 usage( argv[0] );
 }
 merge_vertices = true;
 ++i;
 break;
 default:
 std::cerr << "Unrecognized flag: '" << argv[f][j] << "'" << std::endl;
 usage( argv[0] );
 break;
 }
 }
 }
 else if( input_file && output_file )
 {
 std::cerr << "Extra argument: " << argv[i] << std::endl;
 usage( argv[0] );
 }
 else if( input_file )
 {
 output_file = argv[i++];
 }
 else
 {
 input_file = argv[i++];
 }
 }
 
 if( !input_file )
 {
 std::cerr << "No input file specified" << std::endl;
 usage( argv[0] );
 }
 
 ErrorCode result;
 Core mbimpl;
 Interface* iface = &mbimpl;
 
 if( print_perf )
 {
 double tmp_time1, tmp_mem1;
 get_time_mem( tmp_time1, tmp_mem1 );
 std::cout << "Before reading: cpu time = " << tmp_time1 << ", memory = " << tmp_mem1 / 1.0e6 << "MB."
 << std::endl;
 }
 
 // read input file
 result = iface->load_mesh( input_file );
 if( MB_SUCCESS != result )
 {
 std::cerr << "Failed to load \"" << input_file << "\"." << std::endl;
 return 2;
 }
 std::cerr << "Read \"" << input_file << "\"" << std::endl;
 if( print_perf )
 {
 double tmp_time2, tmp_mem2;
 get_time_mem( tmp_time2, tmp_mem2 );
 std::cout << "After reading: cpu time = " << tmp_time2 << ", memory = " << tmp_mem2 / 1.0e6 << "MB."
 << std::endl;
 }
 
 if( merge_vertices )
 {
 if( merge_epsilon < 0.0 )
 {
 if( MB_SUCCESS != min_edge_length( *iface, merge_epsilon ) )
 {
 std::cerr << "Error determining minimum edge length" << std::endl;
 return 1;
 }
 merge_epsilon /= MIN_EDGE_LEN_DENOM;
 }
 if( MB_SUCCESS != merge_duplicate_vertices( *iface, merge_epsilon ) )
 {
 std::cerr << "Error merging duplicate vertices" << std::endl;
 return 1;
 }
 }
 
 // get entities of largest dimension
 int dim = 4;
 Range entities;
 while( entities.empty() && dim > 1 )
 {
 dim--;
 result = iface->get_entities_by_dimension( 0, dim, entities );CHKERROR( result );
 }
 
 Range skin_ents;
 Tag matset_tag = 0, neuset_tag = 0;
 result = iface->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matset_tag );
 if( MB_SUCCESS != result ) return 1;
 result = iface->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, neuset_tag );
 if( MB_SUCCESS != result ) return 1;
 
 if( matsets.empty() )
 skin_ents = entities;
 else
 {
 // get all entities in the specified blocks
 if( 0 == matset_tag )
 {
 std::cerr << "Couldn't find any material sets in this mesh." << std::endl;
 return 1;
 }
 
 for( std::vector< int >::iterator vit = matsets.begin(); vit != matsets.end(); ++vit )
 {
 int this_matset = *vit;
 const void* this_matset_ptr = &this_matset;
 Range this_range, ent_range;
 result =
 iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &matset_tag, &this_matset_ptr, 1, this_range );
 if( MB_SUCCESS != result )
 {
 std::cerr << "Trouble getting material set #" << *vit << std::endl;
 return 1;
 }
 else if( this_range.empty() )
 {
 std::cerr << "Warning: couldn't find material set " << *vit << std::endl;
 continue;
 }
 
 result = iface->get_entities_by_dimension( *this_range.begin(), dim, ent_range, true );
 if( MB_SUCCESS != result ) continue;
 skin_ents.merge( ent_range );
 }
 }
 
 if( skin_ents.empty() )
 {
 std::cerr << "No entities for which to compute skin; exiting." << std::endl;
 return 1;
 }
 
 if( use_vert_elem_adjs )
 {
 // make a call which we know will generate vert-elem adjs
 Range dum_range;
 result = iface->get_adjacencies( &( *skin_ents.begin() ), 1, 1, false, dum_range );
 if( MB_SUCCESS != result ) return 1;
 }
 
 double tmp_time = 0.0, tmp_mem = 0.0;
 if( print_perf )
 {
 get_time_mem( tmp_time, tmp_mem );
 std::cout << "Before skinning: cpu time = " << tmp_time << ", memory = " << tmp_mem / 1.0e6 << "MB."
 << std::endl;
 }
 
 // skin the mesh
 Range forward_lower, reverse_lower;
 Skinner tool( iface );
 if( use_scd )
 result = tool.find_skin( 0, skin_ents, false, forward_lower, NULL, false, true, true );
 else
 result = tool.find_skin( 0, skin_ents, false, forward_lower, &reverse_lower );
 Range boundary;
 boundary.merge( forward_lower );
 boundary.merge( reverse_lower );
 if( MB_SUCCESS != result || boundary.empty() )
 {
 std::cerr << "Mesh skinning failed." << std::endl;
 return 3;
 }
 
 if( list_skin )
 {
 Range skin_verts;
 result = iface->get_adjacencies( boundary, 0, true, skin_verts, Interface::UNION );
 std::cout << "Skin has ";
 if( skin_ents.num_of_dimension( 3 ) )
 std::cout << boundary.num_of_dimension( 2 ) << " faces and ";
 else if( skin_ents.num_of_dimension( 2 ) )
 std::cout << boundary.num_of_dimension( 1 ) << " edges and ";
 std::cout << skin_verts.size() << " vertices." << std::endl;
 }
 if( write_tag )
 {
 // get tag handle
 Tag tag;
 int zero = 0;
 result = iface->tag_get_handle( fixed_tag, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );CHKERROR( result );
 
 // Set tags
 std::vector< int > ones;
 Range bverts;
 result = iface->get_adjacencies( boundary, 0, false, bverts, Interface::UNION );
 if( MB_SUCCESS != result )
 {
 std::cerr << "Trouble getting vertices on boundary." << std::endl;
 return 1;
 }
 ones.resize( bverts.size(), 1 );
 result = iface->tag_set_data( tag, bverts, &ones[0] );CHKERROR( result );
 }
 
 if( -1 != neuset_num )
 {
 // create a neumann set with these entities
 if( 0 == neuset_tag )
 {
 result = iface->tag_get_handle( "NEUMANN_SET_TAG_NAME", 1, MB_TYPE_INTEGER, neuset_tag,
 MB_TAG_SPARSE | MB_TAG_CREAT );
 if( MB_SUCCESS != result || 0 == neuset_tag ) return 1;
 }
 
 // always create a forward neumann set, assuming we have something in the set
 EntityHandle forward_neuset = 0;
 result = iface->create_meshset( MESHSET_SET, forward_neuset );
 if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
 result = iface->tag_set_data( neuset_tag, &forward_neuset, 1, &neuset_num );
 if( MB_SUCCESS != result ) return 1;
 
 if( !forward_lower.empty() )
 {
 result = iface->add_entities( forward_neuset, forward_lower );
 if( MB_SUCCESS != result ) return 1;
 }
 if( !reverse_lower.empty() )
 {
 EntityHandle reverse_neuset = 1;
 result = iface->create_meshset( MESHSET_SET, reverse_neuset );
 if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
 
 result = iface->add_entities( reverse_neuset, reverse_lower );
 if( MB_SUCCESS != result ) return 1;
 Tag sense_tag;
 int dum_sense = 0;
 result = iface->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
 &dum_sense );
 if( result != MB_SUCCESS ) return 1;
 int sense_val = -1;
 result = iface->tag_set_data( neuset_tag, &reverse_neuset, 1, &sense_val );
 if( MB_SUCCESS != result ) return 0;
 result = iface->add_entities( forward_neuset, &reverse_neuset, 1 );
 if( MB_SUCCESS != result ) return 0;
 }
 }
 
 if( NULL != output_file && write_whole_mesh )
 {
 
 // write output file
 result = iface->write_mesh( output_file );
 if( MB_SUCCESS != result )
 {
 std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
 return 2;
 }
 std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
 }
 else if( NULL != output_file )
 {
 // write only skin; write them as one set
 EntityHandle skin_set;
 result = iface->create_meshset( MESHSET_SET, skin_set );
 if( MB_SUCCESS != result ) return 1;
 result = iface->add_entities( skin_set, forward_lower );
 if( MB_SUCCESS != result ) return 1;
 result = iface->add_entities( skin_set, reverse_lower );
 if( MB_SUCCESS != result ) return 1;
 
 int dum = 10000;
 result = iface->tag_set_data( matset_tag, &skin_set, 1, &dum );
 if( MB_SUCCESS != result ) return 1;
 
 result = iface->write_mesh( output_file, &skin_set, 1 );
 if( MB_SUCCESS != result )
 {
 std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
 return 2;
 }
 std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
 }
 
 if( print_perf )
 {
 double tot_time, tot_mem;
 get_time_mem( tot_time, tot_mem );
 std::cout << "Total cpu time = " << tot_time << " seconds." << std::endl;
 std::cout << "Total skin cpu time = " << tot_time - tmp_time << " seconds." << std::endl;
 std::cout << "Total memory = " << tot_mem / 1024 << " MB." << std::endl;
 std::cout << "Total skin memory = " << ( tot_mem - tmp_mem ) / 1024 << " MB." << std::endl;
 std::cout << "Entities: " << std::endl;
 iface->list_entities( 0, 0 );
 }
 
 return 0;
}
 
#if defined( _MSC_VER ) || defined( __MINGW32__ )
void get_time_mem( double& tot_time, double& tot_mem )
{
 tot_time = (double)clock() / CLOCKS_PER_SEC;
 tot_mem = 0;
}
#else
void get_time_mem( double& tot_time, double& tot_mem )
{
 struct rusage r_usage;
 getrusage( RUSAGE_SELF, &r_usage );
 double utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
 double stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
 tot_time = utime + stime;
 tot_mem = 0;
 if( 0 != r_usage.ru_maxrss )
 {
 tot_mem = (double)r_usage.ru_maxrss;
 }
 else
 {
 // this machine doesn't return rss - try going to /proc
 // print the file name to open
 char file_str[4096], dum_str[4096];
 int file_ptr = open( "/proc/self/stat", O_RDONLY );
 int file_len = read( file_ptr, file_str, sizeof( file_str ) - 1 );
 if( file_len <= 0 )
 {
 close( file_ptr );
 return;
 }
 
 close( file_ptr );
 file_str[file_len] = '\0';
 // read the preceding fields and the ones we really want...
 int dum_int;
 unsigned int dum_uint, vm_size, rss;
 int num_fields = sscanf( file_str,
 "%d " // pid
 "%s " // comm
 "%c " // state
 "%d %d %d %d %d " // ppid, pgrp, session, tty, tpgid
 "%u %u %u %u %u " // flags, minflt, cminflt, majflt, cmajflt
 "%d %d %d %d %d %d " // utime, stime, cutime, cstime, counter, priority
 "%u %u " // timeout, itrealvalue
 "%d " // starttime
 "%u %u", // vsize, rss
 &dum_int, dum_str, dum_str, &dum_int, &dum_int, &dum_int, &dum_int, &dum_int,
 &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_int, &dum_int, &dum_int,
 &dum_int, &dum_int, &dum_int, &dum_uint, &dum_uint, &dum_int, &vm_size, &rss );
 if( num_fields == 24 ) tot_mem = ( (double)vm_size );
 }
}
#endif
 
ErrorCode min_edge_length( Interface& moab, double& result )
{
 double sqr_result = std::numeric_limits< double >::max();
 
 ErrorCode rval;
 Range entities;
 rval = moab.get_entities_by_handle( 0, entities );
 if( MB_SUCCESS != rval ) return rval;
 Range::iterator i = entities.upper_bound( MBVERTEX );
 entities.erase( entities.begin(), i );
 i = entities.lower_bound( MBENTITYSET );
 entities.erase( i, entities.end() );
 
 std::vector< EntityHandle > storage;
 for( i = entities.begin(); i != entities.end(); ++i )
 {
 EntityType t = moab.type_from_handle( *i );
 const EntityHandle* conn;
 int conn_len, indices[2];
 rval = moab.get_connectivity( *i, conn, conn_len, true, &storage );
 if( MB_SUCCESS != rval ) return rval;
 
 int num_edges = CN::NumSubEntities( t, 1 );
 for( int j = 0; j < num_edges; ++j )
 {
 CN::SubEntityVertexIndices( t, 1, j, indices );
 EntityHandle v[2] = { conn[indices[0]], conn[indices[1]] };
 if( v[0] == v[1] ) continue;
 
 double c[6];
 rval = moab.get_coords( v, 2, c );
 if( MB_SUCCESS != rval ) return rval;
 
 c[0] -= c[3];
 c[1] -= c[4];
 c[2] -= c[5];
 double len_sqr = c[0] * c[0] + c[1] * c[1] + c[2] * c[2];
 if( len_sqr < sqr_result ) sqr_result = len_sqr;
 }
 }
 
 result = sqrt( sqr_result );
 return MB_SUCCESS;
}
 
ErrorCode merge_duplicate_vertices( Interface& moab, const double epsilon )
{
 ErrorCode rval;
 Range verts;
 rval = moab.get_entities_by_type( 0, MBVERTEX, verts );
 if( MB_SUCCESS != rval ) return rval;
 
 AdaptiveKDTree tree( &moab );
 EntityHandle root;
 rval = tree.build_tree( verts, &root );
 if( MB_SUCCESS != rval )
 {
 fprintf( stderr, "Failed to build kD-tree.\n" );
 return rval;
 }
 
 std::set< EntityHandle > dead_verts;
 std::vector< EntityHandle > leaves;
 for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
 {
 double coords[3];
 rval = moab.get_coords( &*i, 1, coords );
 if( MB_SUCCESS != rval ) return rval;
 
 leaves.clear();
 ;
 rval = tree.distance_search( coords, epsilon, leaves, epsilon, epsilon );
 if( MB_SUCCESS != rval ) return rval;
 
 Range near;
 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
 {
 Range tmp;
 rval = moab.get_entities_by_type( *j, MBVERTEX, tmp );
 if( MB_SUCCESS != rval ) return rval;
 near.merge( tmp.begin(), tmp.end() );
 }
 
 Range::iterator v = near.find( *i );
 assert( v != near.end() );
 near.erase( v );
 
 EntityHandle merge = 0;
 for( Range::iterator j = near.begin(); j != near.end(); ++j )
 {
 if( *j < *i && dead_verts.find( *j ) != dead_verts.end() ) continue;
 
 double coords2[3];
 rval = moab.get_coords( &*j, 1, coords2 );
 if( MB_SUCCESS != rval ) return rval;
 
 coords2[0] -= coords[0];
 coords2[1] -= coords[1];
 coords2[2] -= coords[2];
 double dsqr = coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2];
 if( dsqr <= epsilon * epsilon )
 {
 merge = *j;
 break;
 }
 }
 
 if( merge )
 {
 dead_verts.insert( *i );
 rval = moab.merge_entities( merge, *i, false, true );
 if( MB_SUCCESS != rval ) return rval;
 }
 }
 
 if( dead_verts.empty() )
 std::cout << "No duplicate/coincident vertices." << std::endl;
 else
 std::cout << "Merged and deleted " << dead_verts.size() << " vertices "
 << "coincident within " << epsilon << std::endl;
 
 return MB_SUCCESS;
}