Sierra Toolkit Version of the Day
Gear.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 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 <cmath>
00010 #include <stdexcept>
00011 #include <limits>
00012 #include <iostream>
00013 
00014 #include <stk_mesh/base/Types.hpp>
00015 #include <stk_mesh/base/MetaData.hpp>
00016 #include <stk_mesh/base/FieldData.hpp>
00017 #include <stk_mesh/base/FieldParallel.hpp>
00018 #include <stk_mesh/base/Entity.hpp>
00019 #include <stk_mesh/base/BulkData.hpp>
00020 #include <stk_mesh/base/BulkModification.hpp>
00021 #include <stk_mesh/base/Selector.hpp>
00022 #include <stk_mesh/base/GetBuckets.hpp>
00023 
00024 #include <stk_mesh/fixtures/Gear.hpp>
00025 
00026 namespace {
00027 
00028 const stk::mesh::EntityRank NODE_RANK = stk::mesh::fem::FEMMetaData::NODE_RANK;
00029 
00030 }
00031 
00032 namespace stk {
00033 namespace mesh {
00034 namespace fixtures {
00035 
00036 Gear::Gear(
00037     fem::FEMMetaData & meta,
00038     BulkData & bulk,
00039     Part & gear,
00040     Part & arg_cylindrical_coord_part,
00041     Part & hex,
00042     Part & wedge,
00043     CartesianField   & arg_cartesian_coord_field,
00044     CartesianField   & arg_displacement_field,
00045     CartesianField   & arg_translation_field,
00046     CylindricalField & arg_cylindrical_coord_field,
00047     double arg_element_size,
00048     double arg_radius_min,
00049     double arg_radius_max,
00050     double arg_height_min,
00051     double arg_height_max
00052     ) :
00053       element_size(arg_element_size)
00054     , rad_min(arg_radius_min)
00055     , rad_max(arg_radius_max)
00056     , height_min(arg_height_min)
00057     , height_max(arg_height_max)
00058 
00059     , angle_num(static_cast<size_t>(TWO_PI/element_size))
00060     , rad_num(static_cast<size_t>(2 +(rad_max - rad_min)/element_size))
00061     , height_num(static_cast<size_t>(2 +(height_max - height_min)/element_size))
00062 
00063     , angle_increment( TWO_PI / static_cast<double>(angle_num))
00064     , rad_increment( (rad_max-rad_min) / static_cast<double>(rad_num-1))
00065     , height_increment( (height_max-height_min) / static_cast<double>(height_num-1))
00066 
00067     , num_elements( angle_num * (rad_num-1) * (height_num-1) )
00068     , num_nodes( angle_num * (rad_num) * (height_num) )
00069 
00070     , meta_data(meta)
00071     , bulk_data(bulk)
00072 
00073     , gear_part(gear)
00074     , cylindrical_coord_part(arg_cylindrical_coord_part)
00075     , hex_part(hex)
00076     , wedge_part(wedge)
00077 
00078     , cartesian_coord_field(arg_cartesian_coord_field)
00079     , displacement_field(arg_displacement_field)
00080     , translation_field(arg_translation_field)
00081     , cylindrical_coord_field(arg_cylindrical_coord_field)
00082 {}
00083 
00084 //
00085 //-----------------------------------------------------------------------------
00086 //
00087 
00088 void Gear::populate_fields(stk::mesh::FieldState state) {
00089 
00090   //setup the cylindrical_coord_field on the hex nodes
00091   for ( size_t ir = 0 ; ir < rad_num-1; ++ir ) {
00092     const double rad = rad_min + rad_increment * ir ;
00093 
00094     for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
00095       const double angle = angle_increment * ia ;
00096 
00097       for ( size_t iz = 0 ; iz < height_num; ++iz ) {
00098         const double height = height_min + height_increment * iz ;
00099 
00100         Entity & node = get_node(iz,ir,ia);
00101 
00102         double * const cylindrical_data = field_data( cylindrical_coord_field , node );
00103         double * const translation_data = field_data( translation_field , node );
00104         double * const cartesian_data = field_data( cartesian_coord_field , node );
00105         double * const displacement_data = field_data( displacement_field.field_of_state(state) , node );
00106 
00107         cylindrical_data[0] = rad ;
00108         cylindrical_data[1] = angle ;
00109         cylindrical_data[2] = height;
00110 
00111         translation_data[0] = 0;
00112         translation_data[1] = 0;
00113         translation_data[2] = 0;
00114 
00115         displacement_data[0] = 0;
00116         displacement_data[1] = 0;
00117         displacement_data[2] = 0;
00118 
00119         cartesian_data[0] = rad * std::cos(angle);
00120         cartesian_data[1] = rad * std::sin(angle);
00121         cartesian_data[2] = height;
00122       }
00123     }
00124   }
00125 
00126   //setup the cylindrical_coord_field on the wedge nodes
00127   {
00128     size_t ir = rad_num-1;
00129     //const double rad = rad_min + rad_increment * ir ;
00130     const double rad = 1.1*rad_max;
00131 
00132     for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
00133       const double angle = angle_increment * (ia + ia +1.0)/2.0;
00134 
00135       for ( size_t iz = 0 ; iz < height_num; ++iz ) {
00136         const double height = height_min + height_increment * iz ;
00137 
00138         Entity & node = get_node(iz,ir,ia);
00139 
00140         double * const cylindrical_data = field_data( cylindrical_coord_field , node );
00141         double * const translation_data = field_data( translation_field , node );
00142         double * const cartesian_data = field_data( cartesian_coord_field , node );
00143         double * const displacement_data = field_data( displacement_field.field_of_state(state) , node );
00144 
00145         cylindrical_data[0] = rad ;
00146         cylindrical_data[1] = angle ;
00147         cylindrical_data[2] = height;
00148 
00149         translation_data[0] = 0;
00150         translation_data[1] = 0;
00151         translation_data[2] = 0;
00152 
00153         displacement_data[0] = 0;
00154         displacement_data[1] = 0;
00155         displacement_data[2] = 0;
00156 
00157         cartesian_data[0] = rad * std::cos(angle);
00158         cartesian_data[1] = rad * std::sin(angle);
00159         cartesian_data[2] = height;
00160       }
00161     }
00162   }
00163 }
00164 
00165 
00166 //
00167 //-----------------------------------------------------------------------------
00168 //
00169 
00170 void Gear::generate_gear()
00171 {
00172   const stk::mesh::EntityRank element_rank = meta_data.element_rank();
00173 
00174   std::vector<size_t> requests(meta_data.entity_rank_count(), 0);
00175   requests[NODE_RANK]     = num_nodes;
00176   requests[element_rank] = num_elements;
00177 
00178   // Parallel collective call:
00179   bulk_data.generate_new_entities(requests, gear_entities);
00180 
00181   //setup hex elements
00182   {
00183     std::vector<Part*> add_parts, remove_parts ;
00184     add_parts.push_back( & cylindrical_coord_part );
00185     add_parts.push_back( & gear_part );
00186     add_parts.push_back( & hex_part );
00187 
00188     for ( size_t ir = 0 ; ir < rad_num -2 ; ++ir ) {
00189       for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
00190         for ( size_t iz = 0 ; iz < height_num -1 ; ++iz ) {
00191 
00192           Entity & elem = get_element(iz, ir, ia);
00193           bulk_data.change_entity_parts(elem, add_parts, remove_parts);
00194 
00195           const size_t ia_1 = ( ia + 1 ) % angle_num ;
00196           const size_t ir_1 = ir + 1 ;
00197           const size_t iz_1 = iz + 1 ;
00198 
00199           Entity * node[8] ;
00200 
00201           node[0] = & get_node(iz  , ir  , ia_1 );
00202           node[1] = & get_node(iz  , ir  , ia   );
00203           node[2] = & get_node(iz_1, ir  , ia   );
00204           node[3] = & get_node(iz_1, ir  , ia_1 );
00205           node[4] = & get_node(iz  , ir_1, ia_1 );
00206           node[5] = & get_node(iz  , ir_1, ia   );
00207           node[6] = & get_node(iz_1, ir_1, ia   );
00208           node[7] = & get_node(iz_1, ir_1, ia_1 );
00209 
00210           for ( size_t j = 0 ; j < 8 ; ++j ) {
00211             bulk_data.declare_relation( elem , * node[j] , j );
00212           }
00213         }
00214       }
00215     }
00216   }
00217 
00218   //setup wedges elements
00219   {
00220     std::vector<Part*> add_parts, remove_parts ;
00221     add_parts.push_back( & cylindrical_coord_part );
00222     add_parts.push_back( & gear_part );
00223     add_parts.push_back( & wedge_part );
00224 
00225     size_t ir = rad_num-2 ;
00226     for ( size_t ia = 0 ; ia < angle_num; ++ia ) {
00227       for ( size_t iz = 0 ; iz < height_num -1 ; ++iz ) {
00228 
00229         Entity & elem = get_element(iz, ir, ia);
00230         bulk_data.change_entity_parts(elem, add_parts, remove_parts);
00231 
00232         const size_t ia_1 = ( ia + 1 ) % angle_num ;
00233         const size_t ir_1 = ir + 1 ;
00234         const size_t iz_1 = iz + 1 ;
00235 
00236         Entity * node[6] ;
00237 
00238         node[0] = & get_node(iz  , ir  , ia_1 );
00239         node[1] = & get_node(iz  , ir  , ia   );
00240         node[2] = & get_node(iz  , ir_1, ia   );
00241         node[3] = & get_node(iz_1, ir  , ia_1 );
00242         node[4] = & get_node(iz_1, ir  , ia   );
00243         node[5] = & get_node(iz_1, ir_1, ia   );
00244 
00245         for ( size_t j = 0 ; j < 6 ; ++j ) {
00246           bulk_data.declare_relation( elem , * node[j] , j );
00247         }
00248       }
00249     }
00250   }
00251 
00252   //cylindrical and cartesian coordinates
00253   //are 2 state fields.  Need to update
00254   //both states at construction
00255   populate_fields(stk::mesh::StateOld);
00256   populate_fields(stk::mesh::StateNew);
00257 
00258 }
00259 
00260 //
00261 //-----------------------------------------------------------------------------
00262 //
00263 
00264 void Gear::move( const GearMovement & data) {
00265 
00266   enum { Node = 0 };
00267 
00268   Selector select = gear_part & cylindrical_coord_part & (meta_data.locally_owned_part() | meta_data.globally_shared_part());
00269 
00270   BucketVector all_node_buckets = bulk_data.buckets(NODE_RANK);
00271 
00272   BucketVector node_buckets;
00273 
00274   get_buckets(
00275       select,
00276       all_node_buckets,
00277       node_buckets
00278       );
00279 
00280   for (BucketVector::iterator b_itr = node_buckets.begin();
00281       b_itr != node_buckets.end();
00282       ++b_itr)
00283   {
00284     Bucket & b = **b_itr;
00285     BucketArray<CylindricalField> cylindrical_data( cylindrical_coord_field, b);  // ONE STATE
00286     BucketArray<CartesianField>   translation_data( translation_field, b); // ONE STATE
00287     const BucketArray<CartesianField>   old_coordinate_data( cartesian_coord_field, b); // ONE STATE
00288     BucketArray<CartesianField>   new_displacement_data( displacement_field.field_of_state(stk::mesh::StateNew), b); // TWO STATE
00289 
00290     double new_coordinate_data[3] = {0,0,0};
00291     for (size_t i = 0; i < b.size(); ++i) {
00292       int index = i;
00293 
00294       const double   radius = cylindrical_data(0,index);
00295             double & angle  = cylindrical_data(1,index);
00296       const double   height = cylindrical_data(2,index);
00297 
00298 
00299       angle += data.rotation;
00300 
00301       if ( angle < 0.0) {
00302         angle += TWO_PI;
00303       }
00304       else if ( angle > TWO_PI) {
00305         angle -= TWO_PI;
00306       }
00307 
00308       translation_data(0,index) += data.x;
00309       translation_data(1,index) += data.y;
00310       translation_data(2,index) += data.z;
00311 
00312       new_coordinate_data[0] = translation_data(0,index) + radius * std::cos(angle);
00313       new_coordinate_data[1] = translation_data(1,index) + radius * std::sin(angle);
00314       new_coordinate_data[2] = translation_data(2,index) + height;
00315 
00316       new_displacement_data(0,index) = new_coordinate_data[0] - old_coordinate_data(0,index);
00317       new_displacement_data(1,index) = new_coordinate_data[1] - old_coordinate_data(1,index);
00318       new_displacement_data(2,index) = new_coordinate_data[2] - old_coordinate_data(2,index);
00319 
00320     }
00321   }
00322 }
00323 
00324 } // fixtures
00325 } // mesh
00326 } // stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends