Sierra Toolkit Version of the Day
Gears.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010, 2011 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #include <stk_io/util/Gears.hpp>
00010 
00011 #include <math.h>
00012 #include <iostream>
00013 #include <limits>
00014 #include <stdexcept>
00015 
00016 #include <stk_util/parallel/ParallelComm.hpp>
00017 #include <stk_io/IossBridge.hpp>
00018 
00019 #include <stk_mesh/base/BulkData.hpp>
00020 #include <stk_mesh/base/FieldData.hpp>
00021 #include <stk_mesh/base/Comm.hpp>
00022 
00023 #include <stk_mesh/fem/Stencils.hpp>
00024 
00025 #include <Shards_BasicTopologies.hpp>
00026 
00027 
00028 //----------------------------------------------------------------------
00029 //----------------------------------------------------------------------
00030 
00031 namespace stk {
00032 namespace io {
00033 namespace util {
00034 
00035 //----------------------------------------------------------------------
00036 
00037 GearFields::GearFields( stk::mesh::fem::FEMMetaData & S )
00038   : gear_coord(          S.get_meta_data(S).declare_field<CylindricalField>( std::string("gear_coordinates") ) ),
00039     model_coord(         S.get_meta_data(S).declare_field<CartesianField>( std::string("coordinates") ) )
00040 {
00041   const stk::mesh::Part & universe = S.get_meta_data(S).universal_part();
00042 
00043   stk::mesh::put_field( gear_coord    , stk::mesh::fem::FEMMetaData::NODE_RANK , universe , SpatialDimension );
00044   stk::mesh::put_field( model_coord   , stk::mesh::fem::FEMMetaData::NODE_RANK , universe , SpatialDimension );
00045 }
00046 
00047 //----------------------------------------------------------------------
00048 
00049 namespace {
00050 
00051 stk::mesh::EntityId
00052 identifier( size_t nthick ,  // Number of entities through the thickness
00053         size_t nradius , // Number of entities through the radius
00054         size_t iz ,      // Thickness index
00055         size_t ir ,      // Radial index
00056         size_t ia )      // Angle index
00057 {
00058   return static_cast<stk::mesh::EntityId>(iz + nthick * ( ir + nradius * ia ));
00059 }
00060 
00061 }
00062 
00063 
00064 Gear::Gear( stk::mesh::fem::FEMMetaData & S ,
00065             const std::string & name ,
00066             const GearFields & gear_fields ,
00067             const double center[] ,
00068             const double rad_min ,
00069             const double rad_max ,
00070             const size_t rad_num ,
00071             const double z_min ,
00072             const double z_max ,
00073             const size_t z_num ,
00074             const size_t angle_num ,
00075             const int      turn_direction )
00076   : m_mesh_fem_meta_data( &S ),
00077     m_mesh_meta_data( S.get_meta_data(S) ),
00078     m_mesh( NULL ),
00079     m_gear( S.declare_part(std::string("Gear_").append(name), m_mesh_fem_meta_data->element_rank()) ),
00080     m_surf( S.declare_part(std::string("Surf_").append(name), m_mesh_fem_meta_data->side_rank()) ),
00081     m_gear_coord( gear_fields.gear_coord ),
00082     m_model_coord(gear_fields.model_coord )
00083 {
00084   typedef shards::Hexahedron<> Hex ;
00085   typedef shards::Quadrilateral<> Quad ;
00086   enum { SpatialDimension = GearFields::SpatialDimension };
00087 
00088   stk::io::put_io_part_attribute(m_gear);
00089   stk::io::put_io_part_attribute(m_surf);
00090   stk::mesh::fem::CellTopology hex_top (shards::getCellTopologyData<shards::Hexahedron<8> >());
00091   stk::mesh::fem::CellTopology quad_top(shards::getCellTopologyData<shards::Quadrilateral<4> >());
00092 
00093   stk::mesh::fem::set_cell_topology( m_gear, hex_top );
00094   stk::mesh::fem::set_cell_topology( m_surf, quad_top );
00095 
00096   // Meshing parameters for this gear:
00097 
00098   const double TWO_PI = 2.0 * acos( static_cast<double>(-1.0) );
00099 
00100   m_center[0] = center[0] ;
00101   m_center[1] = center[1] ;
00102   m_center[2] = center[2] ;
00103 
00104   m_z_min     = z_min ;
00105   m_z_max     = z_max ;
00106   m_z_inc     = (z_max - z_min) / static_cast<double>(z_num - 1);
00107 
00108   m_rad_min   = rad_min ;
00109   m_rad_max   = rad_max ;
00110   m_rad_inc   = (rad_max - rad_min) / static_cast<double>(rad_num - 1);
00111 
00112   m_ang_inc   = TWO_PI / static_cast<double>(angle_num) ;
00113 
00114   m_rad_num   = rad_num ;
00115   m_z_num     = z_num ;
00116   m_angle_num = angle_num ;
00117   m_turn_dir  = turn_direction ;
00118 }
00119 
00120 
00121 //----------------------------------------------------------------------
00122 
00123 stk::mesh::Entity &Gear::create_node(const std::vector<stk::mesh::Part*> & parts ,
00124                                      stk::mesh::EntityId node_id_base ,
00125                                      size_t iz ,
00126                                      size_t ir ,
00127                                      size_t ia ) const
00128 {
00129   const double angle     = m_ang_inc * ia ;
00130   const double cos_angle = cos( angle );
00131   const double sin_angle = sin( angle );
00132 
00133   const double radius = m_rad_min + m_rad_inc * ir ;
00134   const double x = m_center[0] + radius * cos_angle ;
00135   const double y = m_center[1] + radius * sin_angle ;
00136   const double z = m_center[2] + m_z_min + m_z_inc * iz ;
00137 
00138   // Create the node and set the model_coordinates
00139 
00140   stk::mesh::EntityId id_gear = identifier( m_z_num, m_rad_num, iz, ir, ia );
00141   stk::mesh::EntityId id = node_id_base + id_gear ;
00142 
00143   stk::mesh::Entity & node = m_mesh->declare_entity( stk::mesh::fem::FEMMetaData::NODE_RANK, id , parts );
00144 
00145   double * const gear_data    = field_data( m_gear_coord , node );
00146   double * const model_data   = field_data( m_model_coord , node );
00147 
00148   gear_data[0] = radius ;
00149   gear_data[1] = angle ;
00150   gear_data[2] = z - m_center[2] ;
00151 
00152   model_data[0] = x ;
00153   model_data[1] = y ;
00154   model_data[2] = z ;
00155 
00156   return node;
00157 }
00158 
00159 //----------------------------------------------------------------------
00160 
00161 void Gear::mesh( stk::mesh::BulkData & M )
00162 {
00163   stk::mesh::EntityRank element_rank;
00164   stk::mesh::EntityRank side_rank    ;
00165   if (m_mesh_fem_meta_data) {
00166     element_rank = m_mesh_fem_meta_data->element_rank();
00167     side_rank    = m_mesh_fem_meta_data->side_rank();
00168   }
00169   else {
00170     stk::mesh::fem::FEMMetaData &fem = stk::mesh::fem::FEMMetaData::get(M);
00171     element_rank = fem.element_rank();
00172     side_rank    = fem.side_rank();
00173   }
00174 
00175   M.modification_begin();
00176 
00177   m_mesh = & M ;
00178 
00179   const unsigned p_size = M.parallel_size();
00180   const unsigned p_rank = M.parallel_rank();
00181 
00182   std::vector<size_t> counts ;
00183   stk::mesh::comm_mesh_counts(M, counts);
00184 
00185   // max_id is no longer available from comm_mesh_stats.
00186   // If we assume uniform numbering from 1.., then max_id
00187   // should be equal to counts...
00188   const stk::mesh::EntityId node_id_base = counts[ stk::mesh::fem::FEMMetaData::NODE_RANK ] + 1 ;
00189   const stk::mesh::EntityId elem_id_base = counts[ element_rank ] + 1 ;
00190 
00191   const unsigned long elem_id_gear_max =
00192     m_angle_num * ( m_rad_num - 1 ) * ( m_z_num - 1 );
00193 
00194   std::vector<stk::mesh::Part*> elem_parts ;
00195   std::vector<stk::mesh::Part*> face_parts ;
00196   std::vector<stk::mesh::Part*> node_parts ;
00197 
00198   {
00199     stk::mesh::Part * const p_gear = & m_gear ;
00200     stk::mesh::Part * const p_surf = & m_surf ;
00201 
00202     elem_parts.push_back( p_gear );
00203     face_parts.push_back( p_surf );
00204   }
00205 
00206   for ( unsigned ia = 0 ; ia < m_angle_num ; ++ia ) {
00207     for ( unsigned ir = 0 ; ir < m_rad_num - 1 ; ++ir ) {
00208       for ( unsigned iz = 0 ; iz < m_z_num - 1 ; ++iz ) {
00209 
00210         stk::mesh::EntityId elem_id_gear = identifier( m_z_num-1 , m_rad_num-1 , iz , ir , ia );
00211 
00212         if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
00213 
00214           stk::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
00215 
00216           // Create the node and set the model_coordinates
00217 
00218           const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
00219           const size_t ir_1 = ir + 1 ;
00220           const size_t iz_1 = iz + 1 ;
00221 
00222           stk::mesh::Entity * node[8] ;
00223 
00224           node[0] = &create_node( node_parts, node_id_base, iz  , ir  , ia_1 );
00225           node[1] = &create_node( node_parts, node_id_base, iz_1, ir  , ia_1 );
00226           node[2] = &create_node( node_parts, node_id_base, iz_1, ir  , ia   );
00227           node[3] = &create_node( node_parts, node_id_base, iz  , ir  , ia   );
00228           node[4] = &create_node( node_parts, node_id_base, iz  , ir_1, ia_1 );
00229           node[5] = &create_node( node_parts, node_id_base, iz_1, ir_1, ia_1 );
00230           node[6] = &create_node( node_parts, node_id_base, iz_1, ir_1, ia   );
00231           node[7] = &create_node( node_parts, node_id_base, iz  , ir_1, ia   );
00232 #if 0 /* VERIFY_CENTROID */
00233 
00234           // Centroid of the element for verification
00235 
00236           const double TWO_PI = 2.0 * acos( (double) -1.0 );
00237           const double angle = m_ang_inc * ( 0.5 + (double) ia );
00238           const double z = m_center[2] + m_z_min + m_z_inc * (0.5 + (double)iz);
00239 
00240           double c[3] = { 0 , 0 , 0 };
00241 
00242           for ( size_t j = 0 ; j < 8 ; ++j ) {
00243             double * const coord_data = field_data( m_model_coord , *node[j] );
00244             c[0] += coord_data[0] ;
00245             c[1] += coord_data[1] ;
00246             c[2] += coord_data[2] ;
00247           }
00248           c[0] /= 8 ; c[1] /= 8 ; c[2] /= 8 ;
00249           c[0] -= m_center[0] ;
00250           c[1] -= m_center[1] ;
00251 
00252           double val_a = atan2( c[1] , c[0] );
00253           if ( val_a < 0 ) { val_a += TWO_PI ; }
00254           const double err_a = angle - val_a ;
00255           const double err_z = z - c[2] ;
00256 
00257           const double eps = 100 * std::numeric_limits<double>::epsilon();
00258 
00259           if ( err_z < - eps || eps < err_z ||
00260                err_a < - eps || eps < err_a ) {
00261             std::string msg ;
00262             msg.append("problem setup element centroid error" );
00263             throw std::logic_error( msg );
00264           }
00265 #endif
00266 
00267           stk::mesh::Entity & elem =
00268             M.declare_entity( element_rank, elem_id, elem_parts );
00269 
00270           for ( size_t j = 0 ; j < 8 ; ++j ) {
00271             M.declare_relation( elem , * node[j] ,
00272                                 static_cast<unsigned>(j) );
00273           }
00274         }
00275       }
00276     }
00277   }
00278 
00279   // Array of faces on the surface
00280 
00281   {
00282     const size_t ir = m_rad_num - 1 ;
00283 
00284     for ( size_t ia = 0 ; ia < m_angle_num ; ++ia ) {
00285       for ( size_t iz = 0 ; iz < m_z_num - 1 ; ++iz ) {
00286 
00287         stk::mesh::EntityId elem_id_gear =
00288           identifier( m_z_num-1 , m_rad_num-1 , iz , ir-1 , ia );
00289 
00290         if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
00291 
00292           stk::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
00293 
00294           unsigned face_ord = 5 ;
00295           stk::mesh::EntityId face_id = elem_id * 10 + face_ord + 1;
00296 
00297           stk::mesh::Entity * node[4] ;
00298 
00299           const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
00300           const size_t iz_1 = iz + 1 ;
00301 
00302           node[0] = &create_node( node_parts, node_id_base, iz  , ir  , ia_1 );
00303           node[1] = &create_node( node_parts, node_id_base, iz_1, ir  , ia_1 );
00304           node[2] = &create_node( node_parts, node_id_base, iz_1, ir  , ia   );
00305           node[3] = &create_node( node_parts, node_id_base, iz  , ir  , ia   );
00306 
00307           stk::mesh::Entity & face =
00308             M.declare_entity( side_rank, face_id, face_parts );
00309 
00310           for ( size_t j = 0 ; j < 4 ; ++j ) {
00311             M.declare_relation( face , * node[j] ,
00312                                 static_cast<unsigned>(j) );
00313           }
00314 
00315           stk::mesh::Entity & elem = * M.get_entity(element_rank, elem_id);
00316 
00317           M.declare_relation( elem , face , face_ord );
00318         }
00319       }
00320     }
00321   }
00322   M.modification_begin();
00323 }
00324 
00325 //----------------------------------------------------------------------
00326 // Iterate nodes and turn them by the angle
00327 
00328 void Gear::turn( double /* turn_angle */ ) const
00329 {
00330 #if 0
00331   const unsigned Length = 3 ;
00332 
00333   const std::vector<stk::mesh::Bucket*> & ks = m_mesh->buckets( stk::mesh::Node );
00334   const std::vector<stk::mesh::Bucket*>::const_iterator ek = ks.end();
00335   std::vector<stk::mesh::Bucket*>::const_iterator ik = ks.begin();
00336   for ( ; ik != ek ; ++ik ) {
00337     stk::mesh::Bucket & k = **ik ;
00338     if ( k.has_superset( m_gear ) ) {
00339       const size_t n = k.size();
00340       double * const bucket_gear_data    = stk::mesh::field_data( m_gear_coord,    k.begin() );
00341       double * const bucket_model_data   = stk::mesh::field_data( m_model_coord,   k.begin() );
00342 
00343       for ( size_t i = 0 ; i < n ; ++i ) {
00344         double * const gear_data    = bucket_gear_data    + i * Length ;
00345         double * const model_data   = bucket_model_data   + i * Length ;
00346         double * const current_data = bucket_current_data + i * Length ;
00347         double * const disp_data    = bucket_disp_data    + i * Length ;
00348 
00349         const double radius = gear_data[0] ;
00350         const double angle  = gear_data[1] + turn_angle * m_turn_dir ;
00351 
00352         current_data[0] = m_center[0] + radius * cos( angle );
00353         current_data[1] = m_center[1] + radius * sin( angle );
00354         current_data[2] = m_center[2] + gear_data[2] ;
00355 
00356         disp_data[0] = current_data[0] - model_data[0] ;
00357         disp_data[1] = current_data[1] - model_data[1] ;
00358         disp_data[2] = current_data[2] - model_data[2] ;
00359       }
00360     }
00361   }
00362 #endif
00363 }
00364 
00365 //----------------------------------------------------------------------
00366 
00367 }
00368 }
00369 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines