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_selector_restrictions(),
00068   m_initial_value(NULL),
00069   m_initial_value_num_bytes(0)
00070 {
00071   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::FieldBaseImpl", LOG_FIELD, m_ordinal);
00072 
00073   FieldBase * const pzero = NULL ;
00074   const shards::ArrayDimTag * const dzero = NULL ;
00075   Copy<MaximumFieldStates>(    m_field_states , pzero );
00076   Copy<MaximumFieldDimension>( m_dim_tags ,     dzero );
00077 
00078   for ( unsigned i = 0 ; i < arg_rank ; ++i ) {
00079     m_dim_tags[i] = arg_dim_tags[i];
00080   }
00081 }
00082 
00083 //----------------------------------------------------------------------
00084 FieldBaseImpl::~FieldBaseImpl()
00085 {
00086   if (state() == StateNone) {
00087     void*& init_val = m_initial_value;
00088 
00089     delete [] reinterpret_cast<char*>(init_val);
00090     init_val = NULL;
00091   }
00092 }
00093 
00094 //----------------------------------------------------------------------
00095 const FieldRestrictionVector & FieldBaseImpl::restrictions() const
00096 { return m_field_states[0]->m_impl.m_dim_map ; }
00097 
00098 const FieldRestrictionVector & FieldBaseImpl::selector_restrictions() const
00099 { return m_field_states[0]->m_impl.m_selector_restrictions ; }
00100 
00101 FieldRestrictionVector & FieldBaseImpl::restrictions()
00102 { return m_field_states[0]->m_impl.m_dim_map ; }
00103 
00104 FieldRestrictionVector & FieldBaseImpl::selector_restrictions()
00105 { return m_field_states[0]->m_impl.m_selector_restrictions ; }
00106 
00107 
00108 //----------------------------------------------------------------------
00109 
00110 // Setting the dimension for one field sets the dimension
00111 // for the corresponding fields of the FieldState array.
00112 // If subset exists then replace it.
00113 // If exists or superset exists then do nothing.
00114 
00115 void FieldBaseImpl::insert_restriction(
00116   const char     * arg_method ,
00117   EntityRank       arg_entity_rank ,
00118   const Part     & arg_part ,
00119   const unsigned * arg_stride,
00120   const void*      arg_init_value )
00121 {
00122   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::insert_restriction", LOG_FIELD, m_ordinal);
00123 
00124   FieldRestriction tmp( arg_entity_rank , arg_part.mesh_meta_data_ordinal() );
00125 
00126   {
00127     unsigned i = 0 ;
00128     if ( m_field_rank ) {
00129       for ( i = 0 ; i < m_field_rank ; ++i ) { tmp.stride(i) = arg_stride[i] ; }
00130     }
00131     else { // Scalar field is 0 == m_field_rank
00132       i = 1 ;
00133       tmp.stride(0) = 1 ;
00134     }
00135     // Remaining dimensions are 1, no change to stride
00136     for ( ; i < MaximumFieldDimension ; ++i ) {
00137       tmp.stride(i) = tmp.stride(i-1) ;
00138     }
00139 
00140     for ( i = 1 ; i < m_field_rank ; ++i ) {
00141       const bool bad_stride = 0 == tmp.stride(i) ||
00142                               0 != tmp.stride(i) % tmp.stride(i-1);
00143       ThrowErrorMsgIf( bad_stride,
00144           arg_method << " FAILED for " << *this <<
00145           " WITH BAD STRIDE " <<
00146           print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));;
00147     }
00148   }
00149 
00150   if (arg_init_value != NULL) {
00151     //insert_restriction can be called multiple times for the same field, giving
00152     //the field different lengths on different mesh-parts.
00153     //We will only store one initial-value array, we need to store the one with
00154     //maximum length for this field so that it can be used to initialize data
00155     //for all field-restrictions. For the parts on which the field is shorter,
00156     //a subset of the initial-value array will be used.
00157     //
00158     //We want to end up storing the longest arg_init_value array for this field.
00159     //
00160     //Thus, we call set_initial_value only if the current length is longer
00161     //than what's already been stored.
00162 
00163     //length in bytes is num-scalars X sizeof-scalar:
00164 
00165     size_t num_scalars = 1;
00166     //if rank > 0, then field is not a scalar field, so num-scalars is
00167     //obtained from the stride array:
00168     if (m_field_rank > 0) num_scalars = tmp.stride(m_field_rank-1);
00169 
00170     size_t sizeof_scalar = m_data_traits.size_of;
00171     size_t nbytes = sizeof_scalar * num_scalars;
00172 
00173     size_t old_nbytes = 0;
00174     if (get_initial_value() != NULL) {
00175       old_nbytes = get_initial_value_num_bytes();
00176     }
00177 
00178     if (nbytes > old_nbytes) {
00179       set_initial_value(arg_init_value, num_scalars, nbytes);
00180     }
00181   }
00182 
00183   {
00184     FieldRestrictionVector & restrs = restrictions();
00185 
00186     FieldRestrictionVector::iterator restr = restrs.begin();
00187     FieldRestrictionVector::iterator last_restriction = restrs.end();
00188 
00189     restr = std::lower_bound(restr,last_restriction,tmp);
00190 
00191     const bool new_restriction = ( ( restr == last_restriction ) || !(*restr == tmp) );
00192 
00193     if ( new_restriction ) {
00194       // New field restriction, verify we are not committed:
00195       ThrowRequireMsg(!m_meta_data->is_commit(), "mesh MetaData has been committed.");
00196       unsigned num_subsets = 0;
00197       for(FieldRestrictionVector::iterator i=restrs.begin(), iend=restrs.end(); i!=iend; ++i) {
00198         if (i->entity_rank() != arg_entity_rank) continue;
00199 
00200         const Part& partI = *m_meta_data->get_parts()[i->part_ordinal()];
00201         bool found_subset = contain(arg_part.subsets(), partI);
00202         if (found_subset) {
00203           ThrowErrorMsgIf( i->not_equal_stride(tmp),
00204             arg_method << " FAILED for " << *this << " " <<
00205             print_restriction( *i, arg_entity_rank, arg_part, m_field_rank ) <<
00206             " WITH INCOMPATIBLE REDECLARATION " <<
00207             print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
00208           *i = tmp;
00209           ++num_subsets;
00210         }
00211 
00212         bool found_superset = contain(arg_part.supersets(), partI);
00213         if (found_superset) {
00214           ThrowErrorMsgIf( i->not_equal_stride(tmp),
00215             arg_method << " FAILED for " << *this << " " <<
00216             print_restriction( *i, arg_entity_rank, arg_part, m_field_rank ) <<
00217             " WITH INCOMPATIBLE REDECLARATION " <<
00218             print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
00219           //if there's already a restriction for a superset of this part, then 
00220           //there's nothing to do and we're out of here..
00221           return;
00222         }
00223       }
00224       if (num_subsets == 0) {
00225         restrs.insert( restr , tmp );
00226       }
00227       else {
00228         //if subsets were found, we replaced them with the new restriction. so now we need
00229         //to sort and unique the vector, and trim it to remove any duplicates:
00230         std::sort(restrs.begin(), restrs.end());
00231         FieldRestrictionVector::iterator it = std::unique(restrs.begin(), restrs.end());
00232         restrs.resize(it - restrs.begin());
00233       }
00234     }
00235     else {
00236       ThrowErrorMsgIf( restr->not_equal_stride(tmp),
00237           arg_method << " FAILED for " << *this << " " <<
00238           print_restriction( *restr, arg_entity_rank, arg_part, m_field_rank ) <<
00239           " WITH INCOMPATIBLE REDECLARATION " <<
00240           print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
00241     }
00242   }
00243 }
00244 
00245 void FieldBaseImpl::insert_restriction(
00246   const char     * arg_method ,
00247   EntityRank       arg_entity_rank ,
00248   const Selector & arg_selector ,
00249   const unsigned * arg_stride,
00250   const void*      arg_init_value )
00251 {
00252   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::insert_restriction", LOG_FIELD, m_ordinal);
00253 
00254   FieldRestriction tmp( arg_entity_rank , arg_selector );
00255 
00256   {
00257     unsigned i = 0 ;
00258     if ( m_field_rank ) {
00259       for ( i = 0 ; i < m_field_rank ; ++i ) { tmp.stride(i) = arg_stride[i] ; }
00260     }
00261     else { // Scalar field is 0 == m_field_rank
00262       i = 1 ;
00263       tmp.stride(0) = 1 ;
00264     }
00265     // Remaining dimensions are 1, no change to stride
00266     for ( ; i < MaximumFieldDimension ; ++i ) {
00267       tmp.stride(i) = tmp.stride(i-1) ;
00268     }
00269 
00270     for ( i = 1 ; i < m_field_rank ; ++i ) {
00271       const bool bad_stride = 0 == tmp.stride(i) ||
00272                               0 != tmp.stride(i) % tmp.stride(i-1);
00273       ThrowErrorMsgIf( bad_stride,
00274           arg_method << " FAILED for " << *this <<
00275           " WITH BAD STRIDE!");
00276     }
00277   }
00278 
00279   if (arg_init_value != NULL) {
00280     //insert_restriction can be called multiple times for the same field, giving
00281     //the field different lengths on different mesh-parts.
00282     //We will only store one initial-value array, we need to store the one with
00283     //maximum length for this field so that it can be used to initialize data
00284     //for all field-restrictions. For the parts on which the field is shorter,
00285     //a subset of the initial-value array will be used.
00286     //
00287     //We want to end up storing the longest arg_init_value array for this field.
00288     //
00289     //Thus, we call set_initial_value only if the current length is longer
00290     //than what's already been stored.
00291 
00292     //length in bytes is num-scalars X sizeof-scalar:
00293 
00294     size_t num_scalars = 1;
00295     //if rank > 0, then field is not a scalar field, so num-scalars is
00296     //obtained from the stride array:
00297     if (m_field_rank > 0) num_scalars = tmp.stride(m_field_rank-1);
00298 
00299     size_t sizeof_scalar = m_data_traits.size_of;
00300     size_t nbytes = sizeof_scalar * num_scalars;
00301 
00302     size_t old_nbytes = 0;
00303     if (get_initial_value() != NULL) {
00304       old_nbytes = get_initial_value_num_bytes();
00305     }
00306 
00307     if (nbytes > old_nbytes) {
00308       set_initial_value(arg_init_value, num_scalars, nbytes);
00309     }
00310   }
00311 
00312   {
00313     FieldRestrictionVector & srvec = selector_restrictions();
00314 
00315     bool restriction_already_exists = false;
00316     for(FieldRestrictionVector::const_iterator it=srvec.begin(), it_end=srvec.end();
00317         it!=it_end; ++it) {
00318       if (tmp == *it) {
00319         restriction_already_exists = true;
00320         if (tmp.not_equal_stride(*it)) {
00321           ThrowErrorMsg("Incompatible selector field-restrictions!");
00322         }
00323       }
00324     }
00325 
00326     if ( !restriction_already_exists ) {
00327       // New field restriction, verify we are not committed:
00328       ThrowRequireMsg(!m_meta_data->is_commit(), "mesh MetaData has been committed.");
00329       srvec.push_back( tmp );
00330     }
00331   }
00332 }
00333 
00334 void FieldBaseImpl::verify_and_clean_restrictions(
00335   const char       * arg_method ,
00336   const Part& superset,
00337   const Part& subset,
00338   const PartVector & arg_all_parts )
00339 {
00340   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::verify_and_clean_restrictions", LOG_FIELD, m_ordinal);
00341 
00342   FieldRestrictionVector & restrs = restrictions();
00343 
00344   //Check whether both 'superset' and 'subset' are in this field's restrictions.
00345   //If they are, make sure they are compatible and remove the subset restriction.
00346   FieldRestrictionVector::iterator superset_restriction = restrs.end();
00347   FieldRestrictionVector::iterator subset_restriction = restrs.end();
00348   for (FieldRestrictionVector::iterator i = restrs.begin() ; i != restrs.end() ; ++i ) {
00349     if (i->part_ordinal() == superset.mesh_meta_data_ordinal()) {
00350       superset_restriction = i;
00351       if (subset_restriction != restrs.end() && subset_restriction->entity_rank() == superset_restriction->entity_rank()) break;
00352     }
00353     if (i->part_ordinal() == subset.mesh_meta_data_ordinal()) {
00354       subset_restriction = i;
00355       if (superset_restriction != restrs.end() && subset_restriction->entity_rank() == superset_restriction->entity_rank()) break;
00356     }
00357   }
00358 
00359   if (superset_restriction != restrs.end() && subset_restriction != restrs.end() &&
00360       superset_restriction->entity_rank() == subset_restriction->entity_rank()) {
00361     ThrowErrorMsgIf( superset_restriction->not_equal_stride(*subset_restriction),
00362       "Incompatible field restrictions for parts "<<superset.name()<<" and "<<subset.name());
00363 
00364     restrs.erase(subset_restriction);
00365   }
00366 }
00367 
00368 const void* FieldBaseImpl::get_initial_value() const
00369 {
00370   return m_field_states[0]->m_impl.m_initial_value;
00371 }
00372 
00373 void* FieldBaseImpl::get_initial_value() {
00374   return m_field_states[0]->m_impl.m_initial_value;
00375 }
00376 
00377 unsigned FieldBaseImpl::get_initial_value_num_bytes() const {
00378   return m_field_states[0]->m_impl.m_initial_value_num_bytes;
00379 }
00380 
00381 void FieldBaseImpl::set_initial_value(const void* new_initial_value, unsigned num_scalars, unsigned num_bytes) {
00382   void*& init_val = m_field_states[0]->m_impl.m_initial_value;
00383 
00384   delete [] reinterpret_cast<char*>(init_val);
00385   init_val = new char[num_bytes];
00386 
00387   m_field_states[0]->m_impl.m_initial_value_num_bytes = num_bytes;
00388 
00389   m_data_traits.copy(init_val, new_initial_value, num_scalars);
00390 }
00391 
00392 
00393 //----------------------------------------------------------------------
00394 //----------------------------------------------------------------------
00395 // This part or any superset of this part
00396 
00397 const FieldRestriction &
00398 FieldBaseImpl::restriction( unsigned entity_rank , const Part & part ) const
00399 {
00400   static const FieldRestriction empty ;
00401 
00402   const FieldRestrictionVector & rMap = restrictions();
00403   const FieldRestrictionVector::const_iterator ie = rMap.end() ;
00404         FieldRestrictionVector::const_iterator i ;
00405 
00406   const PartVector::const_iterator ipe = part.supersets().end();
00407         PartVector::const_iterator ip  = part.supersets().begin() ;
00408 
00409   // Start with this part:
00410   //(putting static here helps performance significantly but is NOT THREAD SAFE !!!)
00411   static FieldRestriction restr;
00412   restr.set_entity_rank( entity_rank );
00413   restr.set_part_ordinal( part.mesh_meta_data_ordinal() );
00414 
00415   while ( ie == ( i = find( rMap , restr ) ) && ipe != ip ) {
00416     // Not found try another superset part:
00417     restr.set_entity_rank( entity_rank );
00418     restr.set_part_ordinal( (*ip)->mesh_meta_data_ordinal() );
00419     ++ip ;
00420   }
00421 
00422   return ie == i ? empty : *i ;
00423 }
00424 
00425 unsigned FieldBaseImpl::max_size( unsigned entity_rank ) const
00426 {
00427   unsigned max = 0 ;
00428 
00429   const FieldRestrictionVector & rMap = restrictions();
00430   const FieldRestrictionVector::const_iterator ie = rMap.end() ;
00431         FieldRestrictionVector::const_iterator i = rMap.begin();
00432 
00433   for ( ; i != ie ; ++i ) {
00434     if ( i->entity_rank() == entity_rank ) {
00435       const unsigned len = m_field_rank ? i->stride( m_field_rank - 1 ) : 1 ;
00436       if ( max < len ) { max = len ; }
00437     }
00438   }
00439 
00440   return max ;
00441 }
00442 
00443 void FieldBaseImpl::set_field_states( FieldBase ** field_states)
00444 {
00445   TraceIfWatching("stk::mesh::impl::FieldBaseImpl::set_field_states", LOG_FIELD, m_ordinal);
00446 
00447   for (unsigned i = 0; i < m_num_states; ++i) {
00448     m_field_states[i] = field_states[i];
00449   }
00450 }
00451 
00452 //----------------------------------------------------------------------
00453 
00454 //----------------------------------------------------------------------
00455 
00456 std::ostream & operator << ( std::ostream & s , const FieldBaseImpl & field )
00457 {
00458   s << "FieldBaseImpl<" ;
00459   s << field.data_traits().name ;
00460   for ( unsigned i = 0 ; i < field.rank() ; ++i ) {
00461     s << "," << field.dimension_tags()[i]->name();
00462   }
00463   s << ">" ;
00464 
00465   s << "[ name = \"" ;
00466   s << field.name() ;
00467   s << "\" , #states = " ;
00468   s << field.number_of_states();
00469   s << " ]" ;
00470   return s ;
00471 }
00472 
00473 std::ostream & print( std::ostream & s ,
00474                       const char * const b ,
00475                       const FieldBase & field )
00476 {
00477   const PartVector & all_parts = MetaData::get(field).get_parts();
00478   const std::vector<FieldBase::Restriction> & rMap = field.restrictions();
00479   s << field.name() ;
00480   s << " {" ;
00481   for ( FieldBase::RestrictionVector::const_iterator
00482         i = rMap.begin() ; i != rMap.end() ; ++i ) {
00483     s << std::endl << b << "  " ;
00484     i->print( s, i->entity_rank(), * all_parts[ i->part_ordinal() ], field.rank() );
00485     s << std::endl;
00486   }
00487   s << std::endl << b << "}" ;
00488   return s ;
00489 }
00490 
00491 //----------------------------------------------------------------------
00492 
00493 
00494 
00495 
00496 } // namespace impl
00497 } // namespace mesh
00498 } // namespace stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines