Sierra Toolkit  Version of the Day
FindRestriction.hpp
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 
10 #ifndef stk_mesh_base_FindRestriction_hpp
11 #define stk_mesh_base_FindRestriction_hpp
12 
13 #include <stk_mesh/base/Types.hpp>
14 #include <stk_mesh/base/FieldBase.hpp>
15 #include <stk_mesh/base/FieldRestriction.hpp>
16 #include <stk_mesh/base/Selector.hpp>
17 
18 namespace stk_classic {
19 namespace mesh {
20 
21 namespace {
22 inline
23 const FieldBase::Restriction & empty_field_restriction()
24 {
25  //NOT THREAD SAFE
26  static const FieldBase::Restriction empty ;
27  return empty ;
28 }
29 }//namespace anonymous
30 
31 //Given a field and a range of parts, determine whether the field has a restriction
32 //for the part range. (Common usage is to provide the part-range from a bucket; i.e.,
33 //determine whether the field should be allocated in the bucket.)
34 template<typename PartIterator, class Compare>
35 const FieldBase::Restriction& find_restriction(const FieldBase& field,
36  EntityRank erank,
37  PartIterator pbegin,
38  PartIterator pend,
39  Compare comp)
40 {
41  GetPartIterOrdinal<PartIterator> get_part_ordinal;
42 
43  const FieldBase::Restriction & empty = empty_field_restriction();
44  const FieldBase::Restriction * restriction = & empty ;
45 
46  const std::vector<FieldBase::Restriction> & restr_vec = field.restrictions();
47  const std::vector<FieldBase::Restriction>::const_iterator iend = restr_vec.end();
48  std::vector<FieldBase::Restriction>::const_iterator ibeg = restr_vec.begin();
49 
50  //NOT THREAD SAFE
51  static FieldRestriction restr;
52 
53  for ( PartIterator it = pbegin; it != pend && iend != ibeg ; ++it ) {
54 
55  restr.set_entity_rank(erank);
56  restr.set_part_ordinal(get_part_ordinal(it));
57 
58  //lower_bound returns an iterator to either the insertion point for the
59  //'restr' argument, or to a matching restriction.
60  //It only returns the 'end' iterator if 'restr' is past the end of the
61  //vector of restrictions being searched.
62  //This depends on the input part ordinals being sorted, and on the restriction
63  //vector being sorted by part ordinal.
64 
65  ibeg = std::lower_bound( ibeg , iend , restr );
66 
67  if ( (iend != ibeg) && (*ibeg == restr) ) {
68  if ( restriction == & empty ) { restriction = & *ibeg ; }
69 
70  if ( ibeg->not_equal_stride(*restriction) ) {
71 
72  Part & p_old = MetaData::get(field).get_part( ibeg->part_ordinal() );
73  Part & p_new = MetaData::get(field).get_part( restriction->part_ordinal() );
74 
75  std::ostringstream msg ;
76  msg << " FAILED WITH INCOMPATIBLE DIMENSIONS FOR " ;
77  msg << field ;
78  msg << " Part[" << p_old.name() ;
79  msg << "] and Part[" << p_new.name() ;
80  msg << "]" ;
81 
82  ThrowErrorMsg( msg.str() );
83  }
84  }
85  }
86 
87  const std::vector<FieldBase::Restriction> & sel_res = field.selector_restrictions();
88  std::pair<PartIterator,PartIterator> part_range = std::make_pair(pbegin, pend);
89  for(std::vector<FieldBase::Restriction>::const_iterator it=sel_res.begin(), it_end=sel_res.end(); it != it_end; ++it) {
90  const Selector& selector = it->selector();
91  if (it->entity_rank() == erank && selector.apply(part_range, comp)) {
92  if (restriction == &empty) {
93  restriction = &*it;
94  }
95  if (it->not_equal_stride(*restriction)) {
96  ThrowErrorMsg("find_restriction calculation failed with different field-restriction selectors giving incompatible sizes.");
97  }
98  }
99  }
100 
101  return *restriction ;
102 }
103 
104 } // namespace mesh
105 } // namespace stk_classic
106 
107 #endif // stk_mesh_base_FindRestriction_hpp
Sierra Toolkit.