Sierra Toolkit  Version of the Day
Gear.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 <cmath>
10 #include <stdexcept>
11 #include <limits>
12 #include <iostream>
13 
14 #include <stk_mesh/base/Types.hpp>
15 #include <stk_mesh/base/MetaData.hpp>
16 #include <stk_mesh/base/FieldData.hpp>
17 #include <stk_mesh/base/FieldParallel.hpp>
18 #include <stk_mesh/base/Entity.hpp>
19 #include <stk_mesh/base/BulkData.hpp>
20 #include <stk_mesh/base/BulkModification.hpp>
21 #include <stk_mesh/base/Selector.hpp>
22 #include <stk_mesh/base/GetBuckets.hpp>
23 
24 #include <stk_mesh/fixtures/Gear.hpp>
25 
26 namespace {
27 
28 const stk_classic::mesh::EntityRank NODE_RANK = stk_classic::mesh::fem::FEMMetaData::NODE_RANK;
29 
30 }
31 
32 namespace stk_classic {
33 namespace mesh {
34 namespace fixtures {
35 
36 Gear::Gear(
37  fem::FEMMetaData & meta,
38  BulkData & bulk,
39  Part & gear,
40  Part & arg_cylindrical_coord_part,
41  Part & hex,
42  Part & wedge,
43  CartesianField & arg_cartesian_coord_field,
44  CartesianField & arg_displacement_field,
45  CartesianField & arg_translation_field,
46  CylindricalField & arg_cylindrical_coord_field,
47  double arg_element_size,
48  double arg_radius_min,
49  double arg_radius_max,
50  double arg_height_min,
51  double arg_height_max
52  ) :
53  element_size(arg_element_size)
54  , rad_min(arg_radius_min)
55  , rad_max(arg_radius_max)
56  , height_min(arg_height_min)
57  , height_max(arg_height_max)
58 
59  , angle_num(static_cast<size_t>(TWO_PI/element_size))
60  , rad_num(static_cast<size_t>(2 +(rad_max - rad_min)/element_size))
61  , height_num(static_cast<size_t>(2 +(height_max - height_min)/element_size))
62 
63  , angle_increment( TWO_PI / static_cast<double>(angle_num))
64  , rad_increment( (rad_max-rad_min) / static_cast<double>(rad_num-1))
65  , height_increment( (height_max-height_min) / static_cast<double>(height_num-1))
66 
67  , num_elements( angle_num * (rad_num-1) * (height_num-1) )
68  , num_nodes( angle_num * (rad_num) * (height_num) )
69 
70  , meta_data(meta)
71  , bulk_data(bulk)
72 
73  , gear_part(gear)
74  , cylindrical_coord_part(arg_cylindrical_coord_part)
75  , hex_part(hex)
76  , wedge_part(wedge)
77 
78  , cartesian_coord_field(arg_cartesian_coord_field)
79  , displacement_field(arg_displacement_field)
80  , translation_field(arg_translation_field)
81  , cylindrical_coord_field(arg_cylindrical_coord_field)
82 {}
83 
84 //
85 //-----------------------------------------------------------------------------
86 //
87 
88 void Gear::populate_fields(stk_classic::mesh::FieldState state) {
89 
90  //setup the cylindrical_coord_field on the hex nodes
91  for ( size_t ir = 0 ; ir < rad_num-1; ++ir ) {
92  const double rad = rad_min + rad_increment * ir ;
93 
94  for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
95  const double angle = angle_increment * ia ;
96 
97  for ( size_t iz = 0 ; iz < height_num; ++iz ) {
98  const double height = height_min + height_increment * iz ;
99 
100  Entity & node = get_node(iz,ir,ia);
101 
102  double * const cylindrical_data = field_data( cylindrical_coord_field , node );
103  double * const translation_data = field_data( translation_field , node );
104  double * const cartesian_data = field_data( cartesian_coord_field , node );
105  double * const displacement_data = field_data( displacement_field.field_of_state(state) , node );
106 
107  cylindrical_data[0] = rad ;
108  cylindrical_data[1] = angle ;
109  cylindrical_data[2] = height;
110 
111  translation_data[0] = 0;
112  translation_data[1] = 0;
113  translation_data[2] = 0;
114 
115  displacement_data[0] = 0;
116  displacement_data[1] = 0;
117  displacement_data[2] = 0;
118 
119  cartesian_data[0] = rad * std::cos(angle);
120  cartesian_data[1] = rad * std::sin(angle);
121  cartesian_data[2] = height;
122  }
123  }
124  }
125 
126  //setup the cylindrical_coord_field on the wedge nodes
127  {
128  size_t ir = rad_num-1;
129  //const double rad = rad_min + rad_increment * ir ;
130  const double rad = 1.1*rad_max;
131 
132  for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
133  const double angle = angle_increment * (ia + ia +1.0)/2.0;
134 
135  for ( size_t iz = 0 ; iz < height_num; ++iz ) {
136  const double height = height_min + height_increment * iz ;
137 
138  Entity & node = get_node(iz,ir,ia);
139 
140  double * const cylindrical_data = field_data( cylindrical_coord_field , node );
141  double * const translation_data = field_data( translation_field , node );
142  double * const cartesian_data = field_data( cartesian_coord_field , node );
143  double * const displacement_data = field_data( displacement_field.field_of_state(state) , node );
144 
145  cylindrical_data[0] = rad ;
146  cylindrical_data[1] = angle ;
147  cylindrical_data[2] = height;
148 
149  translation_data[0] = 0;
150  translation_data[1] = 0;
151  translation_data[2] = 0;
152 
153  displacement_data[0] = 0;
154  displacement_data[1] = 0;
155  displacement_data[2] = 0;
156 
157  cartesian_data[0] = rad * std::cos(angle);
158  cartesian_data[1] = rad * std::sin(angle);
159  cartesian_data[2] = height;
160  }
161  }
162  }
163 }
164 
165 
166 //
167 //-----------------------------------------------------------------------------
168 //
169 
170 void Gear::generate_gear()
171 {
172  const stk_classic::mesh::EntityRank element_rank = meta_data.element_rank();
173 
174  std::vector<size_t> requests(meta_data.entity_rank_count(), 0);
175  requests[NODE_RANK] = num_nodes;
176  requests[element_rank] = num_elements;
177 
178  // Parallel collective call:
179  bulk_data.generate_new_entities(requests, gear_entities);
180 
181  //setup hex elements
182  {
183  std::vector<Part*> add_parts, remove_parts ;
184  add_parts.push_back( & cylindrical_coord_part );
185  add_parts.push_back( & gear_part );
186  add_parts.push_back( & hex_part );
187 
188  for ( size_t ir = 0 ; ir < rad_num -2 ; ++ir ) {
189  for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
190  for ( size_t iz = 0 ; iz < height_num -1 ; ++iz ) {
191 
192  Entity & elem = get_element(iz, ir, ia);
193  bulk_data.change_entity_parts(elem, add_parts, remove_parts);
194 
195  const size_t ia_1 = ( ia + 1 ) % angle_num ;
196  const size_t ir_1 = ir + 1 ;
197  const size_t iz_1 = iz + 1 ;
198 
199  Entity * node[8] ;
200 
201  node[0] = & get_node(iz , ir , ia_1 );
202  node[1] = & get_node(iz , ir , ia );
203  node[2] = & get_node(iz_1, ir , ia );
204  node[3] = & get_node(iz_1, ir , ia_1 );
205  node[4] = & get_node(iz , ir_1, ia_1 );
206  node[5] = & get_node(iz , ir_1, ia );
207  node[6] = & get_node(iz_1, ir_1, ia );
208  node[7] = & get_node(iz_1, ir_1, ia_1 );
209 
210  for ( size_t j = 0 ; j < 8 ; ++j ) {
211  bulk_data.declare_relation( elem , * node[j] , j );
212  }
213  }
214  }
215  }
216  }
217 
218  //setup wedges elements
219  {
220  std::vector<Part*> add_parts, remove_parts ;
221  add_parts.push_back( & cylindrical_coord_part );
222  add_parts.push_back( & gear_part );
223  add_parts.push_back( & wedge_part );
224 
225  size_t ir = rad_num-2 ;
226  for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
227  for ( size_t iz = 0 ; iz < height_num -1 ; ++iz ) {
228 
229  Entity & elem = get_element(iz, ir, ia);
230  bulk_data.change_entity_parts(elem, add_parts, remove_parts);
231 
232  const size_t ia_1 = ( ia + 1 ) % angle_num ;
233  const size_t ir_1 = ir + 1 ;
234  const size_t iz_1 = iz + 1 ;
235 
236  Entity * node[6] ;
237 
238  node[0] = & get_node(iz , ir , ia_1 );
239  node[1] = & get_node(iz , ir , ia );
240  node[2] = & get_node(iz , ir_1, ia );
241  node[3] = & get_node(iz_1, ir , ia_1 );
242  node[4] = & get_node(iz_1, ir , ia );
243  node[5] = & get_node(iz_1, ir_1, ia );
244 
245  for ( size_t j = 0 ; j < 6 ; ++j ) {
246  bulk_data.declare_relation( elem , * node[j] , j );
247  }
248  }
249  }
250  }
251 
252  //cylindrical and cartesian coordinates
253  //are 2 state fields. Need to update
254  //both states at construction
255  populate_fields(stk_classic::mesh::StateOld);
256  populate_fields(stk_classic::mesh::StateNew);
257 
258 }
259 
260 //
261 //-----------------------------------------------------------------------------
262 //
263 
264 void Gear::move( const GearMovement & data) {
265 
266  enum { Node = 0 };
267 
268  Selector select = gear_part & cylindrical_coord_part & (meta_data.locally_owned_part() | meta_data.globally_shared_part());
269 
270  BucketVector all_node_buckets = bulk_data.buckets(NODE_RANK);
271 
272  BucketVector node_buckets;
273 
274  get_buckets(
275  select,
276  all_node_buckets,
277  node_buckets
278  );
279 
280  for (BucketVector::iterator b_itr = node_buckets.begin();
281  b_itr != node_buckets.end();
282  ++b_itr)
283  {
284  Bucket & b = **b_itr;
285  BucketArray<CylindricalField> cylindrical_data( cylindrical_coord_field, b); // ONE STATE
286  BucketArray<CartesianField> translation_data( translation_field, b); // ONE STATE
287  const BucketArray<CartesianField> old_coordinate_data( cartesian_coord_field, b); // ONE STATE
288  BucketArray<CartesianField> new_displacement_data( displacement_field.field_of_state(stk_classic::mesh::StateNew), b); // TWO STATE
289 
290  double new_coordinate_data[3] = {0,0,0};
291  for (size_t i = 0; i < b.size(); ++i) {
292  int index = i;
293 
294  const double radius = cylindrical_data(0,index);
295  double & angle = cylindrical_data(1,index);
296  const double height = cylindrical_data(2,index);
297 
298 
299  angle += data.rotation;
300 
301  if ( angle < 0.0) {
302  angle += TWO_PI;
303  }
304  else if ( angle > TWO_PI) {
305  angle -= TWO_PI;
306  }
307 
308  translation_data(0,index) += data.x;
309  translation_data(1,index) += data.y;
310  translation_data(2,index) += data.z;
311 
312  new_coordinate_data[0] = translation_data(0,index) + radius * std::cos(angle);
313  new_coordinate_data[1] = translation_data(1,index) + radius * std::sin(angle);
314  new_coordinate_data[2] = translation_data(2,index) + height;
315 
316  new_displacement_data(0,index) = new_coordinate_data[0] - old_coordinate_data(0,index);
317  new_displacement_data(1,index) = new_coordinate_data[1] - old_coordinate_data(1,index);
318  new_displacement_data(2,index) = new_coordinate_data[2] - old_coordinate_data(2,index);
319 
320  }
321  }
322 }
323 
324 } // fixtures
325 } // mesh
326 } // stk
FieldTraits< field_type >::data_type * field_data(const field_type &f, const Bucket::iterator i)
Pointer to the field data array.
Definition: FieldData.hpp:116
Newest state of a field with two states.
Definition: FieldState.hpp:36
Sierra Toolkit.
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
FieldState
Enumeration of states for multi-state fields.
Definition: FieldState.hpp:34