Sierra Toolkit Version of the Day
GeomDecomp.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2002, 2010, 2011 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 // Copyright 2001, 2002 Sandia Corporation, Albuquerque, NM.
00010 
00011 #include <limits>
00012 #include <cmath>
00013 #include <vector>
00014 #include <string>
00015 #include <cstdlib>
00016 #include <stdexcept>
00017 
00018 #include <stk_mesh/base/Types.hpp>
00019 #include <stk_mesh/base/Entity.hpp>
00020 #include <stk_mesh/base/Field.hpp>
00021 #include <stk_mesh/base/FieldData.hpp>
00022 
00023 #include <stk_util/environment/ReportHandler.hpp>
00024 
00025 #include <stk_util/parallel/Parallel.hpp>
00026 
00027 #include <stk_rebalance/GeomDecomp.hpp>
00028 
00029 static const size_t NODE_RANK = stk::mesh::fem::FEMMetaData::NODE_RANK;
00030 
00031 namespace stk {
00032 namespace rebalance {
00033 
00034 std::vector<const mesh::Entity *> GeomDecomp::entity_coordinates(const mesh::Entity                 & entity,
00035                                                                   const VectorField            & nodal_coor,
00036                                                                   std::vector<std::vector<double> >  & coordinates)
00037 {
00038   coordinates.clear();
00039   std::vector<const mesh::Entity *> mesh_nodes;
00040 
00041   const mesh::EntityRank enttype   = entity.entity_rank();
00042   if ( enttype == NODE_RANK )
00043   {
00044     throw std::runtime_error("GeomDecomp::entity_coordinates Error: Can not be called for nodal entities.");
00045   } else {
00046 
00047     // Loop over node relations in mesh entities
00048     mesh::PairIterRelation nr   = entity.relations( NODE_RANK );
00049 
00050     for ( ; nr.first != nr.second; ++nr.first )
00051     {
00052       const mesh::Relation &rel = *nr.first;
00053       if (rel.entity_rank() ==  NODE_RANK) { // %fixme: need to check for USES relation
00054         const mesh::Entity *nent = rel.entity();
00055         const unsigned ndim(field_data_size(nodal_coor, *nent)/sizeof(double)); // TODO - is there a better way to get this info?
00056         double * coor = mesh::field_data(nodal_coor, *nent);
00057         if (!coor) {
00058           throw std::runtime_error("GeomDecomp::entity_coordinates Error: The coordinate field does not exist.");
00059         }
00060         std::vector<double> temp(ndim);
00061         for ( unsigned i = 0; i < ndim; ++i ) { temp[i] = coor[i]; }
00062         coordinates.push_back(temp);
00063         mesh_nodes.push_back(nent);
00064       }
00065     }
00066   }
00067   return mesh_nodes;
00068 }
00069 
00070 std::vector<std::vector<double> > GeomDecomp::compute_entity_centroid(const mesh::Entity & entity,
00071                                                                    const VectorField & nodal_coor_ref,
00072                                                                    std::vector<double>   & centroid)
00073 {
00074   std::vector<std::vector<double> > coordinates;
00075   entity_coordinates(entity, nodal_coor_ref, coordinates);
00076 
00077   const int ndim      = coordinates.front().size();
00078   const int num_nodes = coordinates.size();
00079 
00080   centroid.resize(ndim);
00081   for (int i=0; i<ndim; ++i) { centroid[i] = 0; }
00082   for ( int j = 0; j < num_nodes; ++j ) {
00083     for ( int i = 0; i < ndim; ++i ) { centroid[i] += coordinates[j][i]; }
00084   }
00085   if (1 != num_nodes) {
00086     for (int i=0; i<ndim; ++i) { centroid[i] /= num_nodes; }
00087   }
00088   return coordinates;
00089 }
00090 
00091 namespace {
00092 void apply_rotation (std::vector<double> &coor)
00093 {
00094   // Apply slight transformation to "disalign" RCB coordinates
00095   // from the model coordinates.  This causes the RCB axis cuts
00096   // to be "disaligned" from straight lines of the model.
00097 
00098   static const double tS = 0.0001 ; /* sin( angle / 2 ), angle = 0.012 deg */
00099   static const double tC = sqrt( (double)( 1.0 - tS * tS ) );
00100   static const double tQ = tS / sqrt( (double) 3.0 );
00101   static const double t1 = tC * tC - tQ * tQ ;
00102   static const double t2 =  2.0 * tQ * ( tC + tQ );
00103   static const double t3 = -2.0 * tQ * ( tC - tQ );
00104 
00105 
00106   std::vector<double> temp(coor);
00107   const size_t nd = temp.size();
00108 
00109   // Apply minute transformation to the coordinate
00110   // to rotate the RCB axis slightly away from the model axis.
00111 
00112   if ( nd == 3 ) {
00113     coor[0] = t1 * temp[0] + t3 * temp[1] + t2 * temp[2] ;
00114     coor[1] = t2 * temp[0] + t1 * temp[1] + t3 * temp[2] ;
00115     coor[2] = t3 * temp[0] + t2 * temp[1] + t1 * temp[2] ;
00116   }
00117   else if ( nd == 2 ) {
00118     coor[0] = tC * temp[0] - tS * temp[1] ;
00119     coor[1] = tS * temp[0] + tC * temp[1] ;
00120   }
00121   else if ( nd == 1 ) {
00122     coor[0] = temp[0] ;
00123   }
00124   else {
00125     ThrowRequireMsg(false, "Spatial Dimention not 1, 2, or 3, can not apply rotation."); // Should never make it here
00126   }
00127   return;
00128 }
00129 }
00130 
00131 //: Convert a mesh entity to a single point
00132 //: in cartesian coordinates (x,y,z)
00133 void GeomDecomp::entity_to_point (const mesh::Entity            & entity,
00134                                const VectorField & nodeCoord,
00135                                std::vector<double>           & coor)
00136 {
00137   compute_entity_centroid(entity, nodeCoord, coor);
00138   apply_rotation (coor);
00139 }
00140 } // namespace rebalance
00141 } // namespace sierra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends