OctTree.cpp

Go to the documentation of this file.
00001 /*------------------------------------------------------------------------*/
00002 /*      stk : Parallel Heterogneous Dynamic unstructured Mesh         */
00003 /*                Copyright (2007) Sandia Corporation                     */
00004 /*                                                                        */
00005 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00006 /*  license for use of this work by or on behalf of the U.S. Government.  */
00007 /*                                                                        */
00008 /*  This library is free software; you can redistribute it and/or modify  */
00009 /*  it under the terms of the GNU Lesser General Public License as        */
00010 /*  published by the Free Software Foundation; either version 2.1 of the  */
00011 /*  License, or (at your option) any later version.                       */
00012 /*                                                                        */
00013 /*  This library is distributed in the hope that it will be useful,       */
00014 /*  but WITHOUT ANY WARRANTY; without even the implied warranty of        */
00015 /*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU     */
00016 /*  Lesser General Public License for more details.                       */
00017 /*                                                                        */
00018 /*  You should have received a copy of the GNU Lesser General Public      */
00019 /*  License along with this library; if not, write to the Free Software   */
00020 /*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307   */
00021 /*  USA                                                                   */
00022 /*------------------------------------------------------------------------*/
00029 #include <sstream>
00030 #include <ostream>
00031 #include <stdexcept>
00032 #include <stk_search/OctTree.hpp>
00033 
00034 namespace std {
00035 
00036 ostream & operator << ( ostream & os , const stk::OctTreeKey & otk )
00037 {
00038   unsigned j = 0 ;
00039 
00040   while ( j < stk::OctTreeKey::MaxDepth ) {
00041     os << otk.index(++j);
00042   }
00043   
00044   return os ;
00045 }
00046 
00047 }
00048 
00049 namespace stk {
00050 
00051 enum { OK = StaticAssert< OctTreeKey::NWord == 2 >::OK };
00052 
00053 namespace {
00054   void throw_index( const unsigned d , const unsigned i )
00055   {
00056     std::ostringstream msg ;
00057     msg << "OctTree.cpp: index[" << d << "] = " << i << " is out of range [1..8]" ;
00058     throw std::range_error( msg.str() );
00059   }
00060   
00061   void throw_depth( const unsigned min_depth, const unsigned depth )
00062   {
00063     const unsigned m = stk::OctTreeKey::MaxDepth ;
00064     std::ostringstream msg ;
00065     msg << "OctTree.cpp: depth = " << depth << " is out of range [" << min_depth << ".." << m << "]" ;
00066     throw std::range_error( msg.str() );
00067   }
00068 }
00069   
00070 unsigned OctTreeKey::depth() const
00071 {
00072   const unsigned which = m_value[1] ? 1 : 0 ;
00073   const unsigned val   = m_value[which] ;
00074 
00075   int d = IndexPerWord ;
00076   while ( d-- && ( val & ( MaskIndex << ( BitsPerIndex * d ) ) ) );
00077   return ( which + 1 ) * IndexPerWord - ( d + 1 );
00078 }
00079 
00080 unsigned OctTreeKey::index( const unsigned Depth ) const
00081 {
00082   if ( Depth < 1 || MaxDepth < Depth ) { throw_depth( 1, Depth ); }
00083 
00084   const unsigned which = ( Depth - 1 ) / IndexPerWord ;
00085   const unsigned shift = BitsPerWord -
00086                          BitsPerIndex * ( Depth % IndexPerWord ) ;
00087 
00088   return ( m_value[ which ] >> shift ) & MaskIndex ;
00089 }
00090 
00091 OctTreeKey & OctTreeKey::clear_index( const unsigned Depth )
00092 {
00093   if ( Depth < 1 || MaxDepth < Depth ) { throw_depth( 1, Depth ); }
00094 
00095   const value_type m = MaskIndex ;
00096   const unsigned which = ( Depth - 1 ) / IndexPerWord ;
00097   const unsigned shift = BitsPerWord -
00098                          BitsPerIndex * ( Depth % IndexPerWord ) ;
00099   const unsigned mask = ~( m << shift );
00100 
00101   m_value[ which ] &= mask ;
00102 
00103   return *this ;
00104 }
00105 
00106 OctTreeKey &
00107 OctTreeKey::set_index( const unsigned Depth , const unsigned Index )
00108 {
00109   if ( Depth < 1 || MaxDepth < Depth ) { throw_depth( 1, Depth ); }
00110   if ( Index < 1 || 8        < Index ) { throw_index( Depth , Index ); }
00111 
00112   const value_type m = MaskIndex ;
00113   const unsigned which = ( Depth - 1 ) / IndexPerWord ;
00114   const unsigned shift = BitsPerWord -
00115                          BitsPerIndex * ( Depth % IndexPerWord ) ;
00116 
00117   ( m_value[which] &= ~( m << shift ) ) |= Index << shift ;
00118 
00119   return *this ;
00120 }
00121 
00122 OctTreeKey &
00123 OctTreeKey::set_value( const unsigned * val )
00124 {
00125   Copy<NWord>( m_value , 0u );
00126 
00127   for ( unsigned d = 1 ; d <= MaxDepth ; ++d ) {
00128     const unsigned which = ( d - 1 ) / IndexPerWord ;
00129     const unsigned shift = BitsPerWord - BitsPerIndex * ( d % IndexPerWord ) ;
00130     const unsigned index = ( val[ which ] >> shift ) & MaskIndex ;
00131 
00132     if ( 8 < index ) { throw_index( d , index ); }
00133 
00134     m_value[ which ] |= index << shift ;
00135   }
00136   return *this ;
00137 }
00138 
00139 bool OctTreeKey::intersect( const OctTreeKey & k ) const
00140 {
00141   const unsigned which = m_value[0] != k.m_value[0] ? 0 : (
00142                          m_value[1] != k.m_value[1] ? 1 : 2 );
00143 
00144   bool result = which == 2 ;
00145 
00146   if ( ! result ) {
00147     const value_type lval =   m_value[which];
00148     const value_type rval = k.m_value[which];
00149 
00150     result = true ;
00151     for ( unsigned d = IndexPerWord ; result && d ; ) {
00152       --d ;
00153       const value_type mask = MaskIndex << ( BitsPerIndex * d );
00154       const value_type l = lval & mask ;
00155       const value_type r = rval & mask ;
00156       if ( l && r ) { result = l == r ; }
00157       else          { d = 0 ; }
00158     }
00159   }
00160   return result ;
00161 }
00162 
00163 //----------------------------------------------------------------------
00164 //----------------------------------------------------------------------
00165 
00166 namespace {
00167 
00168 const unsigned tree_size[11] = {
00169   OctTreeSize< 0>::value ,
00170   OctTreeSize< 1>::value ,
00171   OctTreeSize< 2>::value ,
00172   OctTreeSize< 3>::value ,
00173   OctTreeSize< 4>::value ,
00174   OctTreeSize< 5>::value ,
00175   OctTreeSize< 6>::value ,
00176   OctTreeSize< 7>::value ,
00177   OctTreeSize< 8>::value ,
00178   OctTreeSize< 9>::value ,
00179   OctTreeSize<10>::value
00180 };
00181 
00182 }
00183 
00184 unsigned oct_tree_size( const unsigned Depth )
00185 {
00186   if ( stk::OctTreeKey::MaxDepth < Depth )
00187     { throw_depth( 0, Depth ); }
00188 
00189   return tree_size[ Depth ];
00190 }
00191 
00192 unsigned oct_tree_offset( const unsigned Depth , const OctTreeKey & k )
00193 {
00194   if ( stk::OctTreeKey::MaxDepth < Depth )
00195     { throw_depth( 0, Depth ); }
00196 
00197   unsigned index = 0 ;
00198 
00199   unsigned d = Depth ;
00200 
00201   while ( d ) {
00202     --d ;
00203     const unsigned i = k.index( Depth - d );
00204     if ( i ) { index += 1 + ( i - 1 ) * tree_size[ d ] ; }
00205     else     { d = 0 ; } // Stop at a zero index
00206   }
00207 
00208   return index ;
00209 }
00210 
00211 //----------------------------------------------------------------------
00212 //----------------------------------------------------------------------
00213 
00214 OctTreeKey hsfc3d( const unsigned Depth , const unsigned * const coord )
00215 {
00216   // Gray & Inverse Gray coding for 3D octree
00217 
00218 /*
00219   const unsigned gray_fwd[8] = { 0 , 1 , 3 , 2 , 6 , 7 , 5 , 4 };
00220 */
00221   const unsigned gray_inv[8] = { 0 , 1 , 3 , 2 , 7 , 6 , 4 , 5 };
00222 
00223   // Three dimensional HSFC rotation data
00224 
00225   const unsigned hsfc3d_rotation_perm[8][3] =
00226     { /* 0 */ { 2 , 1 , 0 } ,
00227       /* 1 */ { 0 , 2 , 1 } ,
00228       /* 2 */ { 0 , 1 , 2 } ,
00229       /* 3 */ { 2 , 0 , 1 } ,
00230       /* 4 */ { 2 , 0 , 1 } ,
00231       /* 5 */ { 0 , 1 , 2 } ,
00232       /* 6 */ { 0 , 2 , 1 } ,
00233       /* 7 */ { 2 , 1 , 0 } };
00234 
00235   const unsigned hsfc3d_rotation_flip[8][3] =
00236       { /* 0 */ { 0 , 0 , 0 } ,
00237         /* 1 */ { 0 , 0 , 0 } ,
00238         /* 2 */ { 0 , 0 , 0 } ,
00239         /* 3 */ { 1 , 1 , 0 } ,
00240         /* 4 */ { 0 , 1 , 1 } ,
00241         /* 5 */ { 0 , 0 , 0 } ,
00242         /* 6 */ { 0 , 1 , 1 } ,
00243         /* 7 */ { 1 , 0 , 1 } };
00244 
00245   OctTreeKey key ;
00246 
00247   // Initial rotation
00248 
00249   unsigned axis[3] ;
00250 
00251   axis[0] = 0 << 1 ;
00252   axis[1] = 1 << 1 ;
00253   axis[2] = 2 << 1 ;
00254 
00255   for ( unsigned d = 1 ; d <= Depth ; ++d ) {
00256 
00257     const unsigned s = OctTreeKey::BitsPerWord - d ;
00258 
00259     // The ordinal for this depth
00260 
00261     const unsigned ord = gray_inv[
00262       (((( coord[ axis[0] >> 1 ] >> s ) ^ axis[0] ) & 01 ) << 0 ) |
00263       (((( coord[ axis[1] >> 1 ] >> s ) ^ axis[1] ) & 01 ) << 1 ) |
00264       (((( coord[ axis[2] >> 1 ] >> s ) ^ axis[2] ) & 01 ) << 2 ) ];
00265 
00266     key.set_index( d , ord + 1 );
00267 
00268     // Determine the recursive rotation for the next ordinal
00269     {
00270       const unsigned * const p = hsfc3d_rotation_perm[ ord ] ;
00271       const unsigned * const f = hsfc3d_rotation_flip[ ord ] ;
00272 
00273       unsigned tmp[3] ;
00274 
00275       tmp[0] = axis[0] ;
00276       tmp[1] = axis[1] ;
00277       tmp[2] = axis[2] ;
00278 
00279       axis[0] = tmp[ p[0] ] ^ f[0] ;
00280       axis[1] = tmp[ p[1] ] ^ f[1] ;
00281       axis[2] = tmp[ p[2] ] ^ f[2] ;
00282     }
00283   }
00284 
00285   return key ;
00286 }
00287 
00288 
00289 }
00290 
00291 

Generated on Tue Jul 13 09:27:32 2010 for Sierra Toolkit by  doxygen 1.4.7