Sierra Toolkit  Version of the Day
GeomDecomp.cpp
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2002, 2010, 2011 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 // Copyright 2001, 2002 Sandia Corporation, Albuquerque, NM.
10 
11 #include <limits>
12 #include <cmath>
13 #include <vector>
14 #include <string>
15 #include <cstdlib>
16 #include <stdexcept>
17 
18 #include <stk_mesh/base/Types.hpp>
19 #include <stk_mesh/base/Entity.hpp>
20 #include <stk_mesh/base/Field.hpp>
21 #include <stk_mesh/base/FieldData.hpp>
22 
23 #include <stk_util/environment/ReportHandler.hpp>
24 
25 #include <stk_util/parallel/Parallel.hpp>
26 
28 
29 static const size_t NODE_RANK = stk_classic::mesh::fem::FEMMetaData::NODE_RANK;
30 
31 namespace stk_classic {
32 namespace rebalance {
33 
34 std::vector<const mesh::Entity *> GeomDecomp::entity_coordinates(const mesh::Entity & entity,
35  const VectorField & nodal_coor,
36  std::vector<std::vector<double> > & coordinates)
37 {
38  coordinates.clear();
39  std::vector<const mesh::Entity *> mesh_nodes;
40 
41  const mesh::EntityRank enttype = entity.entity_rank();
42  if ( enttype == NODE_RANK )
43  {
44  throw std::runtime_error("GeomDecomp::entity_coordinates Error: Can not be called for nodal entities.");
45  } else {
46 
47  // Loop over node relations in mesh entities
48  mesh::PairIterRelation nr = entity.relations( NODE_RANK );
49 
50  for ( ; nr.first != nr.second; ++nr.first )
51  {
52  const mesh::Relation &rel = *nr.first;
53  if (rel.entity_rank() == NODE_RANK) { // %fixme: need to check for USES relation
54  const mesh::Entity *nent = rel.entity();
55  const unsigned ndim(field_data_size(nodal_coor, *nent)/sizeof(double)); // TODO - is there a better way to get this info?
56  double * coor = mesh::field_data(nodal_coor, *nent);
57  if (!coor) {
58  throw std::runtime_error("GeomDecomp::entity_coordinates Error: The coordinate field does not exist.");
59  }
60  std::vector<double> temp(ndim);
61  for ( unsigned i = 0; i < ndim; ++i ) { temp[i] = coor[i]; }
62  coordinates.push_back(temp);
63  mesh_nodes.push_back(nent);
64  }
65  }
66  }
67  return mesh_nodes;
68 }
69 
70 std::vector<std::vector<double> > GeomDecomp::compute_entity_centroid(const mesh::Entity & entity,
71  const VectorField & nodal_coor_ref,
72  std::vector<double> & centroid)
73 {
74  std::vector<std::vector<double> > coordinates;
75  entity_coordinates(entity, nodal_coor_ref, coordinates);
76 
77  const int ndim = coordinates.front().size();
78  const int num_nodes = coordinates.size();
79 
80  centroid.resize(ndim);
81  for (int i=0; i<ndim; ++i) { centroid[i] = 0; }
82  for ( int j = 0; j < num_nodes; ++j ) {
83  for ( int i = 0; i < ndim; ++i ) { centroid[i] += coordinates[j][i]; }
84  }
85  if (1 != num_nodes) {
86  for (int i=0; i<ndim; ++i) { centroid[i] /= num_nodes; }
87  }
88  return coordinates;
89 }
90 
91 namespace {
92 void apply_rotation (std::vector<double> &coor)
93 {
94  // Apply slight transformation to "disalign" RCB coordinates
95  // from the model coordinates. This causes the RCB axis cuts
96  // to be "disaligned" from straight lines of the model.
97 
98  static const double tS = 0.0001 ; /* sin( angle / 2 ), angle = 0.012 deg */
99  static const double tC = sqrt( (double)( 1.0 - tS * tS ) );
100  static const double tQ = tS / sqrt( (double) 3.0 );
101  static const double t1 = tC * tC - tQ * tQ ;
102  static const double t2 = 2.0 * tQ * ( tC + tQ );
103  static const double t3 = -2.0 * tQ * ( tC - tQ );
104 
105 
106  std::vector<double> temp(coor);
107  const size_t nd = temp.size();
108 
109  // Apply minute transformation to the coordinate
110  // to rotate the RCB axis slightly away from the model axis.
111 
112  if ( nd == 3 ) {
113  coor[0] = t1 * temp[0] + t3 * temp[1] + t2 * temp[2] ;
114  coor[1] = t2 * temp[0] + t1 * temp[1] + t3 * temp[2] ;
115  coor[2] = t3 * temp[0] + t2 * temp[1] + t1 * temp[2] ;
116  }
117  else if ( nd == 2 ) {
118  coor[0] = tC * temp[0] - tS * temp[1] ;
119  coor[1] = tS * temp[0] + tC * temp[1] ;
120  }
121  else if ( nd == 1 ) {
122  coor[0] = temp[0] ;
123  }
124  else {
125  ThrowRequireMsg(false, "Spatial Dimention not 1, 2, or 3, can not apply rotation."); // Should never make it here
126  }
127  return;
128 }
129 }
130 
131 //: Convert a mesh entity to a single point
132 //: in cartesian coordinates (x,y,z)
134  const VectorField & nodeCoord,
135  std::vector<double> & coor)
136 {
137  compute_entity_centroid(entity, nodeCoord, coor);
138  apply_rotation (coor);
139 }
140 } // namespace rebalance
141 } // namespace sierra
Entity * entity() const
The referenced entity.
Definition: Relation.hpp:86
unsigned field_data_size(const FieldBase &f, const Bucket &k)
Size, in bytes, of the field data for each entity.
Definition: FieldData.hpp:99
bool rebalance(mesh::BulkData &bulk_data, const mesh::Selector &selector, const VectorField *coord_ref, const ScalarField *elem_weight_ref, Partition &partition, const stk_classic::mesh::EntityRank rank=stk_classic::mesh::InvalidEntityRank)
Rebalance with a Partition object.
Definition: Rebalance.cpp:164
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
static std::vector< const mesh::Entity * > entity_coordinates(const mesh::Entity &entity, const VectorField &ref, std::vector< std::vector< double > > &coordinates)
Used to return the nodal entities that compute_entity_centroid averages.
Definition: GeomDecomp.cpp:34
Field with defined data type and multi-dimensions (if any)
Definition: Field.hpp:118
A relation between two mesh entities with a relation identifier and kind .
Definition: Relation.hpp:58
Geometric support for partitioning of mesh entities.
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
Sierra Toolkit.
unsigned entity_rank() const
The rank of the referenced entity.
Definition: Relation.hpp:261
static void entity_to_point(const mesh::Entity &entity, const VectorField &ref, std::vector< double > &coor)
Convert a single mesh entity to a single point.
Definition: GeomDecomp.cpp:133
static std::vector< std::vector< double > > compute_entity_centroid(const mesh::Entity &entity, const VectorField &ref, std::vector< double > &coor)
Returns a vector of vectors containing the coordinates of the nodes that were used to compute the cen...
Definition: GeomDecomp.cpp:70