9 #include <stk_io/util/Gears.hpp> 16 #include <stk_util/parallel/ParallelComm.hpp> 17 #include <stk_io/IossBridge.hpp> 19 #include <stk_mesh/base/BulkData.hpp> 20 #include <stk_mesh/base/FieldData.hpp> 21 #include <stk_mesh/base/Comm.hpp> 23 #include <stk_mesh/fem/Stencils.hpp> 25 #include <Shards_BasicTopologies.hpp> 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") ) )
51 stk_classic::mesh::EntityId
52 identifier(
size_t nthick ,
58 return static_cast<stk_classic::mesh::EntityId
>(iz + nthick * ( ir + nradius * ia ));
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 ,
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) ),
81 m_gear_coord( gear_fields.gear_coord ),
82 m_model_coord(gear_fields.model_coord )
84 typedef shards::Hexahedron<> Hex ;
85 typedef shards::Quadrilateral<> Quad ;
86 enum { SpatialDimension = GearFields::SpatialDimension };
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> >());
93 stk_classic::mesh::fem::set_cell_topology( m_gear, hex_top );
94 stk_classic::mesh::fem::set_cell_topology( m_surf, quad_top );
98 const double TWO_PI = 2.0 * acos( static_cast<double>(-1.0) );
100 m_center[0] = center[0] ;
101 m_center[1] = center[1] ;
102 m_center[2] = center[2] ;
106 m_z_inc = (z_max - z_min) / static_cast<double>(z_num - 1);
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);
112 m_ang_inc = TWO_PI /
static_cast<double>(angle_num) ;
114 m_rad_num = rad_num ;
116 m_angle_num = angle_num ;
117 m_turn_dir = turn_direction ;
124 stk_classic::mesh::EntityId node_id_base ,
129 const double angle = m_ang_inc * ia ;
130 const double cos_angle = cos( angle );
131 const double sin_angle = sin( angle );
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 ;
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 ;
145 double *
const gear_data =
field_data( m_gear_coord , node );
146 double *
const model_data =
field_data( m_model_coord , node );
148 gear_data[0] = radius ;
149 gear_data[1] = angle ;
150 gear_data[2] = z - m_center[2] ;
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();
182 std::vector<size_t> 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 ;
191 const unsigned long elem_id_gear_max =
192 m_angle_num * ( m_rad_num - 1 ) * ( m_z_num - 1 );
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 ;
202 elem_parts.push_back( p_gear );
203 face_parts.push_back( p_surf );
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 ) {
210 stk_classic::mesh::EntityId elem_id_gear = identifier( m_z_num-1 , m_rad_num-1 , iz , ir , ia );
212 if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
214 stk_classic::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
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 ;
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 );
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);
240 double c[3] = { 0 , 0 , 0 };
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] ;
248 c[0] /= 8 ; c[1] /= 8 ; c[2] /= 8 ;
249 c[0] -= m_center[0] ;
250 c[1] -= m_center[1] ;
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] ;
257 const double eps = 100 * std::numeric_limits<double>::epsilon();
259 if ( err_z < - eps || eps < err_z ||
260 err_a < - eps || eps < err_a ) {
262 msg.append(
"problem setup element centroid error" );
263 throw std::logic_error( msg );
270 for (
size_t j = 0 ; j < 8 ; ++j ) {
272 static_cast<unsigned>(j) );
282 const size_t ir = m_rad_num - 1 ;
284 for (
size_t ia = 0 ; ia < m_angle_num ; ++ia ) {
285 for (
size_t iz = 0 ; iz < m_z_num - 1 ; ++iz ) {
287 stk_classic::mesh::EntityId elem_id_gear =
288 identifier( m_z_num-1 , m_rad_num-1 , iz , ir-1 , ia );
290 if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {
292 stk_classic::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;
294 unsigned face_ord = 5 ;
295 stk_classic::mesh::EntityId face_id = elem_id * 10 + face_ord + 1;
299 const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
300 const size_t iz_1 = iz + 1 ;
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 );
310 for (
size_t j = 0 ; j < 4 ; ++j ) {
312 static_cast<unsigned>(j) );
328 void Gear::turn(
double )
const 331 const unsigned Length = 3 ;
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 ) {
338 if ( k.has_superset( m_gear ) ) {
339 const size_t n = k.
size();
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 ;
349 const double radius = gear_data[0] ;
350 const double angle = gear_data[1] + turn_angle * m_turn_dir ;
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] ;
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] ;
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.
bool comm_mesh_counts(BulkData &M, std::vector< size_t > &counts, bool local_flag)
Global counts for a mesh's entities.
FieldTraits< field_type >::data_type * field_data(const field_type &f, const Bucket::iterator i)
Pointer to the field data array.
void put_io_part_attribute(mesh::Part &part, Ioss::GroupingEntity *entity)
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.
An application-defined subset of a problem domain.
size_t size() const
Number of entities associated with this bucket.
basic_string< char > string
string / wstring
unsigned parallel_size() const
Size of the parallel machine.
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.
iterator begin() const
Beginning of the bucket.
Manager for an integrated collection of entities, entity relations, and buckets of field data...
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
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'...
unsigned parallel_rank() const
Rank of the parallel machine's local processor.
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.
A container for the field data of a homogeneous collection of entities.