Field.cpp

00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00013 #include <cstring>
00014 #include <stdexcept>
00015 #include <iostream>
00016 #include <sstream>
00017 #include <algorithm>
00018 
00019 #include <stk_util/util/SimpleArrayOps.hpp>
00020 #include <stk_util/util/string_case_compare.hpp>
00021 
00022 #include <stk_mesh/base/Types.hpp>
00023 #include <stk_mesh/base/Field.hpp>
00024 #include <stk_mesh/base/Part.hpp>
00025 #include <stk_mesh/base/MetaData.hpp>
00026 
00027 using std::strlen;
00028 
00029 namespace stk {
00030 namespace mesh {
00031 
00032 //----------------------------------------------------------------------
00033 
00034 const char * field_state_name( FieldState s )
00035 {
00036   static const char * name_list[] = {
00037     "StateNew" ,
00038     "StateOld" ,
00039     "StateNM1" ,
00040     "StateNM2" ,
00041     "StateNM3" ,
00042     "StateNM4" ,
00043     "ERROR" };
00044 
00045   unsigned i = s ;
00046   if ( StateNM4 < i ) { i = MaximumFieldStates ; }
00047   return name_list[i] ;
00048 }
00049 
00050 //----------------------------------------------------------------------
00051 
00052 namespace {
00053 
00054 struct RestrictionLess {
00055   bool operator()( const FieldBase::Restriction & lhs ,
00056                    const FieldBase::Restriction & rhs ) const
00057     { return lhs.key < rhs.key ; }
00058 
00059   bool operator()( const FieldBase::Restriction & lhs ,
00060                    const EntityKey & rhs ) const
00061     { return lhs.key < rhs ; }
00062 };
00063 
00064 std::vector<FieldBase::Restriction>::const_iterator
00065 find( const std::vector<FieldBase::Restriction> & v ,
00066       const EntityKey & key )
00067 {
00068   std::vector<FieldBase::Restriction>::const_iterator
00069     i = std::lower_bound( v.begin() , v.end() , key , RestrictionLess() );
00070 
00071   if ( i != v.end() && i->key != key ) { i = v.end(); }
00072 
00073   return i ;
00074 }
00075 
00076 }
00077 
00078 //----------------------------------------------------------------------
00079 
00080 FieldBase::~Field()
00081 { }
00082 
00083 FieldBase::Field(
00084   MetaData * arg_mesh_meta_data ,
00085   unsigned   arg_ordinal ,
00086   const std::string & arg_name ,
00087   const DataTraits & arg_traits ,
00088   unsigned   arg_number_of_states ,
00089   FieldState arg_this_state )
00090 : m_name( arg_name ),
00091   m_attribute(),
00092   m_data_traits( arg_traits ),
00093   m_mesh_meta_data( arg_mesh_meta_data ),
00094   m_mesh_meta_data_ordinal( arg_ordinal ),
00095   m_num_states( arg_number_of_states ),
00096   m_this_state( arg_this_state ),
00097   m_rank( 0 ),
00098   m_dim_map()
00099 {
00100   FieldBase * const pzero = NULL ;
00101   const shards::ArrayDimTag * const dzero = NULL ;
00102   Copy<MaximumFieldStates>(    m_field_states , pzero );
00103   Copy<MaximumFieldDimension>( m_dim_tags ,     dzero );
00104 }
00105 
00106 const std::vector<FieldBase::Restriction> & FieldBase::restrictions() const
00107 { return m_field_states[0]->m_dim_map ; }
00108 
00109 std::vector<FieldBase::Restriction> & FieldBase::restrictions()
00110 { return m_field_states[0]->m_dim_map ; }
00111 
00112 //----------------------------------------------------------------------
00113 //----------------------------------------------------------------------
00114 
00115 void
00116 print_field_type( std::ostream                      & arg_msg ,
00117                   const DataTraits                  & arg_traits ,
00118                   unsigned                            arg_rank ,
00119                   const shards::ArrayDimTag * const * arg_tags )
00120 {
00121   arg_msg << "Field<" ;
00122   arg_msg << arg_traits.name ;
00123   for ( unsigned i = 0 ; i < arg_rank ; ++i ) {
00124     arg_msg << "," << arg_tags[i]->name();
00125   }
00126   arg_msg << ">" ;
00127 }
00128 
00129 namespace {
00130 
00131 void throw_name_suffix( const char * method ,
00132                         const std::string & name ,
00133                         const char * suffix )
00134 {
00135   std::ostringstream msg ;
00136   msg << method ;
00137   msg << " FAILED: name = \"" ;
00138   msg << name ;
00139   msg << "\" CANNOT HAVE THE RESERVED STATE SUFFIX \"" ;
00140   msg << suffix ;
00141   msg << "\"" ;
00142   throw std::runtime_error( msg.str() );
00143 }
00144 
00145 // Check for compatibility:
00146 // 1) Scalar type must match
00147 // 2) Number of states must match
00148 // 3) Dimension must be different by at most one rank,
00149 //    where the tags match for the smaller rank.
00150 void verify_field_type( const char                        * arg_method ,
00151                         const FieldBase                   & arg_field ,
00152                         const DataTraits                  & arg_traits ,
00153                         unsigned                            arg_rank ,
00154                         const shards::ArrayDimTag * const * arg_dim_tags ,
00155                         unsigned                    arg_num_states )
00156 {
00157 
00158   const bool ok_traits = arg_traits.is_void
00159                       || & arg_traits == & arg_field.data_traits();
00160 
00161   const bool ok_number_states =
00162     ! arg_num_states || arg_num_states == arg_field.number_of_states();
00163 
00164   bool ok_dimension = ! arg_rank || arg_rank     == arg_field.rank() ||
00165                                     arg_rank + 1 == arg_field.rank() ||
00166                                     arg_rank - 1 == arg_field.rank() ;
00167 
00168   const unsigned check_rank = arg_rank < arg_field.rank() ?
00169                               arg_rank : arg_field.rank() ;
00170 
00171   for ( unsigned i = 0 ; i < check_rank && ok_dimension ; ++i ) {
00172     ok_dimension = arg_dim_tags[i] == arg_field.dimension_tags()[i] ;
00173   }
00174 
00175   if ( ! ok_traits || ! ok_number_states || ! ok_dimension ) {
00176 
00177     std::ostringstream msg ;             
00178 
00179     msg << arg_method << " FAILED: Existing field = " ;
00180 
00181     print_field_type( msg , arg_field.data_traits() ,
00182                             arg_field.rank() , 
00183                             arg_field.dimension_tags() ); 
00184  
00185     msg << "[ name = \"" << arg_field.name();
00186     msg << "\" , #states = " << arg_field.number_of_states() << " ]" ;
00187     msg << " Expected field info = " ;
00188                                            
00189     print_field_type( msg , arg_traits , arg_rank , arg_dim_tags );
00190     msg << "[ #states = " << arg_num_states << " ]" ;
00191    
00192     throw std::runtime_error( msg.str() );   
00193   }
00194 }
00195 
00196 }
00197 
00198 FieldBase * get_field(
00199   const char                        * arg_method ,
00200   const std::string                 & arg_name ,
00201   const DataTraits                  & arg_traits ,
00202   unsigned                            arg_rank ,
00203   const shards::ArrayDimTag * const * arg_dim_tags ,
00204   unsigned                            arg_num_states ,
00205   const std::vector<FieldBase*>     & arg_meta_data_fields )
00206 {
00207   FieldBase * f = NULL ;
00208 
00209   for ( std::vector<FieldBase*>::const_iterator
00210         j =  arg_meta_data_fields.begin() ; 
00211         j != arg_meta_data_fields.end() && NULL == f ; ++j ) {
00212     if ( equal_case( (*j)->name() , arg_name ) ) {
00213 
00214       f = *j ;
00215 
00216       verify_field_type( arg_method , *f , arg_traits ,
00217                          arg_rank , arg_dim_tags , arg_num_states );
00218     }
00219   }
00220   return f ;
00221 }
00222 
00223 
00224 //----------------------------------------------------------------------
00225 
00226 FieldBase *
00227 FieldBase::declare_field(
00228   const std::string                 & arg_name ,
00229   const DataTraits                  & arg_traits ,
00230   unsigned                            arg_rank ,
00231   const shards::ArrayDimTag * const * arg_dim_tags ,
00232   unsigned                            arg_num_states ,
00233   MetaData                          * arg_meta_data ,
00234   std::vector<FieldBase*>           & arg_meta_data_fields )
00235 {
00236   static const char method[] = "stk::mesh::FieldBase::declare_field" ;
00237 
00238   static const char reserved_state_suffix[6][8] = {
00239     "_OLD" , "_N" , "_NM1" , "_NM2" , "_NM3" , "_NM4" };
00240 
00241   // Check that the name does not have a reserved suffix
00242 
00243   for ( unsigned i = 0 ; i < 6 ; ++i ) {
00244     const int len_name   = arg_name.size();
00245     const int len_suffix = strlen( reserved_state_suffix[i] );
00246     const int offset     = len_name - len_suffix ;
00247     if ( 0 <= offset ) {
00248       const char * const name_suffix = arg_name.c_str() + offset ;
00249       if ( equal_case( name_suffix , reserved_state_suffix[i] ) ) {
00250         throw_name_suffix( method , arg_name , reserved_state_suffix[i] );
00251       }
00252     }
00253   }
00254 
00255   // Check that the field of this name has not already been declared
00256 
00257   FieldBase * f[ MaximumFieldStates ] ;
00258 
00259   f[0] = get_field( method , arg_name , arg_traits ,
00260                     arg_rank , arg_dim_tags , arg_num_states ,
00261                     arg_meta_data_fields );
00262 
00263   if ( NULL != f[0] ) {
00264     for ( unsigned i = 1 ; i < arg_num_states ; ++i ) {
00265       f[i] = f[0]->m_field_states[i] ; 
00266     }
00267   }
00268   else {
00269     // Field does not exist then create it
00270 
00271     std::string field_names[ MaximumFieldStates ];
00272 
00273     field_names[0] = arg_name ;
00274 
00275     if ( 2 == arg_num_states ) {
00276       field_names[1] = arg_name ;
00277       field_names[1].append( reserved_state_suffix[0] );
00278     }
00279     else {
00280       for ( unsigned i = 1 ; i < arg_num_states ; ++i ) {
00281         field_names[i] = arg_name ;
00282         field_names[i].append( reserved_state_suffix[i] );
00283       }
00284     }
00285 
00286     for ( unsigned i = 0 ; i < arg_num_states ; ++i ) {
00287 
00288       f[i] = new FieldBase( arg_meta_data , arg_meta_data_fields.size() ,
00289                             field_names[i] , arg_traits ,
00290                             arg_num_states , static_cast<FieldState>(i) );
00291 
00292       arg_meta_data_fields.push_back( f[i] );
00293     }
00294 
00295     for ( unsigned i = 0 ; i < arg_num_states ; ++i ) {
00296       for ( unsigned j = 0 ; j < arg_num_states ; ++j ) {
00297         f[i]->m_field_states[j] = f[j] ;
00298       }
00299     }
00300   }
00301 
00302   // Update the rank and dimension tags
00303 
00304   const unsigned old_rank = f[0]->m_rank ;
00305 
00306   if ( old_rank < arg_rank ) {
00307     for ( unsigned j = 0 ; j < arg_num_states ; ++j ) {
00308       f[j]->m_rank = arg_rank ;
00309       for ( unsigned i = old_rank ; i < arg_rank ; ++i ) {
00310         f[j]->m_dim_tags[i] = arg_dim_tags[i] ;
00311       }
00312     }
00313   }
00314 
00315   return f[0] ;
00316 }
00317 
00318 //----------------------------------------------------------------------
00319 //----------------------------------------------------------------------
00320 
00321 namespace {
00322 
00323 void print_restriction( std::ostream & os ,
00324                         unsigned type ,
00325                         const Part & part ,
00326                         unsigned rank ,
00327                         const FieldBase::Restriction::size_type * stride )
00328 {
00329   os << "{ entity_type(" << type << ") part(" << part.name() << ") : " ;
00330   os << stride[0] ;
00331   for ( unsigned i = 1 ; i < rank ; ++i ) {
00332     if ( ! stride[i] ) {
00333       os << " , 0 " ;
00334     }
00335     else if ( stride[i] % stride[i-1] ) {
00336       os << " , " << stride[i] << " / " << stride[i-1] ;
00337     }
00338     else {
00339       os << " , " << stride[i] / stride[i-1] ;
00340     }
00341   }
00342   os << " }" ;
00343 }
00344 
00345 }
00346 
00347 // Setting the dimension for one field sets the dimension
00348 // for the corresponding fields of the FieldState array.
00349 // If subset exists then replace it.
00350 // If exists or superset exists then do nothing.
00351 
00352 void FieldBase::insert_restriction(
00353   const char     * arg_method ,
00354   EntityType         arg_entity_type ,
00355   const Part     & arg_part ,
00356   const unsigned * arg_stride )
00357 {
00358   FieldBase::Restriction tmp ;
00359 
00360   tmp.key = EntityKey( arg_entity_type , arg_part.mesh_meta_data_ordinal() );
00361 
00362   {
00363     unsigned i = 0 ;
00364     if ( m_rank ) {
00365       for ( i = 0 ; i < m_rank ; ++i ) { tmp.stride[i] = arg_stride[i] ; }
00366     }
00367     else { // Scalar field is 0 == m_rank
00368       i = 1 ;
00369       tmp.stride[0] = 1 ;
00370     }
00371     // Remaining dimensions are 1, no change to stride
00372     for ( ; i < MaximumFieldDimension ; ++i ) {
00373       tmp.stride[i] = tmp.stride[i-1] ;
00374     }
00375 
00376     for ( i = 1 ; i < m_rank ; ++i ) {
00377       if ( 0 == tmp.stride[i] || 0 != tmp.stride[i] % tmp.stride[i-1] ) {
00378         std::ostringstream msg ;
00379         msg << arg_method << " FAILED for " << *this ;
00380         msg << " WITH BAD STRIDE " ;
00381         print_restriction( msg, arg_entity_type, arg_part, m_rank, tmp.stride);
00382         throw std::runtime_error( msg.str() );
00383       }
00384     }
00385   }
00386 
00387   {
00388     FieldBase::RestrictionVector & rMap = restrictions();
00389 
00390     FieldBase::RestrictionVector::iterator i = rMap.begin(), j = rMap.end();
00391 
00392     i = std::lower_bound(i,j,tmp,RestrictionLess());
00393 
00394     if ( i == j || i->key != tmp.key ) {
00395       rMap.insert( i , tmp );
00396     }
00397     else if ( Compare<MaximumFieldDimension>::not_equal(i->stride,tmp.stride) ){
00398       std::ostringstream msg ;
00399       msg << arg_method << " FAILED for " << *this << " " ;
00400       print_restriction( msg, arg_entity_type, arg_part, m_rank, i->stride );
00401       msg << " WITH INCOMPATIBLE REDECLARATION " ;
00402       print_restriction( msg, arg_entity_type, arg_part, m_rank, tmp.stride );
00403       throw std::runtime_error( msg.str() );
00404     }
00405   }
00406 }
00407 
00408 void FieldBase::verify_and_clean_restrictions(
00409   const char       * arg_method ,
00410   const PartVector & arg_all_parts )
00411 {
00412   const EntityKey invalid_key ;
00413   RestrictionVector & rMap = restrictions();
00414   RestrictionVector::iterator i , j ;
00415 
00416   for ( i = rMap.begin() ; i != rMap.end() ; ++i ) {
00417     if ( i->key != invalid_key ) {
00418       const unsigned typeI = entity_type( i->key );
00419       const Part   & partI = * arg_all_parts[ entity_id( i->key ) ];
00420       bool  found_superset = false ;
00421 
00422       for ( j = i + 1 ; j != rMap.end() && ! found_superset ; ++j ) {
00423         if ( j->key != invalid_key ) {
00424           const unsigned typeJ = entity_type( j->key );
00425           const Part   & partJ = * arg_all_parts[ entity_id( j->key ) ];
00426 
00427           if ( typeI == typeJ ) {
00428             const bool found_subset = contain( partI.subsets() , partJ );
00429             found_superset = ! found_subset &&
00430                              contain( partI.supersets() , partJ );
00431 
00432             if ( found_subset || found_superset ) {
00433               if ( Compare< MaximumFieldDimension >::not_equal( i->stride ,
00434                                                                 j->stride ) ) {
00435                 std::ostringstream msg ;
00436                 msg << arg_method << "[" ;
00437                 msg << *this ;
00438                 msg << "] FAILED: " ;
00439                 print_restriction( msg, typeI, partI, m_rank, i->stride );
00440                 if ( found_subset ) { msg << " INCOMPATIBLE SUBSET " ; }
00441                 else                { msg << " INCOMPATIBLE SUPERSET " ; }
00442                 print_restriction( msg, typeJ, partJ, m_rank, j->stride );
00443                 throw std::runtime_error( msg.str() );
00444               }
00445             }
00446 
00447             if ( found_subset ) { j->key = invalid_key; }
00448           }
00449         }
00450         if ( found_superset ) { i->key = invalid_key; }
00451       }
00452     }
00453   }
00454 
00455   // Clean out redundant entries:
00456 
00457   for ( j = i = rMap.begin() ; j != rMap.end() ; ++j ) {
00458     if ( j->key != invalid_key ) {
00459       if ( i->key == invalid_key ) {
00460         *i = *j ;
00461       }
00462       ++i ;
00463     }
00464   }
00465 
00466   rMap.erase( i , j );
00467 }
00468 
00469 //----------------------------------------------------------------------
00470 //----------------------------------------------------------------------
00471 // This part or any superset of this part
00472 
00473 const FieldBase::Restriction &
00474 FieldBase::restriction( unsigned etype , const Part & part ) const
00475 {
00476   static const FieldBase::Restriction empty ;
00477 
00478   const std::vector<FieldBase::Restriction> & rMap = restrictions();
00479   const std::vector<FieldBase::Restriction>::const_iterator ie = rMap.end() ;
00480         std::vector<FieldBase::Restriction>::const_iterator i ;
00481 
00482   const PartVector::const_iterator ipe = part.supersets().end();
00483         PartVector::const_iterator ip  = part.supersets().begin() ;
00484 
00485   // Start with this part:
00486   EntityKey key = EntityKey( etype , part.mesh_meta_data_ordinal() );
00487 
00488   while ( ie == ( i = find( rMap , key ) ) && ipe != ip ) {
00489     // Not found try another superset part:
00490     key = EntityKey( etype , (*ip)->mesh_meta_data_ordinal() );
00491     ++ip ;
00492   }
00493 
00494   return ie == i ? empty : *i ;
00495 }
00496 
00497 unsigned FieldBase::max_size( unsigned entity_type ) const
00498 {
00499   unsigned max = 0 ;
00500 
00501   const std::vector<FieldBase::Restriction> & rMap = restrictions();
00502   const std::vector<FieldBase::Restriction>::const_iterator ie= rMap.end();
00503         std::vector<FieldBase::Restriction>::const_iterator i = rMap.begin();
00504 
00505   for ( ; i != ie ; ++i ) {
00506     if ( i->type() == entity_type ) {
00507       const unsigned len = m_rank ? i->stride[ m_rank - 1 ] : 1 ;
00508       if ( max < len ) { max = len ; }
00509     }
00510   }
00511 
00512   return max ;
00513 }
00514 
00515 //----------------------------------------------------------------------
00516 
00517 std::ostream & operator << ( std::ostream & s , const FieldBase & field )
00518 {
00519   print_field_type( s , field.data_traits() ,
00520                         field.rank() , field.dimension_tags() );
00521   s << "[ name = \"" ;
00522   s << field.name() ;
00523   s << "\" , #states = " ;
00524   s << field.number_of_states();
00525   s << " ]" ;
00526   return s ;
00527 }
00528 
00529 std::ostream & print( std::ostream & s ,
00530                       const char * const b ,
00531                       const FieldBase & field )
00532 {
00533   const PartVector & all_parts = field.mesh_meta_data().get_parts();
00534   const std::vector<FieldBase::Restriction> & rMap = field.restrictions();
00535   s << field ;
00536   s << " {" ;
00537   for ( std::vector<FieldBase::Restriction>::const_iterator
00538         i = rMap.begin() ; i != rMap.end() ; ++i ) {
00539     s << std::endl << b << "  " ;
00540     print_restriction( s, entity_type( i->key ),
00541                        * all_parts[ entity_id( i->key ) ], 
00542                        field.rank(), i->stride);
00543   }
00544   s << std::endl << b << "}" ;
00545   return s ;
00546 }
00547 
00548 //----------------------------------------------------------------------
00549 
00550 } // namespace mesh
00551 } // namespace stk
00552 

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