Sierra Toolkit  Version of the Day
GearsFixture.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010, 2011 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 <cmath>
10 #include <sstream>
11 #include <stdexcept>
12 #include <limits>
13 #include <iostream>
14 #include <set>
15 
16 #include <stk_mesh/base/Types.hpp>
17 #include <stk_mesh/base/MetaData.hpp>
18 #include <stk_mesh/base/FieldData.hpp>
19 #include <stk_mesh/base/FieldParallel.hpp>
20 #include <stk_mesh/base/Entity.hpp>
21 #include <stk_mesh/base/BulkData.hpp>
22 #include <stk_mesh/base/BulkModification.hpp>
23 #include <stk_mesh/base/Selector.hpp>
24 #include <stk_mesh/base/GetEntities.hpp>
25 #include <stk_mesh/base/GetBuckets.hpp>
26 
27 #include <stk_mesh/fixtures/GearsFixture.hpp>
28 #include <stk_mesh/fixtures/Gear.hpp>
29 
30 #include <stk_mesh/fem/Stencils.hpp>
31 #include <stk_mesh/fem/FEMHelpers.hpp>
32 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
33 
34 #include <Shards_BasicTopologies.hpp>
35 
36 namespace {
37 
38 const stk_classic::mesh::EntityRank NODE_RANK = stk_classic::mesh::fem::FEMMetaData::NODE_RANK;
39 
40 const unsigned ONE_STATE = 1;
41 const unsigned TWO_STATE = 2;
42 
43 typedef shards::Hexahedron<8> Hex8 ;
44 typedef shards::Wedge<6> Wedge6 ;
45 
46 } // namespace
47 
48 namespace stk_classic {
49 namespace mesh {
50 namespace fixtures {
51 
52 GearsFixture::GearsFixture( ParallelMachine pm, size_t num_gears, GearParams gear_params)
53  : NUM_GEARS(num_gears)
54  , meta_data( SpatialDimension )
55  , bulk_data( fem::FEMMetaData::get_meta_data(meta_data) , pm )
56  , element_rank( meta_data.element_rank() )
57  , cylindrical_coord_part( meta_data.declare_part("cylindrical_coord_part", element_rank))
58  , hex_part( fem::declare_part<Hex8>(meta_data, "hex8_part"))
59  , wedge_part( fem::declare_part<Wedge6>(meta_data, "wedge6_part"))
60  , cartesian_coord_field( meta_data.declare_field<CartesianField>("coordinates", ONE_STATE))
61  , displacement_field( meta_data.declare_field<CartesianField>("displacement", TWO_STATE))
62  , translation_field( meta_data.declare_field<CartesianField>("translation", ONE_STATE))
63  , cylindrical_coord_field( meta_data.declare_field<CylindricalField>("cylindrical_coordinates", ONE_STATE))
64  , m_gears()
65  {
66 
67  put_field(
68  cartesian_coord_field,
69  NODE_RANK,
70  meta_data.universal_part(),
71  SpatialDimension
72  );
73 
74  put_field(
75  displacement_field,
76  NODE_RANK,
77  meta_data.universal_part(),
78  SpatialDimension
79  );
80 
81  put_field(
82  translation_field,
83  NODE_RANK,
84  cylindrical_coord_part,
85  SpatialDimension
86  );
87 
88  put_field(
89  cylindrical_coord_field,
90  NODE_RANK,
91  cylindrical_coord_part,
92  SpatialDimension
93  );
94 
95  m_gears.resize(NUM_GEARS);
96 
97  for ( size_t i = 0; i < NUM_GEARS; ++i) {
98  std::ostringstream oss;
99  oss << "Gear_" << i ;
100  m_gears[i] = new Gear (
101  meta_data,
102  bulk_data,
103  meta_data.declare_part(oss.str(),SpatialDimension),
104  cylindrical_coord_part,
105  hex_part,
106  wedge_part,
107  cartesian_coord_field,
108  displacement_field,
109  translation_field,
110  cylindrical_coord_field,
111  gear_params.element_size,
112  gear_params.radius_min,
113  gear_params.radius_max,
114  gear_params.height_min,
115  gear_params.height_max
116  );
117  }
118 }
119 
120 GearsFixture::~GearsFixture()
121 {
122  for( std::vector<Gear *>::iterator i = m_gears.begin();
123  i != m_gears.end();
124  ++i )
125  {
126  delete *i;
127  *i = NULL;
128  }
129 }
130 
131 void GearsFixture::generate_mesh() {
132 
133  // Parallel collective call:
134  bulk_data.modification_begin();
135 
136  const unsigned p_size = bulk_data.parallel_size();
137  const unsigned p_rank = bulk_data.parallel_rank();
138 
139  //create the gears on a line
140  for( size_t i = 0; i < m_gears.size(); ++i) {
141  if (( (i*p_size)/m_gears.size()) == p_rank) {
142  Gear & gear = get_gear(i);
143  gear.generate_gear();
144  } else {
145  // Parallel synchronization:
146  std::vector<size_t> empty_requests(meta_data.entity_rank_count(), 0);
147  EntityVector empty_entities;
148  bulk_data.generate_new_entities(empty_requests, empty_entities);
149  }
150  }
151 
152  // Parallel collective call:
153  bulk_data.modification_end();
154 
155  // Parallel collective call:
156  communicate_model_fields();
157 
158  for ( size_t i = 0 ; i < m_gears.size() ; ++i ) {
159  // Parallel collective call:
160  distribute_gear_across_processors(get_gear(i),cylindrical_coord_field);
161  }
162 
163 }
164 
165 void GearsFixture::communicate_model_fields()
166 {
167  //copy the field data to the aura nodes
168 
169  std::vector< const FieldBase *> fields;
170 
171  fields.push_back(& cartesian_coord_field);
172  fields.push_back(& translation_field);
173  fields.push_back(& cylindrical_coord_field);
174  fields.push_back(& displacement_field.field_of_state(stk_classic::mesh::StateNew));
175  fields.push_back(& displacement_field.field_of_state(stk_classic::mesh::StateOld));
176 
177  // Parallel collective call:
178  communicate_field_data(bulk_data.shared_aura(), fields);
179 }
180 
181 
182 double scale_angle_2pi(double angle) {
183  while ( angle < 0.0 ) {
184  angle += TWO_PI;
185  }
186  while ( angle >= TWO_PI) {
187  angle -= TWO_PI;
188  }
189  return angle;
190 }
191 
192 
193 void select_nodal_data(
194  GearsFixture::CylindricalField & cylindrical_coord_field,
195  Entity & element,
196  double & radius,
197  double & angle,
198  double & height
199  )
200 {
201  radius = 0.0;
202  angle = TWO_PI;
203  height = 0.0;
204  PairIterRelation node_relations = element.relations(NODE_RANK);
205  int numNodes = node_relations.second - node_relations.first;
206  for ( ; node_relations.first != node_relations.second ; ++(node_relations.first) ) {
207  Entity * node = node_relations.first->entity();
208  EntityArray<GearsFixture::CylindricalField> cylindrical_data( cylindrical_coord_field, *node);
209  radius += cylindrical_data(0);
210  angle = std::min(angle,cylindrical_data(1));
211  height += cylindrical_data(2);
212  }
213  radius /= numNodes;
214  height /= numNodes;
215 }
216 
217 // Parallel collective call:
218 void distribute_gear_across_processors(Gear & gear, GearsFixture::CylindricalField & cylindrical_coord_field)
219 {
220  BulkData & bulk_data = gear.bulk_data;
221 
222  const unsigned p_size = bulk_data.parallel_size();
223  const unsigned p_rank = bulk_data.parallel_rank();
224 
225  EntityProcVec elements_to_change_owner;
226 
227  Selector locally_owned = gear.meta_data.locally_owned_part();
228  if (p_rank == 0) {
229  BucketVector all_elements;
230  stk_classic::mesh::get_buckets(locally_owned,bulk_data.buckets(gear.meta_data.element_rank()),all_elements);
231  std::set<Entity *> node_set; // First come first serve nodal movement.
232  for (BucketVector::iterator it = all_elements.begin() ; it != all_elements.end() ; ++it) {
233  Bucket & b = **it;
234  for (size_t i=0 ; i<b.size() ; ++i) {
235  Entity * element = &b[i];
236  double radius = 0.0;
237  double angle = 0.0;
238  double height = 0.0;
239  select_nodal_data(cylindrical_coord_field, *element,radius,angle,height);
240  unsigned destination_processor_rank = destination_processor(gear,radius,angle,height,p_rank,p_size);
241  elements_to_change_owner.push_back(EntityProc(element,destination_processor_rank));
242  // Now add all related nodes to list to move to this processor:
243  PairIterRelation node_relations = element->relations(stk_classic::mesh::fem::FEMMetaData::NODE_RANK);
244  for ( ; node_relations.first != node_relations.second ; ++(node_relations.first) ) {
245  Entity * node = node_relations.first->entity();
246  if (node_set.count(node)==0) {
247  elements_to_change_owner.push_back(EntityProc(node,destination_processor_rank));
248  node_set.insert(node);
249  }
250  }
251  }
252  }
253  }
254 
255  // Parallel collective call:
256  bulk_data.modification_begin();
257  bulk_data.change_entity_owner(elements_to_change_owner);
258  // Parallel collective call:
259  bulk_data.modification_end();
260 
261  // Print out how many ended up on each processor:
262  //{
263  // BucketVector local_elements;
264  // stk_classic::mesh::get_buckets(locally_owned,bulk_data.buckets(Element),local_elements);
265  // BucketVector local_nodes;
266  // stk_classic::mesh::get_buckets(locally_owned,bulk_data.buckets(Node),local_nodes);
267  // size_t element_count = 0;
268  // for (BucketVector::iterator it = local_elements.begin() ; it != local_elements.end() ; ++it) {
269  // Bucket & b = **it;
270  // element_count += b.size();
271  // }
272  // size_t node_count = 0;
273  // for (BucketVector::iterator it = local_nodes.begin() ; it != local_nodes.end() ; ++it) {
274  // Bucket & b = **it;
275  // node_count += b.size();
276  // }
277  // std::cout << "Proc " << p_rank << ": element count = " << element_count << " node count = " << node_count << std::endl;
278  //}
279 
280 }
281 
282 
283 double floor0(double value)
284 {
285  double result = std::floor( std::fabs( value ) );
286  return (value < 0.0) ? -result : result;
287 }
288 
289 
290 void scale_p_rank(unsigned & p_rank, unsigned p_size)
291 {
292  if (p_rank >= p_size) {
293  p_rank = p_size-1;
294  }
295 }
296 
297 unsigned destination_processor(const Gear & gear, double rad, double angle, double height, unsigned p_rank, unsigned p_size)
298 {
299  unsigned result = 0;
300  // Distribute elements across angles: (not working perfectly yet)
301  angle = scale_angle_2pi(angle);
302  result = static_cast<unsigned>(floor0((angle/TWO_PI)*p_size));
303 
304  // Distribute elements across radius:
305  //result = static_cast<unsigned>(floor0((rad-gear.rad_min)/(gear.rad_max-gear.rad_min)*p_size));
306 
307  // Distribute elements across height:
308  //result = static_cast<unsigned>(floor0((height-gear.height_min)/(gear.height_max-gear.height_min)*p_size));
309 
310  // Distribute elements randomly:
311  //result = std::rand() % p_size;
312 
313  // Distribute in round-robin fashion:
314  //static unsigned dest_proc = 0;
315  //result = dest_proc++ % p_size;
316 
317  // Distribute all to processor 2:
318  //result = 1;
319 
320  scale_p_rank(result,p_size);
321 
322  //std::cout << "P"<<p_rank<<"("<<p_size<<"): (r,t,z) = ("<<rad<<", "<<angle<<", "<<height<<"), dest = " << result << std::endl;
323 
324  return result;
325 }
326 
327 } // fixtures
328 } // mesh
329 } // stk
330 
331 
332 
Newest state of a field with two states.
Definition: FieldState.hpp:36
field_type & put_field(field_type &field, EntityRank entity_rank, const Part &part, const void *init_value=NULL)
Declare a field to exist for a given entity type and Part.
std::pair< Entity *, unsigned > EntityProc
Pairing of an entity with a processor rank.
Definition: Types.hpp:111
Part & declare_part(FEMMetaData &meta_data, const std::string &name)
Declare a part with a given cell topology. This is just a convenient function that wraps FEMMetaData&#39;...
Definition: FEMHelpers.hpp:93
Sierra Toolkit.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
Previous state of a field with two states.
Definition: FieldState.hpp:38
AllSelectedBucketsRange get_buckets(const Selector &selector, const BulkData &mesh)
Definition: GetBuckets.cpp:26