Sierra Toolkit Version of the Day
BoundaryAnalysis.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 <vector>
00010 #include <set>
00011 #include <algorithm>
00012 
00013 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
00014 #include <stk_mesh/fem/FEMMetaData.hpp>
00015 #include <stk_mesh/fem/FEMHelpers.hpp>
00016 
00017 #include <stk_mesh/base/BulkData.hpp>
00018 #include <stk_mesh/base/Entity.hpp>
00019 #include <stk_mesh/base/Part.hpp>
00020 
00021 namespace stk {
00022 namespace mesh {
00023 
00024 namespace {
00025 
00026 const EntityRank NODE_RANK = fem::FEMMetaData::NODE_RANK;
00027 
00028 void filter_superimposed_entities(const Entity & entity, EntityVector & entities)
00029 {
00030   // Get the node entities for the nodes that make up the entity, we'll
00031   // use this to check for superimposed entities
00032   PairIterRelation irel = entity.relations(NODE_RANK);
00033   EntityVector entity_nodes;
00034   entity_nodes.reserve(irel.size());
00035   for ( ; !irel.empty(); ++irel ) {
00036     entity_nodes.push_back(irel->entity());
00037   }
00038   std::sort(entity_nodes.begin(), entity_nodes.end());
00039 
00040   // Make sure to remove the all superimposed entities from the list
00041   unsigned num_nodes_in_orig_entity = entity_nodes.size();
00042   EntityVector current_nodes;
00043   current_nodes.resize(num_nodes_in_orig_entity);
00044   EntityVector::iterator itr = entities.begin();
00045   while ( itr != entities.end() ) {
00046     Entity * current_entity = *itr;
00047     PairIterRelation relations = current_entity->relations(NODE_RANK);
00048 
00049     if (current_entity == &entity) {
00050       // Superimposed with self by definition
00051       itr = entities.erase(itr);
00052     }
00053     else if (relations.size() != num_nodes_in_orig_entity) {
00054       // current_entity has a different number of nodes than entity, they
00055       // cannot be superimposed
00056       ++itr;
00057     }
00058     else {
00059       for (unsigned i = 0; relations.first != relations.second;
00060            ++relations.first, ++i ) {
00061         current_nodes[i] = relations.first->entity();
00062       }
00063       std::sort(current_nodes.begin(), current_nodes.end());
00064 
00065       bool entities_are_superimposed = entity_nodes == current_nodes;
00066       if (entities_are_superimposed) {
00067         itr = entities.erase(itr);
00068       }
00069       else {
00070         ++itr;
00071       }
00072     }
00073   }
00074 }
00075 
00088 void get_adjacent_entities( const Entity & entity ,
00089                             EntityRank subcell_rank ,
00090                             unsigned subcell_identifier ,
00091                             std::vector< EntitySideComponent> & adjacent_entities)
00092 {
00093   adjacent_entities.clear();
00094 
00095   // Get nodes that make up the subcell we're looking at
00096   EntityVector subcell_nodes;
00097   const CellTopologyData * subcell_topology = fem::get_subcell_nodes(entity,
00098                                                                 subcell_rank,
00099                                                                 subcell_identifier,
00100                                                                 subcell_nodes);
00101 
00102 
00103   // Given the nodes related to the subcell, find all entities
00104   // with the same rank that have a relation to all of these nodes
00105   EntityVector potentially_adjacent_entities;
00106 
00107   get_entities_through_relations(subcell_nodes,
00108                                  entity.entity_rank(),
00109                                  potentially_adjacent_entities);
00110 
00111   // We don't want to include entities that are superimposed with
00112   // the input entity
00113   filter_superimposed_entities(entity, potentially_adjacent_entities);
00114 
00115   // Add the local ids, from the POV of the adj entitiy, to the return value.
00116   // Reverse the nodes so that the adjacent entity has them in the positive
00117   // orientation
00118   std::reverse(subcell_nodes.begin(),subcell_nodes.end());
00119 
00120   for (EntityVector::const_iterator eitr = potentially_adjacent_entities.begin();
00121        eitr != potentially_adjacent_entities.end(); ++eitr) {
00122     int local_subcell_num = fem::get_entity_subcell_id(**eitr,
00123                                                        subcell_rank,
00124                                                        subcell_topology,
00125                                                        subcell_nodes);
00126     if ( local_subcell_num != -1) {
00127       adjacent_entities.push_back(EntitySideComponent(*eitr, local_subcell_num));
00128     }
00129   }
00130 }
00131 
00132 } // unnamed namespace
00133 
00134 void boundary_analysis(const BulkData& bulk_data,
00135                        const EntityVector & entities_closure,
00136                        EntityRank closure_rank,
00137                        EntitySideVector& boundary)
00138 {
00139   const Selector locally_used = fem::FEMMetaData::get(bulk_data).locally_owned_part()
00140                               | fem::FEMMetaData::get(bulk_data).globally_shared_part();
00141 
00142   // find an iterator that points to the last item in the closure that is of a
00143   // lower-order than the closure_rank
00144   EntityVector::const_iterator itr =
00145     std::lower_bound(entities_closure.begin(),
00146                      entities_closure.end(),
00147                      EntityKey(closure_rank, 0),
00148                      EntityLess());
00149 
00150   // iterate over all the entities in the closure up to the iterator we computed above
00151   for ( ; itr != entities_closure.end() && (*itr)->entity_rank() == closure_rank; ++itr) {
00152     // some temporaries for clarity
00153     std::vector<EntitySideComponent > adjacent_entities;
00154     Entity& curr_entity = **itr;
00155     const CellTopologyData* celltopology = fem::get_cell_topology(curr_entity).getCellTopologyData();
00156     if (celltopology == NULL) {
00157       continue;
00158     }
00159 
00160     unsigned subcell_rank = closure_rank - 1;
00161     PairIterRelation relations = curr_entity.relations(NODE_RANK);
00162 
00163     // iterate over the subcells of the current entity
00164     for (unsigned nitr = 0; nitr < celltopology->subcell_count[subcell_rank]; ++nitr) {
00165       // find the entities (same rank as subcell) adjacent to this subcell
00166       unsigned subcell_identifier = nitr;
00167       get_adjacent_entities(curr_entity,
00168                             closure_rank - 1,
00169                             subcell_identifier,
00170                             adjacent_entities);
00171 
00172       // try to figure out if we want to keep ((curr_entity, subcell_identifier),
00173       //                                       (adjacent_entity, [k]))
00174       // if so, push into boundary
00175 
00176       // it is a keeper if adjacent entities[k] is not in the entities closure
00177       // AND if either curr_entity OR adjacent entities[k] is a member of the
00178       //               bulk_data.locally_used
00179 
00180       if (adjacent_entities.empty()) {
00181         EntitySide keeper;
00182         keeper.inside.entity = &curr_entity;
00183         keeper.inside.side_ordinal = subcell_identifier;
00184         keeper.outside.entity = NULL;
00185         keeper.outside.side_ordinal = 0;
00186         boundary.push_back(keeper);
00187         continue;
00188       }
00189 
00190       // iterate over adjacent entities (our neighbors)
00191       for (std::vector<EntitySideComponent >::const_iterator
00192            adj_itr = adjacent_entities.begin();
00193            adj_itr != adjacent_entities.end(); ++adj_itr) {
00194         // grab a reference to this neighbor for clarity
00195         const Entity& neighbor = *(adj_itr->entity);
00196 
00197         // see if this neighbor is in the closure, if so, not a keeper
00198         bool neighbor_is_in_closure =
00199           std::binary_search(entities_closure.begin(),
00200                              entities_closure.end(),
00201                              neighbor,
00202                              EntityLess());
00203 
00204         if (neighbor_is_in_closure) {
00205           continue;
00206         }
00207 
00208         // if neighbor or curr_entity is locally-used, add it to keeper
00209         if ( locally_used( neighbor.bucket()) || locally_used( curr_entity.bucket() ) ) {
00210           EntitySide keeper;
00211           keeper.inside.entity = &curr_entity;
00212           keeper.inside.side_ordinal = subcell_identifier;
00213           keeper.outside = *adj_itr;
00214 
00215           boundary.push_back(keeper);
00216         }
00217       }
00218     }
00219   }
00220 }
00221 
00222 
00223 }
00224 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends