Sierra Toolkit  Version of the Day
Gears.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 <stk_io/util/Gears.hpp>
10 
11 #include <math.h>
12 #include <iostream>
13 #include <limits>
14 #include <stdexcept>
15 
16 #include <stk_util/parallel/ParallelComm.hpp>
17 #include <stk_io/IossBridge.hpp>
18 
19 #include <stk_mesh/base/BulkData.hpp>
20 #include <stk_mesh/base/FieldData.hpp>
21 #include <stk_mesh/base/Comm.hpp>
22 
23 #include <stk_mesh/fem/Stencils.hpp>
24 
25 #include <Shards_BasicTopologies.hpp>
26 
27 
28 //----------------------------------------------------------------------
29 //----------------------------------------------------------------------
30 
31 namespace stk_classic {
32 namespace io {
33 namespace util {
34 
35 //----------------------------------------------------------------------
36 
37 GearFields::GearFields( stk_classic::mesh::fem::FEMMetaData & S )
38  : gear_coord( S.get_meta_data(S).declare_field<CylindricalField>( std::string("gear_coordinates") ) ),
39  model_coord( S.get_meta_data(S).declare_field<CartesianField>( std::string("coordinates") ) )
40 {
41  const stk_classic::mesh::Part & universe = S.get_meta_data(S).universal_part();
42 
43  stk_classic::mesh::put_field( gear_coord , stk_classic::mesh::fem::FEMMetaData::NODE_RANK , universe , SpatialDimension );
44  stk_classic::mesh::put_field( model_coord , stk_classic::mesh::fem::FEMMetaData::NODE_RANK , universe , SpatialDimension );
45 }
46 
47 //----------------------------------------------------------------------
48 
49 namespace {
50 
51 stk_classic::mesh::EntityId
52 identifier( size_t nthick , // Number of entities through the thickness
53  size_t nradius , // Number of entities through the radius
54  size_t iz , // Thickness index
55  size_t ir , // Radial index
56  size_t ia ) // Angle index
57 {
58  return static_cast<stk_classic::mesh::EntityId>(iz + nthick * ( ir + nradius * ia ));
59 }
60 
61 }
62 
63 
65  const std::string & name ,
66  const GearFields & gear_fields ,
67  const double center[] ,
68  const double rad_min ,
69  const double rad_max ,
70  const size_t rad_num ,
71  const double z_min ,
72  const double z_max ,
73  const size_t z_num ,
74  const size_t angle_num ,
75  const int turn_direction )
76  : m_mesh_fem_meta_data( &S ),
77  m_mesh_meta_data( S.get_meta_data(S) ),
78  m_mesh( NULL ),
79  m_gear( S.declare_part(std::string("Gear_").append(name), m_mesh_fem_meta_data->element_rank()) ),
80  m_surf( S.declare_part(std::string("Surf_").append(name), m_mesh_fem_meta_data->side_rank()) ),
81  m_gear_coord( gear_fields.gear_coord ),
82  m_model_coord(gear_fields.model_coord )
83 {
84  typedef shards::Hexahedron<> Hex ;
85  typedef shards::Quadrilateral<> Quad ;
86  enum { SpatialDimension = GearFields::SpatialDimension };
87 
90  stk_classic::mesh::fem::CellTopology hex_top (shards::getCellTopologyData<shards::Hexahedron<8> >());
91  stk_classic::mesh::fem::CellTopology quad_top(shards::getCellTopologyData<shards::Quadrilateral<4> >());
92 
93  stk_classic::mesh::fem::set_cell_topology( m_gear, hex_top );
94  stk_classic::mesh::fem::set_cell_topology( m_surf, quad_top );
95 
96  // Meshing parameters for this gear:
97 
98  const double TWO_PI = 2.0 * acos( static_cast<double>(-1.0) );
99 
100  m_center[0] = center[0] ;
101  m_center[1] = center[1] ;
102  m_center[2] = center[2] ;
103 
104  m_z_min = z_min ;
105  m_z_max = z_max ;
106  m_z_inc = (z_max - z_min) / static_cast<double>(z_num - 1);
107 
108  m_rad_min = rad_min ;
109  m_rad_max = rad_max ;
110  m_rad_inc = (rad_max - rad_min) / static_cast<double>(rad_num - 1);
111 
112  m_ang_inc = TWO_PI / static_cast<double>(angle_num) ;
113 
114  m_rad_num = rad_num ;
115  m_z_num = z_num ;
116  m_angle_num = angle_num ;
117  m_turn_dir = turn_direction ;
118 }
119 
120 
121 //----------------------------------------------------------------------
122 
123 stk_classic::mesh::Entity &Gear::create_node(const std::vector<stk_classic::mesh::Part*> & parts ,
124  stk_classic::mesh::EntityId node_id_base ,
125  size_t iz ,
126  size_t ir ,
127  size_t ia ) const
128 {
129  const double angle = m_ang_inc * ia ;
130  const double cos_angle = cos( angle );
131  const double sin_angle = sin( angle );
132 
133  const double radius = m_rad_min + m_rad_inc * ir ;
134  const double x = m_center[0] + radius * cos_angle ;
135  const double y = m_center[1] + radius * sin_angle ;
136  const double z = m_center[2] + m_z_min + m_z_inc * iz ;
137 
138  // Create the node and set the model_coordinates
139 
140  stk_classic::mesh::EntityId id_gear = identifier( m_z_num, m_rad_num, iz, ir, ia );
141  stk_classic::mesh::EntityId id = node_id_base + id_gear ;
142 
143  stk_classic::mesh::Entity & node = m_mesh->declare_entity( stk_classic::mesh::fem::FEMMetaData::NODE_RANK, id , parts );
144 
145  double * const gear_data = field_data( m_gear_coord , node );
146  double * const model_data = field_data( m_model_coord , node );
147 
148  gear_data[0] = radius ;
149  gear_data[1] = angle ;
150  gear_data[2] = z - m_center[2] ;
151 
152  model_data[0] = x ;
153  model_data[1] = y ;
154  model_data[2] = z ;
155 
156  return node;
157 }
158 
159 //----------------------------------------------------------------------
160 
161 void Gear::mesh( stk_classic::mesh::BulkData & M )
162 {
163  stk_classic::mesh::EntityRank element_rank;
164  stk_classic::mesh::EntityRank side_rank ;
165  if (m_mesh_fem_meta_data) {
166  element_rank = m_mesh_fem_meta_data->element_rank();
167  side_rank = m_mesh_fem_meta_data->side_rank();
168  }
169  else {
171  element_rank = fem.element_rank();
172  side_rank = fem.side_rank();
173  }
174 
175  M.modification_begin();
176 
177  m_mesh = & M ;
178 
179  const unsigned p_size = M.parallel_size();
180  const unsigned p_rank = M.parallel_rank();
181 
182  std::vector<size_t> counts ;
184 
185  // max_id is no longer available from comm_mesh_stats.
186  // If we assume uniform numbering from 1.., then max_id
187  // should be equal to counts...
188  const stk_classic::mesh::EntityId node_id_base = counts[ stk_classic::mesh::fem::FEMMetaData::NODE_RANK ] + 1 ;
189  const stk_classic::mesh::EntityId elem_id_base = counts[ element_rank ] + 1 ;
190 
191  const unsigned long elem_id_gear_max =
192  m_angle_num * ( m_rad_num - 1 ) * ( m_z_num - 1 );
193 
194  std::vector<stk_classic::mesh::Part*> elem_parts ;
195  std::vector<stk_classic::mesh::Part*> face_parts ;
196  std::vector<stk_classic::mesh::Part*> node_parts ;
197 
198  {
199  stk_classic::mesh::Part * const p_gear = & m_gear ;
200  stk_classic::mesh::Part * const p_surf = & m_surf ;
201 
202  elem_parts.push_back( p_gear );
203  face_parts.push_back( p_surf );
204  }
205 
206  for ( unsigned ia = 0 ; ia < m_angle_num ; ++ia ) {
207  for ( unsigned ir = 0 ; ir < m_rad_num - 1 ; ++ir ) {
208  for ( unsigned iz = 0 ; iz < m_z_num - 1 ; ++iz ) {
209 
210  stk_classic::mesh::EntityId elem_id_gear = identifier( m_z_num-1 , m_rad_num-1 , iz , ir , ia );
211 
212  if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
213 
214  stk_classic::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
215 
216  // Create the node and set the model_coordinates
217 
218  const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
219  const size_t ir_1 = ir + 1 ;
220  const size_t iz_1 = iz + 1 ;
221 
222  stk_classic::mesh::Entity * node[8] ;
223 
224  node[0] = &create_node( node_parts, node_id_base, iz , ir , ia_1 );
225  node[1] = &create_node( node_parts, node_id_base, iz_1, ir , ia_1 );
226  node[2] = &create_node( node_parts, node_id_base, iz_1, ir , ia );
227  node[3] = &create_node( node_parts, node_id_base, iz , ir , ia );
228  node[4] = &create_node( node_parts, node_id_base, iz , ir_1, ia_1 );
229  node[5] = &create_node( node_parts, node_id_base, iz_1, ir_1, ia_1 );
230  node[6] = &create_node( node_parts, node_id_base, iz_1, ir_1, ia );
231  node[7] = &create_node( node_parts, node_id_base, iz , ir_1, ia );
232 #if 0 /* VERIFY_CENTROID */
233 
234  // Centroid of the element for verification
235 
236  const double TWO_PI = 2.0 * acos( (double) -1.0 );
237  const double angle = m_ang_inc * ( 0.5 + (double) ia );
238  const double z = m_center[2] + m_z_min + m_z_inc * (0.5 + (double)iz);
239 
240  double c[3] = { 0 , 0 , 0 };
241 
242  for ( size_t j = 0 ; j < 8 ; ++j ) {
243  double * const coord_data = field_data( m_model_coord , *node[j] );
244  c[0] += coord_data[0] ;
245  c[1] += coord_data[1] ;
246  c[2] += coord_data[2] ;
247  }
248  c[0] /= 8 ; c[1] /= 8 ; c[2] /= 8 ;
249  c[0] -= m_center[0] ;
250  c[1] -= m_center[1] ;
251 
252  double val_a = atan2( c[1] , c[0] );
253  if ( val_a < 0 ) { val_a += TWO_PI ; }
254  const double err_a = angle - val_a ;
255  const double err_z = z - c[2] ;
256 
257  const double eps = 100 * std::numeric_limits<double>::epsilon();
258 
259  if ( err_z < - eps || eps < err_z ||
260  err_a < - eps || eps < err_a ) {
261  std::string msg ;
262  msg.append("problem setup element centroid error" );
263  throw std::logic_error( msg );
264  }
265 #endif
266 
268  M.declare_entity( element_rank, elem_id, elem_parts );
269 
270  for ( size_t j = 0 ; j < 8 ; ++j ) {
271  M.declare_relation( elem , * node[j] ,
272  static_cast<unsigned>(j) );
273  }
274  }
275  }
276  }
277  }
278 
279  // Array of faces on the surface
280 
281  {
282  const size_t ir = m_rad_num - 1 ;
283 
284  for ( size_t ia = 0 ; ia < m_angle_num ; ++ia ) {
285  for ( size_t iz = 0 ; iz < m_z_num - 1 ; ++iz ) {
286 
287  stk_classic::mesh::EntityId elem_id_gear =
288  identifier( m_z_num-1 , m_rad_num-1 , iz , ir-1 , ia );
289 
290  if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
291 
292  stk_classic::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
293 
294  unsigned face_ord = 5 ;
295  stk_classic::mesh::EntityId face_id = elem_id * 10 + face_ord + 1;
296 
297  stk_classic::mesh::Entity * node[4] ;
298 
299  const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
300  const size_t iz_1 = iz + 1 ;
301 
302  node[0] = &create_node( node_parts, node_id_base, iz , ir , ia_1 );
303  node[1] = &create_node( node_parts, node_id_base, iz_1, ir , ia_1 );
304  node[2] = &create_node( node_parts, node_id_base, iz_1, ir , ia );
305  node[3] = &create_node( node_parts, node_id_base, iz , ir , ia );
306 
308  M.declare_entity( side_rank, face_id, face_parts );
309 
310  for ( size_t j = 0 ; j < 4 ; ++j ) {
311  M.declare_relation( face , * node[j] ,
312  static_cast<unsigned>(j) );
313  }
314 
315  stk_classic::mesh::Entity & elem = * M.get_entity(element_rank, elem_id);
316 
317  M.declare_relation( elem , face , face_ord );
318  }
319  }
320  }
321  }
322  M.modification_begin();
323 }
324 
325 //----------------------------------------------------------------------
326 // Iterate nodes and turn them by the angle
327 
328 void Gear::turn( double /* turn_angle */ ) const
329 {
330 #if 0
331  const unsigned Length = 3 ;
332 
333  const std::vector<stk_classic::mesh::Bucket*> & ks = m_mesh->buckets( stk_classic::mesh::Node );
334  const std::vector<stk_classic::mesh::Bucket*>::const_iterator ek = ks.end();
335  std::vector<stk_classic::mesh::Bucket*>::const_iterator ik = ks.begin();
336  for ( ; ik != ek ; ++ik ) {
337  stk_classic::mesh::Bucket & k = **ik ;
338  if ( k.has_superset( m_gear ) ) {
339  const size_t n = k.size();
340  double * const bucket_gear_data = stk_classic::mesh::field_data( m_gear_coord, k.begin() );
341  double * const bucket_model_data = stk_classic::mesh::field_data( m_model_coord, k.begin() );
342 
343  for ( size_t i = 0 ; i < n ; ++i ) {
344  double * const gear_data = bucket_gear_data + i * Length ;
345  double * const model_data = bucket_model_data + i * Length ;
346  double * const current_data = bucket_current_data + i * Length ;
347  double * const disp_data = bucket_disp_data + i * Length ;
348 
349  const double radius = gear_data[0] ;
350  const double angle = gear_data[1] + turn_angle * m_turn_dir ;
351 
352  current_data[0] = m_center[0] + radius * cos( angle );
353  current_data[1] = m_center[1] + radius * sin( angle );
354  current_data[2] = m_center[2] + gear_data[2] ;
355 
356  disp_data[0] = current_data[0] - model_data[0] ;
357  disp_data[1] = current_data[1] - model_data[1] ;
358  disp_data[2] = current_data[2] - model_data[2] ;
359  }
360  }
361  }
362 #endif
363 }
364 
365 //----------------------------------------------------------------------
366 
367 }
368 }
369 }
void declare_relation(Entity &e_from, Entity &e_to, const RelationIdentifier local_id)
Declare a relation and its converse between entities in the same mesh.
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
bool comm_mesh_counts(BulkData &M, std::vector< size_t > &counts, bool local_flag)
Global counts for a mesh&#39;s entities.
Definition: Comm.cpp:26
EntityRank side_rank() const
Returns the side rank which changes depending on spatial dimension.
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
EntityRank element_rank() const
Returns the element rank which is always equal to spatial dimension.
void put_io_part_attribute(mesh::Part &part, Ioss::GroupingEntity *entity)
Definition: IossBridge.cpp:570
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.
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
Definition: BulkData.hpp:211
An application-defined subset of a problem domain.
Definition: Part.hpp:49
size_t size() const
Number of entities associated with this bucket.
Definition: Bucket.hpp:119
Part & universal_part() const
Universal subset for the problem domain. All other parts are a subset of the universal part...
Definition: MetaData.hpp:88
basic_string< char > string
string / wstring
unsigned parallel_size() const
Size of the parallel machine.
Definition: BulkData.hpp:82
bool modification_begin()
Begin a modification phase during which the mesh bulk data could become parallel inconsistent. This is a parallel synchronous call. The first time this method is called the mesh meta data is verified to be committed and parallel consistent. An exception is thrown if this verification fails.
Definition: BulkData.cpp:172
iterator begin() const
Beginning of the bucket.
Definition: Bucket.hpp:113
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
static MetaData & get_meta_data(FEMMetaData &fem_meta)
Getter for MetaData off of a FEMMetaData object.
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
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.
unsigned parallel_rank() const
Rank of the parallel machine&#39;s local processor.
Definition: BulkData.hpp:85
Entity & declare_entity(EntityRank ent_rank, EntityId ent_id, const PartVector &parts)
Create or retrieve a locally owned entity of a given rank and id.
Definition: BulkData.cpp:215
static FEMMetaData & get(const MetaData &meta)
Getter for FEMMetaData off of a MetaData object.
A container for the field data of a homogeneous collection of entities.
Definition: Bucket.hpp:94