TopologyHelpers.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 
00013 #include <algorithm>
00014 #include <stdexcept>
00015 #include <sstream>
00016 #include <cassert>
00017 
00018 #include <stk_mesh/base/MetaData.hpp>
00019 #include <stk_mesh/base/Part.hpp>
00020 #include <stk_mesh/base/BulkData.hpp>
00021 #include <stk_mesh/base/Entity.hpp>
00022 #include <stk_mesh/fem/EntityTypes.hpp>
00023 #include <stk_mesh/fem/TopologyHelpers.hpp>
00024 
00025 #include <stk_util/util/StaticAssert.hpp>
00026 
00027 namespace stk {
00028 namespace mesh {
00029 
00030 //----------------------------------------------------------------------
00031 
00032 const CellTopologyData * get_cell_topology( const Part & p )
00033 { return p.attribute<CellTopologyData>(); }
00034 
00035 void set_cell_topology( Part & p , const CellTopologyData * singleton )
00036 {
00037   static const char method[] = "stk::mesh::set_cell_topology" ;
00038 
00039   MetaData & m = p.mesh_meta_data();
00040 
00041   const CellTopologyData * t = NULL ;
00042 
00043   if ( p.mesh_meta_data().entity_type_count() <= p.primary_entity_type() ||
00044        singleton == NULL ||
00045        singleton != ( t = m.declare_attribute_no_delete(p,singleton) ) ) {
00046     std::ostringstream msg ;
00047     msg << method << "( " << p.name();
00048     msg << " entity_type(" << p.primary_entity_type() << ") , " ;
00049     if ( singleton ) { msg << singleton->name ; }
00050     else             { msg << "NULL" ; }
00051     msg << " ) ERROR" ;
00052     if ( t ) { msg << "Existing topology = " << t->name ; }
00053     throw std::runtime_error( msg.str() );
00054   }
00055 }
00056 
00057 //----------------------------------------------------------------------
00058 
00059 
00060 const CellTopologyData * get_cell_topology( const Bucket & bucket )
00061 {
00062   const CellTopologyData * top = NULL ;
00063   PartVector parts ;
00064   bucket.supersets( parts );
00065 
00066   PartVector::iterator i = parts.begin() ;
00067 
00068   for ( ; NULL == top && i != parts.end() ; ++i ) {
00069     if ( bucket.entity_type() == (**i).primary_entity_type() ) {
00070       top = get_cell_topology( **i );
00071     }
00072   }
00073 
00074   bool ok = true ;
00075   
00076   for ( ; ok && i != parts.end() ; ++i ) {
00077     if ( bucket.entity_type() == (**i).primary_entity_type() ) {
00078       const CellTopologyData * const tmp = get_cell_topology( **i );
00079       ok = ((tmp == NULL) || (tmp == top)) ;
00080     }
00081   }
00082   
00083   if ( ! ok ) {
00084     std::ostringstream msg ;
00085     msg << "stk::mesh::get_cell_topology( Bucket[" ;
00086     for ( i = parts.begin() ; i != parts.end() ; ++i ) {
00087       if ( bucket.entity_type() == (**i).primary_entity_type() ) {
00088         const CellTopologyData * const tmp = get_cell_topology( **i );
00089         msg << " " << (*i)->name();
00090         if ( tmp ) { msg << "->" << tmp->name ; }
00091         msg << " ] ) FAILED WITH MULTIPLE LOCAL TOPOLOGIES" ;
00092         throw std::runtime_error( msg.str() );
00093       }
00094     }
00095   }
00096   return top ;
00097 }
00098 
00099 const CellTopologyData * get_cell_topology( const Entity & entity )
00100 { return get_cell_topology( entity.bucket() ); }
00101 
00102 //----------------------------------------------------------------------
00103 
00104 Entity & declare_element_side(
00105   BulkData & mesh ,
00106   const stk::mesh::EntityId global_side_id ,
00107   Entity & elem ,
00108   const unsigned local_side_id ,
00109   Part * part )
00110 {
00111   static const char method[] = "stk::mesh::declare_element_side" ;
00112 
00113   const CellTopologyData * const elem_top = get_cell_topology( elem );
00114 
00115   const CellTopologyData * const side_top =
00116     ( elem_top && local_side_id < elem_top->side_count )
00117     ? elem_top->side[ local_side_id ].topology : NULL ;
00118 
00119   if ( NULL == side_top ) {
00120      std::ostringstream msg ;
00121      msg << method << "( mesh , "
00122          << global_side_id
00123          << " , " ;
00124      print_entity_key( msg , mesh.mesh_meta_data() , elem.key() );
00125      msg << " , "
00126          << local_side_id
00127          << " ) FAILED" ;
00128      if ( NULL == elem_top ) {
00129        msg << " Cannot discern element topology" ;
00130      }
00131      else {
00132        msg << " Cell side id exceeds " ;
00133        msg << elem_top->name ;
00134        msg << ".side_count = " ;
00135        msg << elem_top->side_count ;
00136      }
00137      throw std::runtime_error( msg.str() );
00138    }
00139 
00140   const unsigned * const side_node_map = elem_top->side[ local_side_id ].node ;
00141 
00142   // This is dangerous if the unsigned enums are changed. Try to catch at compile time...
00143   enum { DimensionMappingAssumption_OK =
00144            StaticAssert< stk::mesh::Edge == 1 && stk::mesh::Face == 2 >::OK };
00145 
00146   PartVector add_parts ;
00147 
00148   if ( part ) { add_parts.push_back( part ); }
00149 
00150   //\TODO refactor: is 'dimension' the right thing to use for EntityType here???
00151   Entity & side = mesh.declare_entity( side_top->dimension, global_side_id, add_parts );
00152 
00153   mesh.declare_relation( elem , side , local_side_id );
00154 
00155   PairIterRelation rel = elem.relations( Node );
00156 
00157   for ( unsigned i = 0 ; i < side_top->node_count ; ++i ) {
00158     Entity & node = * rel[ side_node_map[i] ].entity();
00159     mesh.declare_relation( side , node , i );
00160   }
00161 
00162   return side ;
00163 }
00164 
00165 //----------------------------------------------------------------------
00166 
00167 bool element_side_polarity( const Entity & elem ,
00168                             const Entity & side , int local_side_id )
00169 {
00170   static const char method[] = "stk::mesh::element_side_polarity" ;
00171 
00172   const bool is_side = side.entity_type() != Edge ;
00173 
00174   const CellTopologyData * const elem_top = get_cell_topology( elem );
00175 
00176   const unsigned side_count = ! elem_top ? 0 : (
00177                                 is_side ? elem_top->side_count
00178                                         : elem_top->edge_count );
00179 
00180   if ( NULL == elem_top ||
00181        local_side_id < 0 ||
00182        static_cast<int>(side_count) <= local_side_id ) {
00183     const MetaData & meta_data = elem.bucket().mesh().mesh_meta_data();
00184     std::ostringstream msg ;
00185     msg << method ;
00186     msg << "( Element[" << elem.identifier() << "]" ;
00187     msg << " , " << meta_data.entity_type_names()[ side.entity_type() ];
00188     msg << "[" << side.identifier() << "]" ;
00189     msg << " , local_side_id = " << local_side_id << " ) FAILED: " ;
00190     if ( NULL == elem_top ) {
00191       msg << " Element has no defined topology" ;
00192     }
00193     else {
00194       msg << " Unsupported local_side_id" ;
00195     }
00196     throw std::runtime_error( msg.str() );
00197   }
00198 
00199   const CellTopologyData * const side_top =
00200     is_side ? elem_top->side[ local_side_id ].topology
00201             : elem_top->edge[ local_side_id ].topology ;
00202 
00203   const unsigned * const side_map =
00204     is_side ? elem_top->side[ local_side_id ].node 
00205             : elem_top->edge[ local_side_id ].node ;
00206 
00207   const PairIterRelation elem_nodes = elem.relations( Node );
00208   const PairIterRelation side_nodes = side.relations( Node );
00209 
00210   bool good = true ;
00211   for ( unsigned j = 0 ; good && j < side_top->node_count ; ++j ) {
00212     good = side_nodes[j].entity() == elem_nodes[ side_map[j] ].entity();
00213   }
00214   return good ;
00215 }
00216 
00217 void get_adjacent_entities( const Entity & entity ,
00218                             unsigned subcell_rank ,
00219                             unsigned subcell_identifier ,
00220                             std::vector<std::pair<Entity*, unsigned> > & adjacent_entities )
00221 {
00222   adjacent_entities.clear();
00223 
00224   // get cell topology
00225   const CellTopologyData* celltopology = get_cell_topology(entity);
00226   if (celltopology == NULL) {
00227     return;
00228   }
00229 
00230   // valid ranks fall within the dimension of the cell topology
00231   bool bad_rank = subcell_rank >= celltopology->dimension;
00232 
00233   // local id should be < number of entities of the desired type
00234   // (if you have 4 edges, their ids should be 0-3)
00235   bool bad_id = false;
00236   if (!bad_rank) {
00237     bad_id = subcell_identifier >= celltopology->subcell_count[subcell_rank];
00238   }
00239 
00240   if (bad_rank || bad_id) {
00241     std::ostringstream msg;
00242     //parallel consisent throw
00243     if (bad_rank) {
00244       msg << "stk::mesh::get_adjacent_entities( const Entity& entity, unsigned subcell_rank, ... ) subcell_rank is >= celltopology dimension\n";
00245     }
00246     else if (bad_id) {
00247       msg << "stk::mesh::get_adjacent_entities( const Entity& entity, unsigned subcell_rank, unsigned subcell_identifier, ... ) subcell_identifier is >= subcell count\n";
00248     }
00249 
00250     throw std::runtime_error(msg.str());
00251   }
00252 
00253   // For the potentially common subcell, get it's nodes and num_nodes
00254   const unsigned* nodes = celltopology->subcell[subcell_rank][subcell_identifier].node;
00255   unsigned num_nodes = celltopology->subcell[subcell_rank][subcell_identifier].topology->node_count;
00256 
00257   // Get all the nodal relationships for this entity. We are guaranteed
00258   // that, if we make it this far, the entity is guaranteed to have
00259   // some relationship to nodes (we know it is a higher-order entity
00260   // than Node).
00261   PairIterRelation relations = entity.relations(Node);
00262 
00263   // Get the node entities that are related to entity
00264   std::vector<Entity*> node_entities;
00265   for (unsigned itr = 0; itr < num_nodes; ++itr) {
00266     node_entities.push_back(relations[nodes[itr]].entity());
00267   }
00268 
00269   // Given the nodes related to the original entity, find all entities
00270   // of similar rank that have some relation to one or more of these nodes
00271   std::vector<Entity*> elements;
00272   get_entities_through_relations(node_entities, entity.entity_type(), elements);
00273 
00274   // Make sure to remove the original entity from the list
00275   bool found = false;
00276   for (std::vector<Entity*>::iterator itr = elements.begin();
00277        itr != elements.end(); ++itr) {
00278     if (*itr == &entity) {
00279       elements.erase(itr);
00280       found = true;
00281       break;
00282     }
00283   }
00284   // The original entity should be related to the nodes of its subcells
00285   assert(found);
00286 
00287   // Add the local ids, from the POV of the adj entitiy, to the return value
00288   for (std::vector<Entity*>::const_iterator itr = elements.begin();
00289        itr != elements.end(); ++itr) {
00290     unsigned local_side_num = element_local_side_id(**itr, node_entities);
00291     adjacent_entities.push_back(std::pair<Entity*, unsigned>(*itr, local_side_num));
00292   }
00293 }
00294 
00295 int element_local_side_id( const Entity & elem ,
00296                            const Entity & side )
00297 {
00298   return -1;
00299 }
00300 
00301 int element_local_side_id( const Entity & elem ,
00302                            const std::vector<Entity*>& entity_nodes )
00303 {
00304   // sort the input nodes
00305   std::vector<Entity*> sorted_entity_nodes(entity_nodes);
00306   std::sort(sorted_entity_nodes.begin(), sorted_entity_nodes.end(), EntityLess());
00307 
00308   // get topology of elem
00309   const CellTopologyData* celltopology = get_cell_topology(elem);
00310   if (celltopology == NULL) {
00311     return -1;
00312   }
00313 
00314   // get nodal relations
00315   PairIterRelation relations = elem.relations(Node);
00316 
00317   const unsigned subcell_rank = celltopology->dimension - 1;
00318 
00319   // Iterate over the subcells of elem
00320   for (unsigned itr = 0; itr < celltopology->subcell_count[subcell_rank]; ++itr) {
00321     // get the nodes for this subcell
00322     const unsigned* nodes = celltopology->subcell[subcell_rank][itr].node;
00323     unsigned num_nodes = celltopology->subcell[subcell_rank][itr].topology->node_count;
00324 
00325     // Get the nodes in the subcell ???
00326     std::vector<Entity*> node_entities;
00327     for (unsigned nitr = 0; nitr < num_nodes; ++nitr) {
00328       node_entities.push_back(relations[nodes[nitr]].entity());
00329     }
00330 
00331     // check to see if this subcell exactly contains the nodes that were passed in
00332     std::sort(node_entities.begin(), node_entities.end(), EntityLess());
00333     if (node_entities == sorted_entity_nodes) {
00334       return itr;
00335     }
00336   }
00337   return -1;
00338 }
00339 
00340 //----------------------------------------------------------------------
00341 
00342 }// namespace mesh
00343 }// namespace stk
00344 

Generated on Tue Jul 13 09:27:32 2010 for Sierra Toolkit by  doxygen 1.4.7