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