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 <stk_mesh/fem/BoundaryAnalysis.hpp>
00010 
00011 #include <stk_mesh/base/BulkData.hpp>
00012 #include <stk_mesh/base/Entity.hpp>
00013 #include <stk_mesh/base/Part.hpp>
00014 #include <stk_mesh/fem/TopologyHelpers.hpp>
00015 
00016 namespace stk {
00017 namespace mesh {
00018 
00019 void boundary_analysis(const BulkData& bulk_data,
00020                        const std::vector< Entity *> & entities_closure,
00021                        unsigned closure_rank,
00022                        EntitySideVector& boundary)
00023 {
00024   Part& locally_used_part = bulk_data.mesh_meta_data().locally_used_part();
00025 
00026   // find an iterator that points to the last item in the closure that is of a 
00027   // lower-order than the closure_rank
00028   std::vector<Entity*>::const_iterator itr = std::lower_bound(entities_closure.begin(),
00029                                                               entities_closure.end(),
00030                                                               EntityKey(closure_rank, 0),
00031                                                               EntityLess());
00032 
00033   // iterate over all the entities in the closure up to the iterator we computed above
00034   for ( ; itr != entities_closure.end() && (*itr)->entity_type() == closure_rank; ++itr) {
00035     // some temporaries for clarity
00036     std::vector<std::pair<Entity*, unsigned> > adjacent_entities;
00037     Entity& curr_entity = **itr;
00038     const CellTopologyData* celltopology = get_cell_topology(curr_entity);
00039     if (celltopology == NULL) {
00040       continue;
00041     }
00042 
00043     unsigned subcell_rank = closure_rank - 1;
00044     PairIterRelation relations = curr_entity.relations(Node);
00045 
00046     // iterate over the subcells of the current entity
00047     for (unsigned nitr = 0; nitr < celltopology->subcell_count[subcell_rank]; ++nitr) {
00048       // find the entities (same rank as subcell) adjacent to this subcell
00049       unsigned subcell_identifier = nitr;
00050       get_adjacent_entities(curr_entity,
00051                             closure_rank - 1,
00052                             subcell_identifier,
00053                             adjacent_entities);
00054       
00055       // try to figure out if we want to keep ((curr_entity, subcell_identifier), 
00056       //                                       (adjacent_entity, [k]))
00057       // if so, push into boundary
00058       
00059       // it is a keeper if adjacent entities[k] is not in the entities closure
00060       // AND if either curr_entity OR adjacent entities[k] is a member of the 
00061       //               bulk_data.locally_used
00062 
00063       if (adjacent_entities.empty()) {
00064         EntitySide keeper;
00065         keeper.first.first = &curr_entity;
00066         keeper.first.second = subcell_identifier;
00067         keeper.second.first = NULL;
00068         keeper.second.second = 0;
00069         boundary.push_back(keeper);
00070         continue;
00071       }
00072 
00073       // iterate over adjacent entities (our neighbors)
00074       for (std::vector<std::pair<Entity*, unsigned> >::const_iterator adj_itr = adjacent_entities.begin(); adj_itr != adjacent_entities.end(); ++adj_itr) {
00075         // grab a reference to this neighbor for clarity
00076         const Entity& neighbor = *(adj_itr->first);
00077 
00078         // see if this neighbor is in the closure, if so, not a keeper
00079         std::vector<Entity*>::const_iterator search_itr = std::lower_bound(entities_closure.begin(), entities_closure.end(), neighbor, EntityLess());
00080         EntityEqual eq;
00081         if (search_itr != entities_closure.end() && eq(**search_itr, neighbor)) {
00082           continue;
00083         }
00084 
00085         // if neighbor or curr_entity is locally-used, add it to keeper
00086         if (has_superset(neighbor.bucket(), locally_used_part) ||
00087             has_superset(curr_entity.bucket(), locally_used_part)) {
00088           EntitySide keeper;
00089           keeper.first.first = &curr_entity;
00090           keeper.first.second = subcell_identifier;
00091           keeper.second = *adj_itr;
00092 
00093           boundary.push_back(keeper);
00094         }
00095       }
00096     }
00097   }
00098 }
00099 
00100 }
00101 }

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