![]() |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include <iostream>
#include <fstream>
#include <sstream>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/OrientedBoxTreeTool.hpp"
#include "moab/CartVect.hpp"
#include "moab/OrientedBox.hpp"
#include "moab/GeomTopoTool.hpp"
#include "cgm2moab.hpp"
Go to the source code of this file.
Classes | |
class | TriCounter |
class | TriStats |
Functions | |
ErrorCode | obbvis_create (GeomTopoTool >t, std::vector< int > &volumes, int grid, std::string &filename) |
static double | std_dev (double sqr, double sum, double count) |
ErrorCode | obbstat_write (GeomTopoTool >t, std::vector< int > &volumes, std::vector< std::string > &properties, std::ostream &out) |
ErrorCode obbstat_write | ( | GeomTopoTool & | gtt, |
std::vector< int > & | volumes, | ||
std::vector< std::string > & | properties, | ||
std::ostream & | out | ||
) |
Definition at line 374 of file obb_analysis.cpp.
References moab::Range::begin(), CHECKERR, moab::Range::end(), moab::GeomTopoTool::entity_by_id(), ErrorCode, moab::Interface::get_child_meshsets(), moab::GeomTopoTool::get_moab_instance(), moab::GeomTopoTool::get_root(), moab::GeomTopoTool::global_id(), moab::GeomTopoTool::is_implicit_complement(), MB_SUCCESS, moab::GeomTopoTool::obb_tree(), moab::OrientedBoxTreeTool::preorder_traverse(), moab::Range::size(), moab::OrientedBoxTreeTool::stats(), and verbose.
{
ErrorCode ret = MB_SUCCESS;
OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
// can assume that volume numbers are valid.
for( std::vector< int >::iterator i = volumes.begin(); i != volumes.end(); ++i )
{
EntityHandle vol_root;
EntityHandle vol = gtt.entity_by_id( 3, *i );
CHECKERR( gtt, ret );
if( vol == 0 )
{
std::cerr << "ERROR: volume " << *i << " has no entity." << std::endl;
continue;
}
ret = gtt.get_root( vol, vol_root );
CHECKERR( gtt, ret );
out << "\nVolume " << *i << " " << std::flush;
if( gtt.is_implicit_complement( vol ) ) out << "(implicit complement) ";
out << std::endl;
// get all surfaces in volume
Range surfs;
ret = gtt.get_moab_instance()->get_child_meshsets( vol, surfs );
CHECKERR( gtt, ret );
out << " with " << surfs.size() << " surfaces" << std::endl;
TriStats ts( gtt.get_moab_instance(), &obbtool, vol_root );
ret = obbtool.preorder_traverse( vol_root, ts );
CHECKERR( gtt, ret );
ts.write_results( out );
if( verbose )
{
out << "Surface list: " << std::flush;
for( Range::iterator j = surfs.begin(); j != surfs.end(); ++j )
{
out << gtt.global_id( *j );
if( j + 1 != surfs.end() ) out << ",";
}
out << std::endl;
ret = obbtool.stats( vol_root, out );
CHECKERR( gtt, ret );
}
out << "\n ------------ " << std::endl;
}
return ret;
}
ErrorCode obbvis_create | ( | GeomTopoTool & | gtt, |
std::vector< int > & | volumes, | ||
int | grid, | ||
std::string & | filename | ||
) |
Definition at line 50 of file obb_analysis.cpp.
References moab::Range::begin(), center(), CHECKERR, moab::Core::create_element(), moab::Core::create_vertices(), moab::Range::end(), moab::GeomTopoTool::entity_by_id(), ErrorCode, moab::Interface::get_child_meshsets(), moab::GeomTopoTool::get_moab_instance(), moab::GeomTopoTool::get_obb(), moab::GeomTopoTool::get_root(), MB_TAG_BYTES, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_EXCL, MB_TYPE_INTEGER, MBHEX, moab::Range::merge(), moab::GeomTopoTool::obb_tree(), moab::OrientedBoxTreeTool::preorder_traverse(), moab::Range::size(), moab::sum(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), verbose, and moab::Core::write_file().
{
OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
CartVect min, max;
EntityHandle vol = gtt.entity_by_id( 3, volumes.front() );
double middle[3];
double axis1[3], axis2[3], axis3[3];
double minPt[3], maxPt[3];
ErrorCode rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
// compute min and max verticies
for( int i = 0; i < 3; i++ )
{
double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
minPt[i] = middle[i] - sum;
maxPt[i] = middle[i] + sum;
}
CHECKERR( gtt, rval );
/* Compute an axis-aligned bounding box of all the requested volumes */
for( std::vector< int >::iterator i = volumes.begin() + 1; i != volumes.end(); ++i )
{
CartVect i_min, i_max;
vol = gtt.entity_by_id( 3, *i );
rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
// compute min and max verticies
for( int j = 0; j < 3; j++ )
{
double sum = fabs( axis1[j] ) + fabs( axis2[j] ) + fabs( axis3[j] );
minPt[j] = middle[j] - sum;
maxPt[j] = middle[j] + sum;
}
for( int j = 0; j < 3; ++j )
{
min[j] = std::min( min[j], i_min[j] );
max[j] = std::max( max[j], i_max[j] );
}
}
// These vectors could be repurposed to describe an OBB without changing the loops below
CartVect center( middle );
CartVect v1( axis1 );
CartVect v2( axis2 );
CartVect v3( axis3 );
/* Compute the vertices of the visualization grid. Calculation points are at the center
of each cell in this grid, so make grid+1 vertices in each direction. */
int numpoints = pow( (double)( grid + 1 ), 3 );
double* pgrid = new double[numpoints * 3];
int idx = 0;
for( int i = 0; i < numpoints * 3; ++i )
pgrid[i] = 0.0;
for( int i = 0; i <= grid; ++i )
{
CartVect x = -v1 + ( ( v1 * 2.0 ) * ( i / (double)grid ) );
for( int j = 0; j <= grid; ++j )
{
CartVect y = -v2 + ( ( v2 * 2.0 ) * ( j / (double)grid ) );
for( int k = 0; k <= grid; ++k )
{
CartVect z = -v3 + ( ( v3 * 2.0 ) * ( k / (double)grid ) );
CartVect p = center + x + y + z;
for( int d = 0; d < 3; ++d )
{
pgrid[idx++] = p[d];
}
}
}
}
/* Create a new MOAB to use for output, and build the vertex grid */
Core mb2;
Range r;
rval = mb2.create_vertices( pgrid, numpoints, r );
CHECKERR( mb2, rval );
Tag lttag;
rval = mb2.tag_get_handle( "LEAFTRIS", sizeof( int ), MB_TYPE_INTEGER, lttag,
MB_TAG_EXCL | MB_TAG_CREAT | MB_TAG_BYTES | MB_TAG_DENSE, 0 );
CHECKERR( mb2, rval );
int row = grid + 1;
int side = row * row;
EntityHandle connect[8];
EntityHandle hex;
// offset from grid corner to grid center
CartVect grid_hex_center_offset = ( v1 + v2 + v3 ) * 2 * ( 1.0 / grid );
// collect all the surfaces from the requested volumes to iterate over --
// this prevents checking a shared surface more than once.
Range surfs;
for( std::vector< int >::iterator it = volumes.begin(); it != volumes.end(); ++it )
{
vol = gtt.entity_by_id( 3, *it );
Range it_surfs;
rval = gtt.get_moab_instance()->get_child_meshsets( vol, it_surfs );
CHECKERR( gtt, rval );
surfs.merge( it_surfs );
}
std::cout << "visualizing " << surfs.size() << " surfaces." << std::endl;
/* Build hexes for any point with 1 or more leaftris */
for( int i = 0; i < grid; ++i )
{
for( int j = 0; j < grid; ++j )
{
for( int k = 0; k < grid; ++k )
{
idx = ( side * k ) + ( row * j ) + i;
assert( idx + 1 + row + side > numpoints - 1 );
CartVect loc = CartVect( ( pgrid + ( idx * 3 ) ) ) + grid_hex_center_offset;
TriCounter tc( gtt.get_moab_instance(), &obbtool, loc );
for( Range::iterator it = surfs.begin(); it != surfs.end(); ++it )
{
EntityHandle surf_tree;
rval = gtt.get_root( *it, surf_tree );
CHECKERR( gtt, rval );
rval = obbtool.preorder_traverse( surf_tree, tc );
CHECKERR( gtt, rval );
}
if( tc.count == 0 ) continue;
connect[0] = r[idx];
connect[1] = r[idx + 1];
connect[2] = r[idx + 1 + row];
connect[3] = r[idx + row];
connect[4] = r[idx + side];
connect[5] = r[idx + 1 + side];
connect[6] = r[idx + 1 + row + side];
connect[7] = r[idx + row + side];
rval = mb2.create_element( MBHEX, connect, 8, hex );
CHECKERR( mb2, rval );
rval = mb2.tag_set_data( lttag, &hex, 1, &( tc.count ) );
CHECKERR( mb2, rval );
}
}
}
if( verbose )
{
std::cout << "Writing " << filename << std::endl;
}
rval = mb2.write_file( filename.c_str() );
CHECKERR( mb2, rval );
return rval;
}
static double std_dev | ( | double | sqr, |
double | sum, | ||
double | count | ||
) | [inline, static] |
Definition at line 215 of file obb_analysis.cpp.
{
sum /= count;
sqr /= count;
return sqrt( sqr - sum * sum );
}