Sierra Toolkit Version of the Day
FieldBaseImpl.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 
00009 #include <stk_mesh/baseImpl/FieldBaseImpl.hpp>
00010 
00011 #include <cstring>
00012 #include <stdexcept>
00013 #include <iostream>
00014 #include <sstream>
00015 #include <algorithm>
00016 
00017 #include <stk_util/util/SimpleArrayOps.hpp>
00018 
00019 #include <stk_mesh/base/Types.hpp>
00020 #include <stk_mesh/base/FieldBase.hpp>
00021 #include <stk_mesh/base/Part.hpp>
00022 #include <stk_mesh/base/MetaData.hpp>
00023 #include <stk_mesh/base/Trace.hpp>
00024 
00025 namespace stk {
00026 namespace mesh {
00027 namespace impl {
00028 
00029 //----------------------------------------------------------------------
00030 
00031 namespace {
00032 
00033 FieldRestrictionVector::const_iterator
00034   find( const FieldRestrictionVector & v , const FieldRestriction & restr )
00035 {
00036   FieldRestrictionVector::const_iterator
00037     i = std::lower_bound( v.begin() , v.end() , restr );
00038 
00039   if ( i != v.end() && !(*i == restr) ) { i = v.end(); }
00040 
00041   return i ;
00042 }
00043 
00044 }
00045 
00046 //----------------------------------------------------------------------
00047 
00048 FieldBaseImpl::FieldBaseImpl(
00049     MetaData                   * arg_mesh_meta_data ,
00050     unsigned                     arg_ordinal ,
00051     const std::string          & arg_name ,
00052     const DataTraits           & arg_traits ,
00053     unsigned                     arg_rank,
00054     const shards::ArrayDimTag  * const * arg_dim_tags,
00055     unsigned                     arg_number_of_states ,
00056     FieldState                   arg_this_state
00057     )
00058 : m_name( arg_name ),
00059   m_attribute(),
00060   m_data_traits( arg_traits ),
00061   m_meta_data( arg_mesh_meta_data ),
00062   m_ordinal( arg_ordinal ),
00063   m_num_states( arg_number_of_states ),
00064   m_this_state( arg_this_state ),
00065   m_field_rank( arg_rank ),
00066   m_dim_map(),
00067   m_initial_value(NULL)
00068 {
00069   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::FieldBaseImpl", LOG_FIELD, m_ordinal);
00070 
00071   FieldBase * const pzero = NULL ;
00072   const shards::ArrayDimTag * const dzero = NULL ;
00073   Copy<MaximumFieldStates>(    m_field_states , pzero );
00074   Copy<MaximumFieldDimension>( m_dim_tags ,     dzero );
00075 
00076   for ( unsigned i = 0 ; i < arg_rank ; ++i ) {
00077     m_dim_tags[i] = arg_dim_tags[i];
00078   }
00079 }
00080 
00081 //----------------------------------------------------------------------
00082 FieldBaseImpl::~FieldBaseImpl()
00083 {
00084   if (state() == StateNone) {
00085     void*& init_val = m_initial_value;
00086 
00087     delete [] reinterpret_cast<char*>(init_val);
00088     init_val = NULL;
00089   }
00090 }
00091 
00092 //----------------------------------------------------------------------
00093 const FieldRestrictionVector & FieldBaseImpl::restrictions() const
00094 { return m_field_states[0]->m_impl.m_dim_map ; }
00095 
00096 FieldRestrictionVector & FieldBaseImpl::restrictions()
00097 { return m_field_states[0]->m_impl.m_dim_map ; }
00098 
00099 
00100 //----------------------------------------------------------------------
00101 
00102 // Setting the dimension for one field sets the dimension
00103 // for the corresponding fields of the FieldState array.
00104 // If subset exists then replace it.
00105 // If exists or superset exists then do nothing.
00106 
00107 void FieldBaseImpl::insert_restriction(
00108   const char     * arg_method ,
00109   EntityRank       arg_entity_rank ,
00110   const Part     & arg_part ,
00111   const unsigned * arg_stride,
00112   const void*      arg_init_value )
00113 {
00114   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::insert_restriction", LOG_FIELD, m_ordinal);
00115 
00116   FieldRestriction tmp( arg_entity_rank , arg_part.mesh_meta_data_ordinal() );
00117 
00118   {
00119     unsigned i = 0 ;
00120     if ( m_field_rank ) {
00121       for ( i = 0 ; i < m_field_rank ; ++i ) { tmp.stride(i) = arg_stride[i] ; }
00122     }
00123     else { // Scalar field is 0 == m_field_rank
00124       i = 1 ;
00125       tmp.stride(0) = 1 ;
00126     }
00127     // Remaining dimensions are 1, no change to stride
00128     for ( ; i < MaximumFieldDimension ; ++i ) {
00129       tmp.stride(i) = tmp.stride(i-1) ;
00130     }
00131 
00132     for ( i = 1 ; i < m_field_rank ; ++i ) {
00133       const bool bad_stride = 0 == tmp.stride(i) ||
00134                               0 != tmp.stride(i) % tmp.stride(i-1);
00135       ThrowErrorMsgIf( bad_stride,
00136           arg_method << " FAILED for " << *this <<
00137           " WITH BAD STRIDE " <<
00138           print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));;
00139     }
00140   }
00141 
00142   if (arg_init_value != NULL) {
00143     //insert_restriction can be called multiple times for the same field, giving
00144     //the field different lengths on different mesh-parts.
00145     //We will only store one initial-value array, we need to store the one with
00146     //maximum length for this field so that it can be used to initialize data
00147     //for all field-restrictions. For the parts on which the field is shorter,
00148     //a subset of the initial-value array will be used.
00149     //
00150     //We want to end up storing the longest arg_init_value array for this field.
00151     //
00152     //Thus, we call set_initial_value only if the current length is longer
00153     //than what's already been stored.
00154 
00155     //length in bytes is num-scalars X sizeof-scalar:
00156 
00157     size_t num_scalars = 1;
00158     //if rank > 0, then field is not a scalar field, so num-scalars is
00159     //obtained from the stride array:
00160     if (m_field_rank > 0) num_scalars = tmp.stride(m_field_rank-1);
00161 
00162     size_t sizeof_scalar = m_data_traits.size_of;
00163     size_t nbytes = sizeof_scalar * num_scalars;
00164 
00165     size_t old_nbytes = 0;
00166     if (get_initial_value() != NULL) {
00167       old_nbytes = get_initial_value_num_bytes();
00168     }
00169 
00170     if (nbytes > old_nbytes) {
00171       set_initial_value(arg_init_value, num_scalars, nbytes);
00172     }
00173   }
00174 
00175   {
00176     FieldRestrictionVector & rMap = restrictions();
00177 
00178     FieldRestrictionVector::iterator restr = rMap.begin();
00179     FieldRestrictionVector::iterator last_restriction = rMap.end();
00180 
00181     restr = std::lower_bound(restr,last_restriction,tmp);
00182 
00183     const bool new_restriction = ( ( restr == last_restriction ) || !(*restr == tmp) );
00184 
00185     if ( new_restriction ) {
00186       // New field restriction, verify we are not committed:
00187       ThrowRequireMsg(!m_meta_data->is_commit(), "mesh MetaData has been committed.");
00188       rMap.insert( restr , tmp );
00189     }
00190     else {
00191       ThrowErrorMsgIf( restr->not_equal_stride(tmp),
00192           arg_method << " FAILED for " << *this << " " <<
00193           print_restriction( *restr, arg_entity_rank, arg_part, m_field_rank ) <<
00194           " WITH INCOMPATIBLE REDECLARATION " <<
00195           print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
00196     }
00197   }
00198 }
00199 
00200 void FieldBaseImpl::verify_and_clean_restrictions(
00201   const char       * arg_method ,
00202   const PartVector & arg_all_parts )
00203 {
00204   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::verify_and_clean_restrictions", LOG_FIELD, m_ordinal);
00205 
00206   const FieldRestriction invalid_restr ;
00207   FieldRestrictionVector & rMap = restrictions();
00208 
00209   // Search for redundant field restrictions. A restriction R is redundant if there exists
00210   // another restriction on a superset of R's part that is for entities of the same rank as R.
00211   for (FieldRestrictionVector::iterator i = rMap.begin() ; i != rMap.end() ; ++i ) {
00212     if ( *i != invalid_restr ) {
00213       const EntityRank rankI = i->entity_rank();
00214       const Part     & partI = * arg_all_parts[ i->part_ordinal() ];
00215       bool    found_superset = false ;
00216 
00217       for (FieldRestrictionVector::iterator j = i + 1 ; j != rMap.end() && ! found_superset ; ++j ) {
00218         if ( *j != invalid_restr ) {
00219           const EntityRank rankJ = j->entity_rank();
00220           const Part     & partJ = * arg_all_parts[ j->part_ordinal() ];
00221 
00222           if ( rankI == rankJ ) {
00223             const bool found_subset = contain( partI.subsets() , partJ );
00224             found_superset = ! found_subset &&
00225                              contain( partI.supersets() , partJ );
00226 
00227             if ( found_subset || found_superset ) {
00228               ThrowErrorMsgIf( i->not_equal_stride(*j),
00229                   arg_method << "[" << *this << "] FAILED: " <<
00230                   print_restriction( *i, rankI, partI, m_field_rank ) <<
00231                   ( found_subset ? " INCOMPATIBLE SUBSET " : " INCOMPATIBLE SUPERSET ") <<
00232                   print_restriction( *j, rankJ, partJ, m_field_rank ));
00233             }
00234 
00235             if ( found_subset ) { *j = invalid_restr; }
00236           }
00237         }
00238         if ( found_superset ) { *i = invalid_restr; }
00239       }
00240     }
00241   }
00242 
00243   // Clean out redundant entries:
00244   FieldRestrictionVector::iterator new_end = std::remove(rMap.begin(), rMap.end(), invalid_restr);
00245   rMap.erase(new_end, rMap.end());
00246 }
00247 
00248 const void* FieldBaseImpl::get_initial_value() const
00249 {
00250   return m_field_states[0]->m_impl.m_initial_value;
00251 }
00252 
00253 void* FieldBaseImpl::get_initial_value() {
00254   return m_field_states[0]->m_impl.m_initial_value;
00255 }
00256 
00257 unsigned FieldBaseImpl::get_initial_value_num_bytes() const {
00258   return m_field_states[0]->m_impl.m_initial_value_num_bytes;
00259 }
00260 
00261 void FieldBaseImpl::set_initial_value(const void* new_initial_value, unsigned num_scalars, unsigned num_bytes) {
00262   void*& init_val = m_field_states[0]->m_impl.m_initial_value;
00263 
00264   delete [] reinterpret_cast<char*>(init_val);
00265   init_val = new char[num_bytes];
00266 
00267   m_field_states[0]->m_impl.m_initial_value_num_bytes = num_bytes;
00268 
00269   m_data_traits.copy(init_val, new_initial_value, num_scalars);
00270 }
00271 
00272 
00273 //----------------------------------------------------------------------
00274 //----------------------------------------------------------------------
00275 // This part or any superset of this part
00276 
00277 const FieldRestriction &
00278 FieldBaseImpl::restriction( unsigned entity_rank , const Part & part ) const
00279 {
00280   static const FieldRestriction empty ;
00281 
00282   const FieldRestrictionVector & rMap = restrictions();
00283   const FieldRestrictionVector::const_iterator ie = rMap.end() ;
00284         FieldRestrictionVector::const_iterator i ;
00285 
00286   const PartVector::const_iterator ipe = part.supersets().end();
00287         PartVector::const_iterator ip  = part.supersets().begin() ;
00288 
00289   // Start with this part:
00290   FieldRestriction restr( entity_rank , part.mesh_meta_data_ordinal() );
00291 
00292   while ( ie == ( i = find( rMap , restr ) ) && ipe != ip ) {
00293     // Not found try another superset part:
00294     restr = FieldRestriction( entity_rank , (*ip)->mesh_meta_data_ordinal() );
00295     ++ip ;
00296   }
00297 
00298   return ie == i ? empty : *i ;
00299 }
00300 
00301 unsigned FieldBaseImpl::max_size( unsigned entity_rank ) const
00302 {
00303   unsigned max = 0 ;
00304 
00305   const FieldRestrictionVector & rMap = restrictions();
00306   const FieldRestrictionVector::const_iterator ie = rMap.end() ;
00307         FieldRestrictionVector::const_iterator i = rMap.begin();
00308 
00309   for ( ; i != ie ; ++i ) {
00310     if ( i->entity_rank() == entity_rank ) {
00311       const unsigned len = m_field_rank ? i->stride( m_field_rank - 1 ) : 1 ;
00312       if ( max < len ) { max = len ; }
00313     }
00314   }
00315 
00316   return max ;
00317 }
00318 
00319 void FieldBaseImpl::set_field_states( FieldBase ** field_states)
00320 {
00321   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::set_field_states", LOG_FIELD, m_ordinal);
00322 
00323   for (unsigned i = 0; i < m_num_states; ++i) {
00324     m_field_states[i] = field_states[i];
00325   }
00326 }
00327 
00328 //----------------------------------------------------------------------
00329 
00330 //----------------------------------------------------------------------
00331 
00332 std::ostream & operator << ( std::ostream & s , const FieldBaseImpl & field )
00333 {
00334   s << "FieldBaseImpl<" ;
00335   s << field.data_traits().name ;
00336   for ( unsigned i = 0 ; i < field.rank() ; ++i ) {
00337     s << "," << field.dimension_tags()[i]->name();
00338   }
00339   s << ">" ;
00340 
00341   s << "[ name = \"" ;
00342   s << field.name() ;
00343   s << "\" , #states = " ;
00344   s << field.number_of_states();
00345   s << " ]" ;
00346   return s ;
00347 }
00348 
00349 std::ostream & print( std::ostream & s ,
00350                       const char * const b ,
00351                       const FieldBase & field )
00352 {
00353   const PartVector & all_parts = MetaData::get(field).get_parts();
00354   const std::vector<FieldBase::Restriction> & rMap = field.restrictions();
00355   s << field.name() ;
00356   s << " {" ;
00357   for ( FieldBase::RestrictionVector::const_iterator
00358         i = rMap.begin() ; i != rMap.end() ; ++i ) {
00359     s << std::endl << b << "  " ;
00360     i->print( s, i->entity_rank(), * all_parts[ i->part_ordinal() ], field.rank() );
00361     s << std::endl;
00362   }
00363   s << std::endl << b << "}" ;
00364   return s ;
00365 }
00366 
00367 //----------------------------------------------------------------------
00368 
00369 
00370 
00371 
00372 } // namespace impl
00373 } // namespace mesh
00374 } // namespace stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends