Sierra Toolkit  Version of the Day
FieldBaseImpl.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 #include <stk_mesh/baseImpl/FieldBaseImpl.hpp>
10 
11 #include <cstring>
12 #include <stdexcept>
13 #include <iostream>
14 #include <sstream>
15 #include <algorithm>
16 
17 #include <stk_util/util/SimpleArrayOps.hpp>
18 
19 #include <stk_mesh/base/Types.hpp>
20 #include <stk_mesh/base/FieldBase.hpp>
21 #include <stk_mesh/base/Part.hpp>
22 #include <stk_mesh/base/MetaData.hpp>
23 #include <stk_mesh/base/Trace.hpp>
24 
25 namespace stk_classic {
26 namespace mesh {
27 namespace impl {
28 
29 //----------------------------------------------------------------------
30 
31 namespace {
32 
33 FieldRestrictionVector::const_iterator
34  find( const FieldRestrictionVector & v , const FieldRestriction & restr )
35 {
36  FieldRestrictionVector::const_iterator
37  i = std::lower_bound( v.begin() , v.end() , restr );
38 
39  if ( i != v.end() && !(*i == restr) ) { i = v.end(); }
40 
41  return i ;
42 }
43 
44 }
45 
46 //----------------------------------------------------------------------
47 
48 FieldBaseImpl::FieldBaseImpl(
49  MetaData * arg_mesh_meta_data ,
50  unsigned arg_ordinal ,
51  const std::string & arg_name ,
52  const DataTraits & arg_traits ,
53  unsigned arg_rank,
54  const shards::ArrayDimTag * const * arg_dim_tags,
55  unsigned arg_number_of_states ,
56  FieldState arg_this_state
57  )
58 : m_name( arg_name ),
59  m_attribute(),
60  m_data_traits( arg_traits ),
61  m_meta_data( arg_mesh_meta_data ),
62  m_ordinal( arg_ordinal ),
63  m_num_states( arg_number_of_states ),
64  m_this_state( arg_this_state ),
65  m_field_rank( arg_rank ),
66  m_dim_map(),
67  m_selector_restrictions(),
68  m_initial_value(NULL),
69  m_initial_value_num_bytes(0)
70 {
71  TraceIfWatching("stk_classic::mesh::impl::FieldBaseImpl::FieldBaseImpl", LOG_FIELD, m_ordinal);
72 
73  FieldBase * const pzero = NULL ;
74  const shards::ArrayDimTag * const dzero = NULL ;
75  Copy<MaximumFieldStates>( m_field_states , pzero );
76  Copy<MaximumFieldDimension>( m_dim_tags , dzero );
77 
78  for ( unsigned i = 0 ; i < arg_rank ; ++i ) {
79  m_dim_tags[i] = arg_dim_tags[i];
80  }
81 }
82 
83 //----------------------------------------------------------------------
84 FieldBaseImpl::~FieldBaseImpl()
85 {
86  if (state() == StateNone) {
87  void*& init_val = m_initial_value;
88 
89  delete [] reinterpret_cast<char*>(init_val);
90  init_val = NULL;
91  }
92 }
93 
94 //----------------------------------------------------------------------
95 const FieldRestrictionVector & FieldBaseImpl::restrictions() const
96 { return m_field_states[0]->m_impl.m_dim_map ; }
97 
98 const FieldRestrictionVector & FieldBaseImpl::selector_restrictions() const
99 { return m_field_states[0]->m_impl.m_selector_restrictions ; }
100 
101 FieldRestrictionVector & FieldBaseImpl::restrictions()
102 { return m_field_states[0]->m_impl.m_dim_map ; }
103 
104 FieldRestrictionVector & FieldBaseImpl::selector_restrictions()
105 { return m_field_states[0]->m_impl.m_selector_restrictions ; }
106 
107 
108 //----------------------------------------------------------------------
109 
110 // Setting the dimension for one field sets the dimension
111 // for the corresponding fields of the FieldState array.
112 // If subset exists then replace it.
113 // If exists or superset exists then do nothing.
114 
115 void FieldBaseImpl::insert_restriction(
116  const char * arg_method ,
117  EntityRank arg_entity_rank ,
118  const Part & arg_part ,
119  const unsigned * arg_stride,
120  const void* arg_init_value )
121 {
122  TraceIfWatching("stk_classic::mesh::impl::FieldBaseImpl::insert_restriction", LOG_FIELD, m_ordinal);
123 
124  FieldRestriction tmp( arg_entity_rank , arg_part.mesh_meta_data_ordinal() );
125 
126  {
127  unsigned i = 0 ;
128  if ( m_field_rank ) {
129  for ( i = 0 ; i < m_field_rank ; ++i ) { tmp.stride(i) = arg_stride[i] ; }
130  }
131  else { // Scalar field is 0 == m_field_rank
132  i = 1 ;
133  tmp.stride(0) = 1 ;
134  }
135  // Remaining dimensions are 1, no change to stride
136  for ( ; i < MaximumFieldDimension ; ++i ) {
137  tmp.stride(i) = tmp.stride(i-1) ;
138  }
139 
140  for ( i = 1 ; i < m_field_rank ; ++i ) {
141  const bool bad_stride = 0 == tmp.stride(i) ||
142  0 != tmp.stride(i) % tmp.stride(i-1);
143  ThrowErrorMsgIf( bad_stride,
144  arg_method << " FAILED for " << *this <<
145  " WITH BAD STRIDE " <<
146  print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));;
147  }
148  }
149 
150  if (arg_init_value != NULL) {
151  //insert_restriction can be called multiple times for the same field, giving
152  //the field different lengths on different mesh-parts.
153  //We will only store one initial-value array, we need to store the one with
154  //maximum length for this field so that it can be used to initialize data
155  //for all field-restrictions. For the parts on which the field is shorter,
156  //a subset of the initial-value array will be used.
157  //
158  //We want to end up storing the longest arg_init_value array for this field.
159  //
160  //Thus, we call set_initial_value only if the current length is longer
161  //than what's already been stored.
162 
163  //length in bytes is num-scalars X sizeof-scalar:
164 
165  size_t num_scalars = 1;
166  //if rank > 0, then field is not a scalar field, so num-scalars is
167  //obtained from the stride array:
168  if (m_field_rank > 0) num_scalars = tmp.stride(m_field_rank-1);
169 
170  size_t sizeof_scalar = m_data_traits.size_of;
171  size_t nbytes = sizeof_scalar * num_scalars;
172 
173  size_t old_nbytes = 0;
174  if (get_initial_value() != NULL) {
175  old_nbytes = get_initial_value_num_bytes();
176  }
177 
178  if (nbytes > old_nbytes) {
179  set_initial_value(arg_init_value, num_scalars, nbytes);
180  }
181  }
182 
183  {
184  FieldRestrictionVector & restrs = restrictions();
185 
186  FieldRestrictionVector::iterator restr = restrs.begin();
187  FieldRestrictionVector::iterator last_restriction = restrs.end();
188 
189  restr = std::lower_bound(restr,last_restriction,tmp);
190 
191  const bool new_restriction = ( ( restr == last_restriction ) || !(*restr == tmp) );
192 
193  if ( new_restriction ) {
194  // New field restriction, verify we are not committed:
195  ThrowRequireMsg(!m_meta_data->is_commit(), "mesh MetaData has been committed.");
196  unsigned num_subsets = 0;
197  for(FieldRestrictionVector::iterator i=restrs.begin(), iend=restrs.end(); i!=iend; ++i) {
198  if (i->entity_rank() != arg_entity_rank) continue;
199 
200  const Part& partI = *m_meta_data->get_parts()[i->part_ordinal()];
201  bool found_subset = contain(arg_part.subsets(), partI);
202  if (found_subset) {
203  ThrowErrorMsgIf( i->not_equal_stride(tmp),
204  arg_method << " FAILED for " << *this << " " <<
205  print_restriction( *i, arg_entity_rank, arg_part, m_field_rank ) <<
206  " WITH INCOMPATIBLE REDECLARATION " <<
207  print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
208  *i = tmp;
209  ++num_subsets;
210  }
211 
212  bool found_superset = contain(arg_part.supersets(), partI);
213  if (found_superset) {
214  ThrowErrorMsgIf( i->not_equal_stride(tmp),
215  arg_method << " FAILED for " << *this << " " <<
216  print_restriction( *i, arg_entity_rank, arg_part, m_field_rank ) <<
217  " WITH INCOMPATIBLE REDECLARATION " <<
218  print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
219  //if there's already a restriction for a superset of this part, then
220  //there's nothing to do and we're out of here..
221  return;
222  }
223  }
224  if (num_subsets == 0) {
225  restrs.insert( restr , tmp );
226  }
227  else {
228  //if subsets were found, we replaced them with the new restriction. so now we need
229  //to sort and unique the vector, and trim it to remove any duplicates:
230  std::sort(restrs.begin(), restrs.end());
231  FieldRestrictionVector::iterator it = std::unique(restrs.begin(), restrs.end());
232  restrs.resize(it - restrs.begin());
233  }
234  }
235  else {
236  ThrowErrorMsgIf( restr->not_equal_stride(tmp),
237  arg_method << " FAILED for " << *this << " " <<
238  print_restriction( *restr, arg_entity_rank, arg_part, m_field_rank ) <<
239  " WITH INCOMPATIBLE REDECLARATION " <<
240  print_restriction( tmp, arg_entity_rank, arg_part, m_field_rank ));
241  }
242  }
243 }
244 
245 void FieldBaseImpl::insert_restriction(
246  const char * arg_method ,
247  EntityRank arg_entity_rank ,
248  const Selector & arg_selector ,
249  const unsigned * arg_stride,
250  const void* arg_init_value )
251 {
252  TraceIfWatching("stk_classic::mesh::impl::FieldBaseImpl::insert_restriction", LOG_FIELD, m_ordinal);
253 
254  FieldRestriction tmp( arg_entity_rank , arg_selector );
255 
256  {
257  unsigned i = 0 ;
258  if ( m_field_rank ) {
259  for ( i = 0 ; i < m_field_rank ; ++i ) { tmp.stride(i) = arg_stride[i] ; }
260  }
261  else { // Scalar field is 0 == m_field_rank
262  i = 1 ;
263  tmp.stride(0) = 1 ;
264  }
265  // Remaining dimensions are 1, no change to stride
266  for ( ; i < MaximumFieldDimension ; ++i ) {
267  tmp.stride(i) = tmp.stride(i-1) ;
268  }
269 
270  for ( i = 1 ; i < m_field_rank ; ++i ) {
271  const bool bad_stride = 0 == tmp.stride(i) ||
272  0 != tmp.stride(i) % tmp.stride(i-1);
273  ThrowErrorMsgIf( bad_stride,
274  arg_method << " FAILED for " << *this <<
275  " WITH BAD STRIDE!");
276  }
277  }
278 
279  if (arg_init_value != NULL) {
280  //insert_restriction can be called multiple times for the same field, giving
281  //the field different lengths on different mesh-parts.
282  //We will only store one initial-value array, we need to store the one with
283  //maximum length for this field so that it can be used to initialize data
284  //for all field-restrictions. For the parts on which the field is shorter,
285  //a subset of the initial-value array will be used.
286  //
287  //We want to end up storing the longest arg_init_value array for this field.
288  //
289  //Thus, we call set_initial_value only if the current length is longer
290  //than what's already been stored.
291 
292  //length in bytes is num-scalars X sizeof-scalar:
293 
294  size_t num_scalars = 1;
295  //if rank > 0, then field is not a scalar field, so num-scalars is
296  //obtained from the stride array:
297  if (m_field_rank > 0) num_scalars = tmp.stride(m_field_rank-1);
298 
299  size_t sizeof_scalar = m_data_traits.size_of;
300  size_t nbytes = sizeof_scalar * num_scalars;
301 
302  size_t old_nbytes = 0;
303  if (get_initial_value() != NULL) {
304  old_nbytes = get_initial_value_num_bytes();
305  }
306 
307  if (nbytes > old_nbytes) {
308  set_initial_value(arg_init_value, num_scalars, nbytes);
309  }
310  }
311 
312  {
313  FieldRestrictionVector & srvec = selector_restrictions();
314 
315  bool restriction_already_exists = false;
316  for(FieldRestrictionVector::const_iterator it=srvec.begin(), it_end=srvec.end();
317  it!=it_end; ++it) {
318  if (tmp == *it) {
319  restriction_already_exists = true;
320  if (tmp.not_equal_stride(*it)) {
321  ThrowErrorMsg("Incompatible selector field-restrictions!");
322  }
323  }
324  }
325 
326  if ( !restriction_already_exists ) {
327  // New field restriction, verify we are not committed:
328  ThrowRequireMsg(!m_meta_data->is_commit(), "mesh MetaData has been committed.");
329  srvec.push_back( tmp );
330  }
331  }
332 }
333 
334 void FieldBaseImpl::verify_and_clean_restrictions(
335  const char * arg_method ,
336  const Part& superset,
337  const Part& subset,
338  const PartVector & arg_all_parts )
339 {
340  TraceIfWatching("stk_classic::mesh::impl::FieldBaseImpl::verify_and_clean_restrictions", LOG_FIELD, m_ordinal);
341 
342  FieldRestrictionVector & restrs = restrictions();
343 
344  //Check whether both 'superset' and 'subset' are in this field's restrictions.
345  //If they are, make sure they are compatible and remove the subset restriction.
346  FieldRestrictionVector::iterator superset_restriction = restrs.end();
347  FieldRestrictionVector::iterator subset_restriction = restrs.end();
348  for (FieldRestrictionVector::iterator i = restrs.begin() ; i != restrs.end() ; ++i ) {
349  if (i->part_ordinal() == superset.mesh_meta_data_ordinal()) {
350  superset_restriction = i;
351  if (subset_restriction != restrs.end() && subset_restriction->entity_rank() == superset_restriction->entity_rank()) break;
352  }
353  if (i->part_ordinal() == subset.mesh_meta_data_ordinal()) {
354  subset_restriction = i;
355  if (superset_restriction != restrs.end() && subset_restriction->entity_rank() == superset_restriction->entity_rank()) break;
356  }
357  }
358 
359  if (superset_restriction != restrs.end() && subset_restriction != restrs.end() &&
360  superset_restriction->entity_rank() == subset_restriction->entity_rank()) {
361  ThrowErrorMsgIf( superset_restriction->not_equal_stride(*subset_restriction),
362  "Incompatible field restrictions for parts "<<superset.name()<<" and "<<subset.name());
363 
364  restrs.erase(subset_restriction);
365  }
366 }
367 
368 const void* FieldBaseImpl::get_initial_value() const
369 {
370  return m_field_states[0]->m_impl.m_initial_value;
371 }
372 
373 void* FieldBaseImpl::get_initial_value() {
374  return m_field_states[0]->m_impl.m_initial_value;
375 }
376 
377 unsigned FieldBaseImpl::get_initial_value_num_bytes() const {
378  return m_field_states[0]->m_impl.m_initial_value_num_bytes;
379 }
380 
381 void FieldBaseImpl::set_initial_value(const void* new_initial_value, unsigned num_scalars, unsigned num_bytes) {
382  void*& init_val = m_field_states[0]->m_impl.m_initial_value;
383 
384  delete [] reinterpret_cast<char*>(init_val);
385  init_val = new char[num_bytes];
386 
387  m_field_states[0]->m_impl.m_initial_value_num_bytes = num_bytes;
388 
389  m_data_traits.copy(init_val, new_initial_value, num_scalars);
390 }
391 
392 
393 //----------------------------------------------------------------------
394 //----------------------------------------------------------------------
395 // This part or any superset of this part
396 
397 const FieldRestriction &
398 FieldBaseImpl::restriction( unsigned entity_rank , const Part & part ) const
399 {
400  static const FieldRestriction empty ;
401 
402  const FieldRestrictionVector & rMap = restrictions();
403  const FieldRestrictionVector::const_iterator ie = rMap.end() ;
404  FieldRestrictionVector::const_iterator i ;
405 
406  const PartVector::const_iterator ipe = part.supersets().end();
407  PartVector::const_iterator ip = part.supersets().begin() ;
408 
409  // Start with this part:
410  //(putting static here helps performance significantly but is NOT THREAD SAFE !!!)
411  static FieldRestriction restr;
412  restr.set_entity_rank( entity_rank );
413  restr.set_part_ordinal( part.mesh_meta_data_ordinal() );
414 
415  while ( ie == ( i = find( rMap , restr ) ) && ipe != ip ) {
416  // Not found try another superset part:
417  restr.set_entity_rank( entity_rank );
418  restr.set_part_ordinal( (*ip)->mesh_meta_data_ordinal() );
419  ++ip ;
420  }
421 
422  return ie == i ? empty : *i ;
423 }
424 
425 unsigned FieldBaseImpl::max_size( unsigned entity_rank ) const
426 {
427  unsigned max = 0 ;
428 
429  const FieldRestrictionVector & rMap = restrictions();
430  const FieldRestrictionVector::const_iterator ie = rMap.end() ;
431  FieldRestrictionVector::const_iterator i = rMap.begin();
432 
433  for ( ; i != ie ; ++i ) {
434  if ( i->entity_rank() == entity_rank ) {
435  const unsigned len = m_field_rank ? i->stride( m_field_rank - 1 ) : 1 ;
436  if ( max < len ) { max = len ; }
437  }
438  }
439 
440  return max ;
441 }
442 
443 void FieldBaseImpl::set_field_states( FieldBase ** field_states)
444 {
445  TraceIfWatching("stk_classic::mesh::impl::FieldBaseImpl::set_field_states", LOG_FIELD, m_ordinal);
446 
447  for (unsigned i = 0; i < m_num_states; ++i) {
448  m_field_states[i] = field_states[i];
449  }
450 }
451 
452 //----------------------------------------------------------------------
453 
454 //----------------------------------------------------------------------
455 
456 std::ostream & operator << ( std::ostream & s , const FieldBaseImpl & field )
457 {
458  s << "FieldBaseImpl<" ;
459  s << field.data_traits().name ;
460  for ( unsigned i = 0 ; i < field.rank() ; ++i ) {
461  s << "," << field.dimension_tags()[i]->name();
462  }
463  s << ">" ;
464 
465  s << "[ name = \"" ;
466  s << field.name() ;
467  s << "\" , #states = " ;
468  s << field.number_of_states();
469  s << " ]" ;
470  return s ;
471 }
472 
473 std::ostream & print( std::ostream & s ,
474  const char * const b ,
475  const FieldBase & field )
476 {
477  const PartVector & all_parts = MetaData::get(field).get_parts();
478  const std::vector<FieldBase::Restriction> & rMap = field.restrictions();
479  s << field.name() ;
480  s << " {" ;
481  for ( FieldBase::RestrictionVector::const_iterator
482  i = rMap.begin() ; i != rMap.end() ; ++i ) {
483  s << std::endl << b << " " ;
484  i->print( s, i->entity_rank(), * all_parts[ i->part_ordinal() ], field.rank() );
485  s << std::endl;
486  }
487  s << std::endl << b << "}" ;
488  return s ;
489 }
490 
491 //----------------------------------------------------------------------
492 
493 
494 
495 
496 } // namespace impl
497 } // namespace mesh
498 } // namespace stk_classic
std::ostream & print(std::ostream &os, const std::string &indent, const Bucket &bucket)
Print the parts and entities of this bucket.
Definition: Bucket.cpp:259
bool contain(const PartVector &v, const Part &part)
Query containment within properly ordered PartVector.
Definition: Part.cpp:108
State of a field with one state.
Definition: FieldState.hpp:35
Sierra Toolkit.
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
Part * find(const PartVector &parts, const std::string &name)
Find a part by name in a collection of parts.
Definition: Part.cpp:22
EntityRank entity_rank(const EntityKey &key)
Given an entity key, return an entity type (rank).
FieldState
Enumeration of states for multi-state fields.
Definition: FieldState.hpp:34