BoundBox.cpp

Go to the documentation of this file.

    1 #include "moab/Interface.hpp"
    2 #include "moab/BoundBox.hpp"
    3  
    4 // this is used for spherical elements, bounding box needs to be updated due to curvature
    5 #include "moab/IntxMesh/IntxUtils.hpp"
    6  
    7 namespace moab
    8 {
    9 ErrorCode BoundBox::update( Interface& iface, const Range& elems, bool spherical, double radius )
   10 {
   11     ErrorCode rval;
   12     bMin = CartVect( HUGE_VAL );
   13     bMax = CartVect( -HUGE_VAL );
   14  
   15     CartVect coords;
   16     EntityHandle const *conn = NULL, *conn2 = NULL;
   17     int len = 0, len2 = 0;
   18     Range::const_iterator i;
   19     CartVect coordverts[27];  // maximum nodes per element supported by MOAB
   20  
   21     // vertices
   22     const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
   23     for( i = elems.begin(); i != elem_begin; ++i )
   24     {
   25         rval = iface.get_coords( &*i, 1, coords.array() );
   26         if( MB_SUCCESS != rval ) return rval;
   27         update_min( coords.array() );
   28         update_max( coords.array() );
   29     }
   30  
   31     // elements with vertex-handle connectivity list
   32     const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
   33     std::vector< EntityHandle > dum_vector;
   34     for( i = elem_begin; i != poly_begin; ++i )
   35     {
   36         rval = iface.get_connectivity( *i, conn, len, true, &dum_vector );
   37         if( MB_SUCCESS != rval ) return rval;
   38  
   39         rval = iface.get_coords( conn, len, &( coordverts[0][0] ) );
   40         if( MB_SUCCESS != rval ) return rval;
   41         for( int j = 0; j < len; ++j )
   42         {
   43             update_min( coordverts[j].array() );
   44             update_max( coordverts[j].array() );
   45         }
   46         // if spherical, use gnomonic projection to improve the computed box
   47         // it is used now for coupling only
   48         if( spherical ) update_box_spherical_elem( coordverts, len, radius );
   49     }
   50  
   51     // polyhedra
   52     const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
   53     for( i = poly_begin; i != set_begin; ++i )
   54     {
   55         rval = iface.get_connectivity( *i, conn, len, true );
   56         if( MB_SUCCESS != rval ) return rval;
   57  
   58         for( int j = 0; j < len; ++j )
   59         {
   60             rval = iface.get_connectivity( conn[j], conn2, len2 );
   61             if( MB_SUCCESS != rval ) return rval;
   62             rval = iface.get_coords( conn2, len2, &( coordverts[0][0] ) );
   63             if( MB_SUCCESS != rval ) return rval;
   64             for( int k = 0; k < len2; ++k )
   65             {
   66                 update_min( coordverts[k].array() );
   67                 update_max( coordverts[k].array() );
   68             }
   69         }
   70     }
   71  
   72     // sets
   73     BoundBox box;
   74     for( i = set_begin; i != elems.end(); ++i )
   75     {
   76         Range tmp_elems;
   77         rval = iface.get_entities_by_handle( *i, tmp_elems );
   78         if( MB_SUCCESS != rval ) return rval;
   79         rval = box.update( iface, tmp_elems, spherical, radius );
   80         if( MB_SUCCESS != rval ) return rval;
   81  
   82         update( box );
   83     }
   84  
   85     return MB_SUCCESS;
   86 }
   87 void BoundBox::update_box_spherical_elem( const CartVect* verts, int len, double R )
   88 {
   89     // decide first the gnomonic plane, based on first coordinate
   90     int plane = -1;
   91     CartVect pos;
   92     IntxUtils::decide_gnomonic_plane( verts[0], plane );
   93     double in_plane_positions[20];  // at most polygons with 10 edges; do we need to revise this?
   94     for( int i = 0; i < len && i < 10; i++ )
   95         IntxUtils::gnomonic_projection( verts[i], R, plane, in_plane_positions[2 * i], in_plane_positions[2 * i + 1] );
   96     // look for points on the intersection between edges in gnomonic plane and coordinate axes
   97     // if yes, reverse projection of intx point, and update box
   98     double oriented_area2 = 0;
   99     // if oriented area is != zero, it means the origin is inside the polygon, so we need to upgrade
  100     // the box with the origin, so we do not miss anything (the top of the dome)
  101     for( int i = 0; i < len; i++ )
  102     {
  103         int i1 = ( i + 1 ) % len;  // next vertex in polygon
  104         // check intx with x axis
  105         double ax = in_plane_positions[2 * i], ay = in_plane_positions[2 * i + 1];
  106         double bx = in_plane_positions[2 * i1], by = in_plane_positions[2 * i1 + 1];
  107         if( ay * by < 0 )  // it intersects with x axis
  108         {
  109             double alfa  = ay / ( ay - by );
  110             double xintx = ax + alfa * ( bx - ax );
  111             IntxUtils::reverse_gnomonic_projection( xintx, 0, R, plane, pos );
  112             update_min( pos.array() );
  113             update_max( pos.array() );
  114         }
  115         if( ax * bx < 0 )  // it intersects with y axis
  116         {
  117             double alfa  = ax / ( ax - bx );
  118             double yintx = ay + alfa * ( by - ay );
  119             IntxUtils::reverse_gnomonic_projection( 0, yintx, R, plane, pos );
  120             update_min( pos.array() );
  121             update_max( pos.array() );
  122         }
  123         oriented_area2 += ( ax * by - ay * bx );
  124     }
  125     if( fabs( oriented_area2 ) > R * R * 1.e-6 )  // origin is in the interior, add the center
  126     {
  127         IntxUtils::reverse_gnomonic_projection( 0, 0, R, plane, pos );
  128         update_min( pos.array() );
  129         update_max( pos.array() );
  130     }
  131 }
  132 }  // namespace moab