Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
BoundBox.cpp
Go to the documentation of this file.
00001 #include "moab/Interface.hpp"
00002 #include "moab/BoundBox.hpp"
00003 
00004 // this is used for spherical elements, bounding box needs to be updated due to curvature
00005 #include "moab/IntxMesh/IntxUtils.hpp"
00006 
00007 namespace moab
00008 {
00009 ErrorCode BoundBox::update( Interface& iface, const Range& elems, bool spherical, double radius )
00010 {
00011     ErrorCode rval;
00012     bMin = CartVect( HUGE_VAL );
00013     bMax = CartVect( -HUGE_VAL );
00014 
00015     CartVect coords;
00016     EntityHandle const *conn = NULL, *conn2 = NULL;
00017     int len = 0, len2 = 0;
00018     Range::const_iterator i;
00019     CartVect coordverts[27];  // maximum nodes per element supported by MOAB
00020 
00021     // vertices
00022     const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
00023     for( i = elems.begin(); i != elem_begin; ++i )
00024     {
00025         rval = iface.get_coords( &*i, 1, coords.array() );
00026         if( MB_SUCCESS != rval ) return rval;
00027         update_min( coords.array() );
00028         update_max( coords.array() );
00029     }
00030 
00031     // elements with vertex-handle connectivity list
00032     const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
00033     std::vector< EntityHandle > dum_vector;
00034     for( i = elem_begin; i != poly_begin; ++i )
00035     {
00036         rval = iface.get_connectivity( *i, conn, len, true, &dum_vector );
00037         if( MB_SUCCESS != rval ) return rval;
00038 
00039         rval = iface.get_coords( conn, len, &( coordverts[0][0] ) );
00040         if( MB_SUCCESS != rval ) return rval;
00041         for( int j = 0; j < len; ++j )
00042         {
00043             update_min( coordverts[j].array() );
00044             update_max( coordverts[j].array() );
00045         }
00046         // if spherical, use gnomonic projection to improve the computed box
00047         // it is used now for coupling only
00048         if( spherical ) update_box_spherical_elem( coordverts, len, radius );
00049     }
00050 
00051     // polyhedra
00052     const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
00053     for( i = poly_begin; i != set_begin; ++i )
00054     {
00055         rval = iface.get_connectivity( *i, conn, len, true );
00056         if( MB_SUCCESS != rval ) return rval;
00057 
00058         for( int j = 0; j < len; ++j )
00059         {
00060             rval = iface.get_connectivity( conn[j], conn2, len2 );
00061             if( MB_SUCCESS != rval ) return rval;
00062             rval = iface.get_coords( conn2, len2, &( coordverts[0][0] ) );
00063             if( MB_SUCCESS != rval ) return rval;
00064             for( int k = 0; k < len2; ++k )
00065             {
00066                 update_min( coordverts[k].array() );
00067                 update_max( coordverts[k].array() );
00068             }
00069         }
00070     }
00071 
00072     // sets
00073     BoundBox box;
00074     for( i = set_begin; i != elems.end(); ++i )
00075     {
00076         Range tmp_elems;
00077         rval = iface.get_entities_by_handle( *i, tmp_elems );
00078         if( MB_SUCCESS != rval ) return rval;
00079         rval = box.update( iface, tmp_elems, spherical, radius );
00080         if( MB_SUCCESS != rval ) return rval;
00081 
00082         update( box );
00083     }
00084 
00085     return MB_SUCCESS;
00086 }
00087 void BoundBox::update_box_spherical_elem( const CartVect* verts, int len, double R )
00088 {
00089     // decide first the gnomonic plane, based on first coordinate
00090     int plane = -1;
00091     CartVect pos;
00092     IntxUtils::decide_gnomonic_plane( verts[0], plane );
00093     double in_plane_positions[20];  // at most polygons with 10 edges; do we need to revise this?
00094     for( int i = 0; i < len && i < 10; i++ )
00095         IntxUtils::gnomonic_projection( verts[i], R, plane, in_plane_positions[2 * i], in_plane_positions[2 * i + 1] );
00096     // look for points on the intersection between edges in gnomonic plane and coordinate axes
00097     // if yes, reverse projection of intx point, and update box
00098     double oriented_area2 = 0;
00099     // if oriented area is != zero, it means the origin is inside the polygon, so we need to upgrade
00100     // the box with the origin, so we do not miss anything (the top of the dome)
00101     for( int i = 0; i < len; i++ )
00102     {
00103         int i1 = ( i + 1 ) % len;  // next vertex in polygon
00104         // check intx with x axis
00105         double ax = in_plane_positions[2 * i], ay = in_plane_positions[2 * i + 1];
00106         double bx = in_plane_positions[2 * i1], by = in_plane_positions[2 * i1 + 1];
00107         if( ay * by < 0 )  // it intersects with x axis
00108         {
00109             double alfa  = ay / ( ay - by );
00110             double xintx = ax + alfa * ( bx - ax );
00111             IntxUtils::reverse_gnomonic_projection( xintx, 0, R, plane, pos );
00112             update_min( pos.array() );
00113             update_max( pos.array() );
00114         }
00115         if( ax * bx < 0 )  // it intersects with y axis
00116         {
00117             double alfa  = ax / ( ax - bx );
00118             double yintx = ay + alfa * ( by - ay );
00119             IntxUtils::reverse_gnomonic_projection( 0, yintx, R, plane, pos );
00120             update_min( pos.array() );
00121             update_max( pos.array() );
00122         }
00123         oriented_area2 += ( ax * by - ay * bx );
00124     }
00125     if( fabs( oriented_area2 ) > R * R * 1.e-6 )  // origin is in the interior, add the center
00126     {
00127         IntxUtils::reverse_gnomonic_projection( 0, 0, R, plane, pos );
00128         update_min( pos.array() );
00129         update_max( pos.array() );
00130     }
00131 }
00132 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines