MOAB: Mesh Oriented datABase  (version 5.4.0)
obb_analysis.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <sstream>
00004 
00005 #include "moab/Core.hpp"
00006 #include "moab/Interface.hpp"
00007 #include "moab/OrientedBoxTreeTool.hpp"
00008 #include "moab/CartVect.hpp"
00009 #include "moab/OrientedBox.hpp"
00010 #include "moab/GeomTopoTool.hpp"
00011 
00012 #include "cgm2moab.hpp"
00013 
00014 using namespace moab;
00015 
00016 class TriCounter : public OrientedBoxTreeTool::Op
00017 {
00018 
00019   public:
00020     int count;
00021     Interface* mbi;
00022     OrientedBoxTreeTool* tool;
00023     const CartVect& pt;
00024 
00025     TriCounter( Interface* mbi_p, OrientedBoxTreeTool* tool_p, const CartVect& p )
00026         : OrientedBoxTreeTool::Op(), count( 0 ), mbi( mbi_p ), tool( tool_p ), pt( p )
00027     {
00028     }
00029 
00030     virtual ErrorCode visit( EntityHandle node, int, bool& descend )
00031     {
00032         OrientedBox box;
00033         ErrorCode rval = tool->box( node, box );
00034 
00035         descend = box.contained( pt, 1e-6 );
00036 
00037         return rval;
00038     }
00039 
00040     virtual ErrorCode leaf( EntityHandle node )
00041     {
00042 
00043         int numtris;
00044         ErrorCode rval = tool->get_moab_instance()->get_number_entities_by_type( node, MBTRI, numtris );
00045         count += numtris;
00046         return rval;
00047     }
00048 };
00049 
00050 ErrorCode obbvis_create( GeomTopoTool& gtt, std::vector< int >& volumes, int grid, std::string& filename )
00051 {
00052     OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
00053 
00054     CartVect min, max;
00055     EntityHandle vol = gtt.entity_by_id( 3, volumes.front() );
00056     double middle[3];
00057     double axis1[3], axis2[3], axis3[3];
00058     double minPt[3], maxPt[3];
00059     ErrorCode rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
00060     // compute min and max verticies
00061     for( int i = 0; i < 3; i++ )
00062     {
00063         double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
00064         minPt[i]   = middle[i] - sum;
00065         maxPt[i]   = middle[i] + sum;
00066     }
00067     CHECKERR( gtt, rval );
00068 
00069     /* Compute an axis-aligned bounding box of all the requested volumes */
00070     for( std::vector< int >::iterator i = volumes.begin() + 1; i != volumes.end(); ++i )
00071     {
00072         CartVect i_min, i_max;
00073         vol  = gtt.entity_by_id( 3, *i );
00074         rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
00075         // compute min and max verticies
00076         for( int j = 0; j < 3; j++ )
00077         {
00078             double sum = fabs( axis1[j] ) + fabs( axis2[j] ) + fabs( axis3[j] );
00079             minPt[j]   = middle[j] - sum;
00080             maxPt[j]   = middle[j] + sum;
00081         }
00082 
00083         for( int j = 0; j < 3; ++j )
00084         {
00085             min[j] = std::min( min[j], i_min[j] );
00086             max[j] = std::max( max[j], i_max[j] );
00087         }
00088     }
00089 
00090     // These vectors could be repurposed to describe an OBB without changing the loops below
00091     CartVect center( middle );
00092     CartVect v1( axis1 );
00093     CartVect v2( axis2 );
00094     CartVect v3( axis3 );
00095 
00096     /* Compute the vertices of the visualization grid.  Calculation points are at the center
00097        of each cell in this grid, so make grid+1 vertices in each direction. */
00098     int numpoints = pow( (double)( grid + 1 ), 3 );
00099     double* pgrid = new double[numpoints * 3];
00100     int idx       = 0;
00101 
00102     for( int i = 0; i < numpoints * 3; ++i )
00103         pgrid[i] = 0.0;
00104 
00105     for( int i = 0; i <= grid; ++i )
00106     {
00107         CartVect x = -v1 + ( ( v1 * 2.0 ) * ( i / (double)grid ) );
00108 
00109         for( int j = 0; j <= grid; ++j )
00110         {
00111             CartVect y = -v2 + ( ( v2 * 2.0 ) * ( j / (double)grid ) );
00112 
00113             for( int k = 0; k <= grid; ++k )
00114             {
00115                 CartVect z = -v3 + ( ( v3 * 2.0 ) * ( k / (double)grid ) );
00116 
00117                 CartVect p = center + x + y + z;
00118                 for( int d = 0; d < 3; ++d )
00119                 {
00120                     pgrid[idx++] = p[d];
00121                 }
00122             }
00123         }
00124     }
00125 
00126     /* Create a new MOAB to use for output, and build the vertex grid */
00127     Core mb2;
00128     Range r;
00129     rval = mb2.create_vertices( pgrid, numpoints, r );
00130     CHECKERR( mb2, rval );
00131 
00132     Tag lttag;
00133     rval = mb2.tag_get_handle( "LEAFTRIS", sizeof( int ), MB_TYPE_INTEGER, lttag,
00134                                MB_TAG_EXCL | MB_TAG_CREAT | MB_TAG_BYTES | MB_TAG_DENSE, 0 );
00135     CHECKERR( mb2, rval );
00136 
00137     int row  = grid + 1;
00138     int side = row * row;
00139     EntityHandle connect[8];
00140     EntityHandle hex;
00141 
00142     // offset from grid corner to grid center
00143     CartVect grid_hex_center_offset = ( v1 + v2 + v3 ) * 2 * ( 1.0 / grid );
00144 
00145     // collect all the surfaces from the requested volumes to iterate over --
00146     // this prevents checking a shared surface more than once.
00147     Range surfs;
00148     for( std::vector< int >::iterator it = volumes.begin(); it != volumes.end(); ++it )
00149     {
00150 
00151         vol = gtt.entity_by_id( 3, *it );
00152         Range it_surfs;
00153         rval = gtt.get_moab_instance()->get_child_meshsets( vol, it_surfs );
00154         CHECKERR( gtt, rval );
00155         surfs.merge( it_surfs );
00156     }
00157     std::cout << "visualizing " << surfs.size() << " surfaces." << std::endl;
00158 
00159     /* Build hexes for any point with 1 or more leaftris */
00160     for( int i = 0; i < grid; ++i )
00161     {
00162         for( int j = 0; j < grid; ++j )
00163         {
00164             for( int k = 0; k < grid; ++k )
00165             {
00166 
00167                 idx = ( side * k ) + ( row * j ) + i;
00168                 assert( idx + 1 + row + side > numpoints - 1 );
00169 
00170                 CartVect loc = CartVect( ( pgrid + ( idx * 3 ) ) ) + grid_hex_center_offset;
00171                 TriCounter tc( gtt.get_moab_instance(), &obbtool, loc );
00172 
00173                 for( Range::iterator it = surfs.begin(); it != surfs.end(); ++it )
00174                 {
00175 
00176                     EntityHandle surf_tree;
00177                     rval = gtt.get_root( *it, surf_tree );
00178                     CHECKERR( gtt, rval );
00179 
00180                     rval = obbtool.preorder_traverse( surf_tree, tc );
00181                     CHECKERR( gtt, rval );
00182                 }
00183 
00184                 if( tc.count == 0 ) continue;
00185 
00186                 connect[0] = r[idx];
00187                 connect[1] = r[idx + 1];
00188                 connect[2] = r[idx + 1 + row];
00189                 connect[3] = r[idx + row];
00190                 connect[4] = r[idx + side];
00191                 connect[5] = r[idx + 1 + side];
00192                 connect[6] = r[idx + 1 + row + side];
00193                 connect[7] = r[idx + row + side];
00194 
00195                 rval = mb2.create_element( MBHEX, connect, 8, hex );
00196                 CHECKERR( mb2, rval );
00197 
00198                 rval = mb2.tag_set_data( lttag, &hex, 1, &( tc.count ) );
00199                 CHECKERR( mb2, rval );
00200             }
00201         }
00202     }
00203 
00204     if( verbose )
00205     {
00206         std::cout << "Writing " << filename << std::endl;
00207     }
00208     rval = mb2.write_file( filename.c_str() );
00209     CHECKERR( mb2, rval );
00210 
00211     return rval;
00212 }
00213 
00214 // stats code borrowed from OrientedBoxTreeTool.cpp
00215 static inline double std_dev( double sqr, double sum, double count )
00216 {
00217     sum /= count;
00218     sqr /= count;
00219     return sqrt( sqr - sum * sum );
00220 }
00221 
00222 class TriStats : public OrientedBoxTreeTool::Op
00223 {
00224 
00225   public:
00226     unsigned min, max, sum, leaves;
00227     double sqr;
00228 
00229     const static unsigned ten_buckets_max = 5;
00230     unsigned ten_buckets[ten_buckets_max];
00231     double ten_buckets_vol[ten_buckets_max];
00232 
00233     Interface* mbi;
00234     OrientedBoxTreeTool* tool;
00235 
00236     double tot_vol;
00237 
00238     TriStats( Interface* mbi_p, OrientedBoxTreeTool* tool_p, EntityHandle root )
00239         : OrientedBoxTreeTool::Op(), sum( 0 ), leaves( 0 ), sqr( 0 ), mbi( mbi_p ), tool( tool_p )
00240     {
00241         min = std::numeric_limits< unsigned >::max();
00242         max = std::numeric_limits< unsigned >::min();
00243 
00244         for( unsigned i = 0; i < ten_buckets_max; ++i )
00245         {
00246             ten_buckets[i]     = 0;
00247             ten_buckets_vol[i] = 0.;
00248         }
00249 
00250         OrientedBox box;
00251         ErrorCode rval = tool->box( root, box );
00252         CHECKERR( mbi, rval );
00253         tot_vol = box.volume();
00254     }
00255 
00256     virtual ErrorCode visit( EntityHandle, int, bool& descend )
00257     {
00258 
00259         descend = true;
00260         return MB_SUCCESS;
00261     }
00262 
00263     virtual ErrorCode leaf( EntityHandle node )
00264     {
00265 
00266         Range tris;
00267         ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
00268         unsigned count = tris.size();
00269 
00270         sum += count;
00271         sqr += ( count * count );
00272         if( min > count ) min = count;
00273         if( max < count ) max = count;
00274 
00275         for( unsigned i = 0; i < ten_buckets_max; ++i )
00276         {
00277             if( count > std::pow( (double)10, (int)( i + 1 ) ) )
00278             {
00279                 ten_buckets[i] += 1;
00280                 OrientedBox box;
00281                 rval = tool->box( node, box );
00282                 CHECKERR( mbi, rval );
00283                 ten_buckets_vol[i] += box.volume();
00284             }
00285         }
00286 
00287         leaves++;
00288 
00289         return rval;
00290     }
00291 
00292     std::string commafy( int num )
00293     {
00294         std::stringstream str;
00295         str << num;
00296         std::string s = str.str();
00297 
00298         int n = s.size();
00299         for( int i = n - 3; i >= 1; i -= 3 )
00300         {
00301             s.insert( i, 1, ',' );
00302             n++;
00303         }
00304 
00305         return s;
00306     }
00307 
00308     void write_results( std::ostream& out )
00309     {
00310 
00311         out << commafy( sum ) << " triangles in " << commafy( leaves ) << " leaves." << std::endl;
00312 
00313         double avg    = sum / (double)leaves;
00314         double stddev = std_dev( sqr, sum, leaves );
00315 
00316         out << "Tris per leaf: Min " << min << ", Max " << max << ", avg " << avg << ", stddev " << stddev << std::endl;
00317 
00318         for( unsigned i = 0; i < ten_buckets_max; ++i )
00319         {
00320             if( ten_buckets[i] )
00321             {
00322                 out << "Leaves exceeding " << std::pow( (double)10, (int)( i + 1 ) )
00323                     << " triangles: " << ten_buckets[i];
00324 
00325                 double frac_total_vol = ten_buckets_vol[i] / tot_vol;
00326                 double avg_ftv        = frac_total_vol / ten_buckets[i];
00327 
00328                 out << " (avg " << avg_ftv * 100.0 << "% of OBB volume)" << std::endl;
00329             }
00330         }
00331     }
00332 };
00333 
00334 /* no longer supported
00335 static std::string make_property_string( DagMC& dag, EntityHandle eh, std::vector<std::string>
00336 &properties )
00337 {
00338   ErrorCode ret;
00339   std::string propstring;
00340   for( std::vector<std::string>::iterator p = properties.begin();
00341     p != properties.end(); ++p )
00342   {
00343     if( dag.has_prop( eh, *p ) ){
00344       std::vector< std::string> vals;
00345       ret = dag.prop_values( eh, *p, vals );
00346       CHECKERR(dag,ret);
00347       propstring += *p;
00348       if( vals.size() == 1 ){
00349         propstring += "=";
00350         propstring += vals[0];
00351       }
00352       else if( vals.size() > 1 ){
00353         // this property has multiple values, list within brackets
00354         propstring += "=[";
00355         for( std::vector<std::string>::iterator i = vals.begin();
00356              i != vals.end(); ++i )
00357         {
00358             propstring += *i;
00359             propstring += ",";
00360         }
00361         // replace the last trailing comma with a close braket
00362         propstring[ propstring.length()-1 ] = ']';
00363       }
00364       propstring += ", ";
00365     }
00366   }
00367   if( propstring.length() ){
00368     propstring.resize( propstring.length() - 2 ); // drop trailing comma
00369   }
00370   return propstring;
00371 }
00372 */
00373 
00374 ErrorCode obbstat_write( GeomTopoTool& gtt,
00375                          std::vector< int >& volumes,
00376                          std::vector< std::string >& properties,
00377                          std::ostream& out )
00378 {
00379 
00380     ErrorCode ret                = MB_SUCCESS;
00381     OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
00382 
00383     // can assume that volume numbers are valid.
00384     for( std::vector< int >::iterator i = volumes.begin(); i != volumes.end(); ++i )
00385     {
00386         EntityHandle vol_root;
00387         EntityHandle vol = gtt.entity_by_id( 3, *i );
00388         CHECKERR( gtt, ret );
00389 
00390         if( vol == 0 )
00391         {
00392             std::cerr << "ERROR: volume " << *i << " has no entity." << std::endl;
00393             continue;
00394         }
00395 
00396         ret = gtt.get_root( vol, vol_root );
00397         CHECKERR( gtt, ret );
00398 
00399         out << "\nVolume " << *i << " " << std::flush;
00400 
00401         if( gtt.is_implicit_complement( vol ) ) out << "(implicit complement) ";
00402         out << std::endl;
00403 
00404         // get all surfaces in volume
00405         Range surfs;
00406         ret = gtt.get_moab_instance()->get_child_meshsets( vol, surfs );
00407         CHECKERR( gtt, ret );
00408 
00409         out << "   with " << surfs.size() << " surfaces" << std::endl;
00410 
00411         TriStats ts( gtt.get_moab_instance(), &obbtool, vol_root );
00412         ret = obbtool.preorder_traverse( vol_root, ts );
00413         CHECKERR( gtt, ret );
00414         ts.write_results( out );
00415 
00416         if( verbose )
00417         {
00418             out << "Surface list: " << std::flush;
00419             for( Range::iterator j = surfs.begin(); j != surfs.end(); ++j )
00420             {
00421                 out << gtt.global_id( *j );
00422                 if( j + 1 != surfs.end() ) out << ",";
00423             }
00424             out << std::endl;
00425             ret = obbtool.stats( vol_root, out );
00426             CHECKERR( gtt, ret );
00427         }
00428 
00429         out << "\n    ------------ " << std::endl;
00430     }
00431 
00432     return ret;
00433 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines