Sierra Toolkit Version of the Day
CreateAdjacentEntities.cpp
00001 #include <map>
00002 #include <set>
00003 #include <algorithm>
00004 
00005 #include <stk_mesh/base/Types.hpp>
00006 #include <stk_mesh/base/BulkData.hpp>
00007 #include <stk_mesh/base/BulkModification.hpp>
00008 #include <stk_mesh/base/Entity.hpp>
00009 #include <stk_mesh/base/GetBuckets.hpp>
00010 #include <stk_mesh/base/MetaData.hpp>
00011 #include <stk_mesh/base/Selector.hpp>
00012 #include <stk_mesh/base/Relation.hpp>
00013 
00014 #include <stk_mesh/fem/FEMMetaData.hpp>
00015 #include <stk_mesh/fem/FEMHelpers.hpp>
00016 #include <stk_mesh/fem/CellTopology.hpp>
00017 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
00018 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
00019 
00020 #include <stk_util/parallel/ParallelComm.hpp>
00021 
00022 namespace stk {
00023 namespace mesh {
00024 
00025 namespace {
00026 
00027 
00028 struct EntitySubcellComponent {
00029   public:
00030     EntitySubcellComponent()
00031       : entity(NULL)
00032       , subcell_rank(0)
00033       , subcell_id(0)
00034   {}
00035 
00036     EntitySubcellComponent(
00037         Entity     * arg_entity,
00038         EntityRank   arg_subcell_rank,
00039         unsigned     arg_subcell_id
00040         )
00041       : entity(arg_entity)
00042       , subcell_rank(arg_subcell_rank)
00043       , subcell_id(arg_subcell_id)
00044   {}
00045 
00046     Entity     * entity;
00047     EntityRank subcell_rank;
00048     unsigned   subcell_id;
00049 };
00050 
00051 
00052 
00053 void get_entities_with_given_subcell(
00054   const CellTopologyData * subcell_topology,
00055   const EntityRank subcell_rank,
00056   const EntityVector & subcell_nodes,
00057   const EntityRank entities_rank,
00058   std::vector< EntitySubcellComponent> & entities_with_subcell
00059                                      )
00060 {
00061   // Get all entities that have relations to all the subcell nodes
00062   EntityVector entities;
00063   get_entities_through_relations(subcell_nodes,
00064                                  entities_rank,
00065                                  entities);
00066 
00067   // For all such entities, add id info for the subcell if the subcell
00068   // nodes compose a valid subcell of the entity
00069   for (EntityVector::const_iterator eitr = entities.begin();
00070        eitr != entities.end(); ++eitr) {
00071     int local_subcell_num = fem::get_entity_subcell_id(
00072       **eitr,
00073       subcell_rank,
00074       subcell_topology,
00075       subcell_nodes);
00076     if ( local_subcell_num != -1) {
00077       entities_with_subcell.push_back(EntitySubcellComponent(*eitr, subcell_rank, local_subcell_num));
00078     }
00079   }
00080 }
00081 
00082 
00083 
00084 // Check if 3d element topology topo is degenerate
00085 bool is_degenerate( const fem::CellTopology & topo)
00086 {
00087   return topo.getSideCount() < 3;
00088 }
00089 
00090 // Check if entity has a specific relation to an entity of subcell_rank
00091 bool relation_exist( const Entity & entity, EntityRank subcell_rank, RelationIdentifier subcell_id )
00092 {
00093   bool found = false;
00094   PairIterRelation relations = entity.relations(subcell_rank);
00095 
00096   for (; !relations.empty(); ++relations) {
00097     if (relations->identifier() == subcell_id) {
00098       found = true;
00099       break;
00100     }
00101   }
00102 
00103   return found;
00104 }
00105 
00106 
00107 
00108 void internal_count_entities_to_create( BulkData & mesh, std::vector<size_t> & entities_to_request) {
00109 
00110   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
00111   const EntityRank element_rank = fem_meta.element_rank();
00112   const EntityRank side_rank = fem_meta.side_rank();
00113   const EntityRank edge_rank = fem_meta.edge_rank();
00114 
00115   Selector select_owned = fem::FEMMetaData::get(mesh).locally_owned_part();
00116 
00117 
00118   BucketVector element_buckets;
00119 
00120   get_buckets( select_owned, mesh.buckets(element_rank), element_buckets);
00121 
00122 
00123   for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
00124     for (BucketVector::iterator bitr = element_buckets.begin();
00125         bitr != element_buckets.end();
00126         ++bitr)
00127     {
00128       Bucket & b = **bitr;
00129       const fem::CellTopology topo = fem::get_cell_topology(b);
00130 
00131       ThrowErrorMsgIf( is_degenerate(topo),
00132           "stk::mesh::create_adjacent_entities(...) does not yet support degenerate topologies (i.e. shells and beams)");
00133 
00134 
00135       if ( !is_degenerate(topo) ) { // don't loop over shell elements
00136 
00137         for (size_t i = 0; i<b.size(); ++i) {
00138 
00139           Entity & elem = b[i];
00140 
00141           const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
00142 
00143           for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
00144 
00145             if ( ! relation_exist( elem, subcell_rank, subcell_id) ) { //
00146 
00147 
00148               EntityVector subcell_nodes;
00149 
00150               const CellTopologyData * subcell_topology =
00151                 fem::get_subcell_nodes(
00152                     elem,
00153                     subcell_rank,
00154                     subcell_id,
00155                     subcell_nodes
00156                     );
00157 
00158               std::vector<EntitySubcellComponent> adjacent_elements;
00159 
00160               get_entities_with_given_subcell(
00161                   subcell_topology,
00162                   subcell_rank,
00163                   subcell_nodes,
00164                   element_rank,
00165                   adjacent_elements
00166                   );
00167 
00168               std::reverse( subcell_nodes.begin(), subcell_nodes.end());
00169 
00170               get_entities_with_given_subcell(
00171                   subcell_topology,
00172                   subcell_rank,
00173                   subcell_nodes,
00174                   element_rank,
00175                   adjacent_elements
00176                   );
00177 
00178               bool current_elem_has_lowest_id = true;
00179               //does this process own the element with the lowest id?
00180 
00181               for (std::vector<EntitySubcellComponent>::iterator adjacent_itr = adjacent_elements.begin();
00182                   adjacent_itr != adjacent_elements.end();
00183                   ++adjacent_itr)
00184               {
00185                 if (adjacent_itr->entity->identifier() < elem.identifier()) {
00186                   current_elem_has_lowest_id = false;
00187                   break;
00188                 }
00189               }
00190 
00191               // This process owns the lowest element so
00192               // needs to generate a request to create
00193               // the subcell
00194               if (current_elem_has_lowest_id) {
00195                 entities_to_request[subcell_rank]++;
00196               }
00197             }
00198           }
00199         }
00200       }
00201     }
00202   }
00203 }
00204 
00205 void request_entities(
00206    BulkData & mesh,
00207    std::vector<size_t> & entities_to_request,
00208    std::vector< EntityVector > & requested_entities)
00209 {
00210   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
00211   const size_t num_ranks = fem_meta.entity_rank_count();
00212 
00213   requested_entities.clear();
00214   requested_entities.resize(num_ranks);
00215 
00216   EntityVector requested_entities_flat_vector;
00217   mesh.generate_new_entities(entities_to_request, requested_entities_flat_vector);
00218 
00219   EntityVector::iterator b_itr = requested_entities_flat_vector.begin();
00220 
00221   for (size_t i=0; i<num_ranks; ++i) {
00222     EntityVector & temp = requested_entities[i];
00223     temp.insert(temp.begin(), b_itr, b_itr + entities_to_request[i]);
00224     b_itr += entities_to_request[i];
00225   }
00226 
00227   ThrowRequire(b_itr == requested_entities_flat_vector.end());
00228 
00229 }
00230 
00231 void internal_create_adjacent_entities( BulkData & mesh, const PartVector & arg_add_parts, std::vector<size_t> & entities_to_request) {
00232 
00233   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
00234   const EntityRank element_rank = fem_meta.element_rank();
00235   const EntityRank side_rank = fem_meta.side_rank();
00236   const EntityRank edge_rank = fem_meta.edge_rank();
00237 
00238   const size_t num_ranks = fem_meta.entity_rank_count();
00239 
00240   Selector select_owned = fem_meta.locally_owned_part();
00241 
00242   BucketVector element_buckets;
00243 
00244   get_buckets( select_owned, mesh.buckets(element_rank), element_buckets);
00245 
00246 
00247   mesh.modification_begin();
00248 
00249 
00250   std::vector< EntityVector > requested_entities;
00251 
00252   request_entities(
00253       mesh,
00254       entities_to_request,
00255       requested_entities
00256       );
00257 
00258   std::vector<size_t> entities_used(num_ranks, 0);
00259 
00260   for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
00261     for (BucketVector::iterator bitr = element_buckets.begin();
00262         bitr != element_buckets.end();
00263         ++bitr)
00264     {
00265       Bucket & b = **bitr;
00266       const fem::CellTopology topo = fem::get_cell_topology(b);
00267 
00268       if ( !is_degenerate(topo) ) { // don't loop over shell elements
00269 
00270         for (size_t i = 0; i<b.size(); ++i) {
00271 
00272           Entity & elem = b[i];
00273 
00274           const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
00275 
00276           for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
00277 
00278             if ( ! relation_exist( elem, subcell_rank, subcell_id) ) { //
00279 
00280 
00281               EntityVector subcell_nodes;
00282 
00283               const CellTopologyData * subcell_topology =
00284                 fem::get_subcell_nodes(
00285                     elem,
00286                     subcell_rank,
00287                     subcell_id,
00288                     subcell_nodes
00289                     );
00290 
00291               std::vector<EntitySubcellComponent> adjacent_elements;
00292 
00293               get_entities_with_given_subcell(
00294                   subcell_topology,
00295                   subcell_rank,
00296                   subcell_nodes,
00297                   element_rank,
00298                   adjacent_elements
00299                   );
00300 
00301               std::reverse( subcell_nodes.begin(), subcell_nodes.end());
00302 
00303               get_entities_with_given_subcell(
00304                   subcell_topology,
00305                   subcell_rank,
00306                   subcell_nodes,
00307                   element_rank,
00308                   adjacent_elements
00309                   );
00310 
00311               bool current_elem_has_lowest_id = true;
00312               //does this process own the element with the lowest id?
00313 
00314               for (std::vector<EntitySubcellComponent>::iterator adjacent_itr = adjacent_elements.begin();
00315                   adjacent_itr != adjacent_elements.end();
00316                   ++adjacent_itr)
00317               {
00318                 if (adjacent_itr->entity->identifier() < elem.identifier()) {
00319                   current_elem_has_lowest_id = false;
00320                   break;
00321                 }
00322               }
00323 
00324               // This process owns the lowest element so
00325               // needs to generate a request to create
00326               // the subcell
00327               if (current_elem_has_lowest_id) {
00328                 Entity & subcell = * requested_entities[subcell_rank][entities_used[subcell_rank]++];
00329 
00330 
00331                 //declare the node relations for this subcell
00332                 for (size_t n = 0; n<subcell_nodes.size(); ++n) {
00333                   Entity & node = *subcell_nodes[n];
00334                   mesh.declare_relation( subcell, node, n);
00335                 }
00336 
00337                 mesh.declare_relation( elem, subcell, subcell_id);
00338 
00339 
00340                 PartVector empty_remove_parts;
00341 
00342                 PartVector add_parts = arg_add_parts;
00343                 add_parts.push_back( & fem_meta.get_cell_topology_root_part(topo.getCellTopologyData(subcell_rank,subcell_id)));
00344 
00345                 mesh.change_entity_parts(subcell, add_parts, empty_remove_parts);
00346 
00347               }
00348             }
00349           }
00350         }
00351       }
00352     }
00353   }
00354 
00355   mesh.modification_end();
00356 }
00357 
00358 void complete_connectivity( BulkData & mesh ) {
00359 
00360   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
00361   const EntityRank element_rank = fem_meta.element_rank();
00362   const EntityRank side_rank = fem_meta.side_rank();
00363   const EntityRank edge_rank = fem_meta.edge_rank();
00364 
00365   Selector select_owned_or_shared = fem_meta.locally_owned_part() | fem_meta.globally_shared_part();
00366 
00367   BucketVector element_buckets;
00368 
00369   // Add the relationship to the correct entities. So far, we've added the
00370   // edges/side to the sharing element w/ lowest id, but all other elements
00371   // that contain that edge/side still need to have the relationship set up.
00372   // We do that below...
00373 
00374   for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
00375 
00376     mesh.modification_begin();
00377     for (EntityRank entity_rank = element_rank; entity_rank > subcell_rank; --entity_rank) {
00378 
00379 
00380       BucketVector entity_buckets;
00381 
00382 
00383       get_buckets(select_owned_or_shared, mesh.buckets(entity_rank),entity_buckets);
00384 
00385       for (BucketVector::iterator bitr = entity_buckets.begin();
00386           bitr != entity_buckets.end();
00387           ++bitr)
00388       {
00389         Bucket & b = **bitr;
00390         const fem::CellTopology topo = fem::get_cell_topology(b);
00391 
00392         ThrowErrorMsgIf( is_degenerate(topo),
00393           "stk::mesh::create_adjacent_entities(...) does not yet support degenerate topologies (i.e. shells and beams)");
00394 
00395         {
00396           for (size_t i = 0; i<b.size(); ++i) {
00397             Entity & entity = b[i];
00398 
00399             const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
00400 
00401             for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
00402 
00403               if ( !relation_exist(entity, subcell_rank, subcell_id) ) {
00404 
00405                 EntityVector subcell_nodes;
00406 
00407                 const CellTopologyData * subcell_topology =
00408                   fem::get_subcell_nodes(
00409                       entity,
00410                       subcell_rank,
00411                       subcell_id,
00412                       subcell_nodes
00413                       );
00414 
00415                 std::vector<EntitySubcellComponent> adjacent_entities;
00416 
00417                 // add polarity information to newly created relations
00418                 // polarity information is required to correctly attached
00419                 // degenerate elements to the correct faces and edges
00420 
00421                 get_entities_with_given_subcell(
00422                     subcell_topology,
00423                     subcell_rank,
00424                     subcell_nodes,
00425                     subcell_rank,
00426                     adjacent_entities
00427                     );
00428 
00429                 std::reverse( subcell_nodes.begin(), subcell_nodes.end());
00430 
00431                 get_entities_with_given_subcell(
00432                     subcell_topology,
00433                     subcell_rank,
00434                     subcell_nodes,
00435                     subcell_rank,
00436                     adjacent_entities
00437                     );
00438 
00439 
00440                 if ( !adjacent_entities.empty()) {
00441 
00442                   mesh.declare_relation( entity, *adjacent_entities[0].entity, subcell_id);
00443                 }
00444               }
00445             }
00446           }
00447         }
00448       }
00449     }
00450     mesh.modification_end();
00451   }
00452 
00453 }
00454 
00455 } // un-named namespace
00456 
00457 void create_adjacent_entities( BulkData & mesh, PartVector & arg_add_parts)
00458 {
00459   ThrowErrorMsgIf(mesh.synchronized_state() == BulkData::MODIFIABLE,
00460                   "stk::mesh::skin_mesh is not SYNCHRONIZED");
00461 
00462   // to handle degenerate topologies we anticipate the following order of operations
00463   //
00464   // complete_connectivity
00465   // count degenerate entities to create
00466   // create degenerate entities
00467   // complete_connectivity
00468   // count non degenerate entities to create
00469   // create non degenerate entities
00470   // complete_connectivity
00471   //
00472   // to complete the connectivity (with degenerate elements) we require that
00473   // polarity information to be stored on each relation
00474 
00475 
00476   complete_connectivity(mesh);
00477 
00478 
00479   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
00480   const size_t num_ranks = fem_meta.entity_rank_count();
00481   std::vector<size_t> entities_to_request(num_ranks, 0);
00482 
00483   internal_count_entities_to_create( mesh, entities_to_request);
00484 
00485   internal_create_adjacent_entities( mesh, arg_add_parts, entities_to_request);
00486 
00487 
00488   complete_connectivity(mesh);
00489 
00490 
00491 }
00492 
00493 }
00494 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends