Sierra Toolkit Version of the Day
SkinMesh.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/GetEntities.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/BoundaryAnalysis.hpp>
00015 #include <stk_mesh/fem/FEMHelpers.hpp>
00016 #include <stk_mesh/fem/FEMMetaData.hpp>
00017 #include <stk_mesh/fem/SkinMesh.hpp>
00018 
00019 #include <stk_util/parallel/ParallelComm.hpp>
00020 
00021 namespace stk {
00022 namespace mesh {
00023 
00024 namespace {
00025 
00026 typedef std::pair< const CellTopologyData *, EntityVector > SideKey;
00027 typedef std::vector< EntitySideComponent >                  SideVector;
00028 typedef std::map< SideKey, SideVector>                      BoundaryMap;
00029 
00030 //Comparator class to sort EntitySideComponent
00031 //first by entity identifier and then by side_ordinal
00032 class EntitySideComponentLess {
00033   public:
00034     bool operator () (const EntitySideComponent & lhs, const EntitySideComponent &rhs)
00035     {
00036       const EntityId lhs_elem_id = lhs.entity->identifier();
00037       const EntityId rhs_elem_id = rhs.entity->identifier();
00038 
00039       return (lhs_elem_id != rhs_elem_id) ?
00040         (lhs_elem_id < rhs_elem_id) :
00041         (lhs.side_ordinal < rhs.side_ordinal);
00042     }
00043 };
00044 
00045 //Convience class to help with communication.
00046 //sort first by element identifier and then by side_ordinal
00047 class ElementIdSide {
00048   public:
00049     EntityId        elem_id;
00050     RelationIdentifier side_ordinal;
00051 
00052     ElementIdSide() :
00053       elem_id(0), side_ordinal(0) {}
00054 
00055     ElementIdSide( EntityId id, RelationIdentifier ordinal) :
00056       elem_id(id), side_ordinal(ordinal) {}
00057 
00058     ElementIdSide( const EntitySideComponent & esc) :
00059       elem_id(esc.entity->identifier()), side_ordinal(esc.side_ordinal) {}
00060 
00061     ElementIdSide & operator = ( const ElementIdSide & rhs) {
00062       elem_id      = rhs.elem_id;
00063       side_ordinal = rhs.side_ordinal;
00064       return *this;
00065     }
00066 
00067     //compare first by elem_id and then by side_ordinal
00068     bool operator < ( const ElementIdSide & rhs) const {
00069       const ElementIdSide & lhs = *this;
00070 
00071       return (lhs.elem_id != rhs.elem_id) ?
00072         (lhs.elem_id < rhs.elem_id) :
00073         (lhs.side_ordinal < rhs.side_ordinal);
00074     }
00075 };
00076 
00077 //a reverse map used in unpacking to determine what side needs
00078 //to be created
00079 typedef std::map< ElementIdSide, SideKey>   ReverseBoundaryMap;
00080 
00081 //a convience class to help with packing the comm buffer
00082 class SideCommHelper {
00083   public:
00084     unsigned       proc_to;
00085     ElementIdSide  creating_elem_id_side; //used as a look up in the reverse boundary map
00086     EntityId       generated_side_id;
00087 
00088     SideCommHelper() :
00089       proc_to(0), creating_elem_id_side(), generated_side_id(0) {};
00090 
00091     SideCommHelper( unsigned p, const ElementIdSide & creating_eid, const EntityId side_id) :
00092       proc_to(p), creating_elem_id_side(creating_eid), generated_side_id(side_id) {}
00093 
00094     SideCommHelper( const SideCommHelper & sch) :
00095       proc_to(sch.proc_to),
00096       creating_elem_id_side(sch.creating_elem_id_side),
00097       generated_side_id(sch.generated_side_id)
00098   {}
00099 
00100     SideCommHelper & operator = (const SideCommHelper & rhs) {
00101       proc_to                = rhs.proc_to;
00102       creating_elem_id_side  = rhs.creating_elem_id_side;
00103       generated_side_id      = rhs.generated_side_id;
00104       return *this;
00105     }
00106 
00107     //compare first by proc_to, then by creating elem_id
00108     bool operator < ( const SideCommHelper & rhs) const {
00109       const SideCommHelper & lhs = *this;
00110 
00111       return (lhs.proc_to != rhs.proc_to) ?
00112         (lhs.proc_to < rhs.proc_to) :
00113         (lhs.creating_elem_id_side < rhs.creating_elem_id_side);
00114     }
00115 };
00116 
00117 // Use permutation that starts with lowest entity id
00118 void ensure_consistent_order(EntityVector & side_entities)
00119 {
00120   ThrowRequire( !side_entities.empty() );
00121 
00122   EntityId lowest_id = side_entities.front()->identifier();
00123   unsigned idx_of_lowest_id = 0;
00124 
00125   for (unsigned idx = 1; idx < side_entities.size(); ++idx) {
00126     EntityId curr_id = side_entities[idx]->identifier();
00127     if (curr_id < lowest_id) {
00128       idx_of_lowest_id = idx;
00129       lowest_id = curr_id;
00130     }
00131   }
00132 
00133   if (idx_of_lowest_id != 0) {
00134     std::rotate(side_entities.begin(),
00135                 side_entities.begin() + idx_of_lowest_id,
00136                 side_entities.end());
00137   }
00138 }
00139 
00140 // populate the side_map with 'owned' sides that need to be created
00141 //
00142 // a side needs to be created if the outside is NULL and the element
00143 // does not have a current relation to a side for the side_ordinal
00144 void add_owned_sides_to_map(
00145     const BulkData & mesh,
00146     const EntityRank element_rank,
00147     const EntitySideVector & boundary,
00148     BoundaryMap & side_map)
00149 {
00150   for (stk::mesh::EntitySideVector::const_iterator itr = boundary.begin();
00151       itr != boundary.end(); ++itr) {
00152     const EntitySideComponent & inside = itr->inside;
00153     const EntitySideComponent & outside = itr->outside;
00154     const RelationIdentifier side_ordinal = inside.side_ordinal;
00155     const Entity& inside_entity = *(inside.entity);
00156 
00157     if ( inside_entity.owner_rank() == mesh.parallel_rank() &&
00158         outside.entity == NULL ) {
00159       // search through existing sides
00160       PairIterRelation existing_sides = inside_entity.relations(element_rank -1);
00161       for (; existing_sides.first != existing_sides.second &&
00162           existing_sides.first->identifier() != side_ordinal ;
00163           ++existing_sides.first);
00164 
00165       // a relation the side was not found
00166       if (existing_sides.first == existing_sides.second) {
00167         //create the side_key
00168         //side_key.first := CellTopologyData * of the side topology
00169         //side_key.second := EntityVector * of the side_nodes with the nodes in the correct
00170         //                                            permutation for the side starting with
00171         //                                            the node with the smallest identifier
00172         SideKey side_key;
00173 
00174         side_key.first = fem::get_subcell_nodes(
00175             inside_entity,
00176             element_rank - 1, // subcell rank
00177             side_ordinal,     // subcell identifier
00178             side_key.second  // subcell nodes
00179             );
00180         ensure_consistent_order(side_key.second);
00181 
00182         //add this side to the side_map
00183         side_map[side_key].push_back(inside);
00184       }
00185     }
00186   }
00187 }
00188 
00189 // populate the side_map with 'non-owned' sides that need to be created
00190 // who's side_key is already present in the side_map.  This process
00191 // may need to communicate with the process which own these elements
00192 //
00193 // a side needs to be created if the outside is NULL and the element
00194 // does not have a current relation to a side for the side_ordinal
00195 void add_non_owned_sides_to_map(
00196     const BulkData & mesh,
00197     const EntityRank element_rank,
00198     const EntitySideVector & boundary,
00199     BoundaryMap & side_map)
00200 {
00201   for (stk::mesh::EntitySideVector::const_iterator itr = boundary.begin();
00202       itr != boundary.end(); ++itr) {
00203     const EntitySideComponent & inside = itr->inside;
00204     const EntitySideComponent & outside = itr->outside;
00205     const RelationIdentifier side_ordinal = inside.side_ordinal;
00206     const Entity& inside_entity = *(inside.entity);
00207 
00208     // If this process does NOT own the inside and the outside entity does not exist
00209     if ( inside_entity.owner_rank() != mesh.parallel_rank() &&
00210         outside.entity == NULL ) {
00211       // search through existing sides
00212       PairIterRelation existing_sides = inside_entity.relations(element_rank -1);
00213       for (; existing_sides.first != existing_sides.second &&
00214           existing_sides.first->identifier() != side_ordinal ;
00215           ++existing_sides.first);
00216 
00217       // a relation to the side was not found
00218       if (existing_sides.first == existing_sides.second) {
00219         // Get the nodes for the inside entity
00220         SideKey side_key;
00221 
00222         side_key.first = fem::get_subcell_nodes(
00223             inside_entity,
00224             element_rank - 1, // subcell rank
00225             side_ordinal,     // subcell identifier
00226             side_key.second  // subcell nodes
00227             );
00228         ensure_consistent_order(side_key.second);
00229 
00230         //only add the side if the side_key currently exist in the map
00231         if ( side_map.find(side_key) != side_map.end()) {
00232           side_map[side_key].push_back(inside);
00233         }
00234       }
00235     }
00236   }
00237 }
00238 
00239 //sort the SideVector for each side_key.
00240 //the process who owns the element that appears
00241 //first is responsible for generating the identifier
00242 //for the side and communicating the identifier to the
00243 //remaining process
00244 //
00245 //the ElementIdSide of the first side
00246 //is enough to uniquely identify a side.
00247 //the reverse_map is used to find the side_key
00248 //from this ElementIdSide.
00249 //
00250 //return the number of sides this process
00251 //is responsible for creating
00252 size_t determine_creating_processes(
00253     const BulkData & mesh,
00254     BoundaryMap & side_map,
00255     ReverseBoundaryMap & reverse_side_map)
00256 {
00257   int num_sides_to_create = 0;
00258   for (BoundaryMap::iterator i = side_map.begin();
00259       i != side_map.end();
00260       ++i)
00261   {
00262     const SideKey & side_key = i->first;
00263     SideVector & side_vector = i->second;
00264 
00265     //sort the side vectors base on entity identifier and side ordinal
00266     std::sort( side_vector.begin(), side_vector.end(), EntitySideComponentLess() );
00267 
00268     const EntitySideComponent & first_side = *side_vector.begin();
00269 
00270     //does this process create the first side
00271     //if so it needs to create a side_id
00272     if (first_side.entity->owner_rank() == mesh.parallel_rank()) {
00273       ++num_sides_to_create;
00274     }
00275     const ElementIdSide elem_id_side(first_side);
00276 
00277     reverse_side_map[elem_id_side] = side_key;
00278   }
00279 
00280   return num_sides_to_create;
00281 }
00282 
00283 } //end un-named namespace
00284 
00285 void skin_mesh( BulkData & mesh, EntityRank element_rank, Part * skin_part) {
00286   ThrowErrorMsgIf( mesh.synchronized_state() == BulkData::MODIFIABLE,
00287                    "mesh is not SYNCHRONIZED" );
00288 
00289   EntityVector owned_elements;
00290 
00291   // select owned
00292   Selector owned = fem::FEMMetaData::get(mesh).locally_owned_part();
00293   get_selected_entities( owned,
00294                          mesh.buckets(element_rank),
00295                          owned_elements);
00296 
00297   reskin_mesh(mesh, element_rank, owned_elements, skin_part);
00298 }
00299 
00300 void reskin_mesh( BulkData & mesh, EntityRank element_rank, EntityVector & owned_elements, Part * skin_part) {
00301   ThrowErrorMsgIf( mesh.synchronized_state() == BulkData::MODIFIABLE,
00302                    "mesh is not SYNCHRONIZED" );
00303 
00304   EntityVector elements_closure;
00305 
00306   // compute owned closure
00307   find_closure( mesh, owned_elements, elements_closure );
00308 
00309   // compute boundary
00310   EntitySideVector boundary;
00311   boundary_analysis( mesh, elements_closure, element_rank, boundary);
00312 
00313   BoundaryMap side_map;
00314 
00315   add_owned_sides_to_map(mesh, element_rank, boundary, side_map);
00316 
00317   add_non_owned_sides_to_map(mesh, element_rank, boundary, side_map);
00318 
00319   ReverseBoundaryMap reverse_side_map;
00320 
00321   size_t num_sides_to_create = determine_creating_processes(mesh, side_map, reverse_side_map);
00322 
00323   //begin modification
00324   mesh.modification_begin();
00325 
00326   // formulate request ids for the new sides
00327   std::vector<size_t> requests(fem::FEMMetaData::get(mesh).entity_rank_count(), 0);
00328   requests[element_rank -1] = num_sides_to_create;
00329 
00330   // create the new sides
00331   EntityVector requested_sides;
00332   mesh.generate_new_entities(requests, requested_sides);
00333 
00334   //set to aid comm packing
00335   std::set<SideCommHelper> side_comm_helper_set;
00336 
00337 
00338   size_t current_side = 0;
00339   for ( BoundaryMap::iterator map_itr = side_map.begin();
00340         map_itr!= side_map.end();
00341         ++map_itr)
00342   {
00343     const SideKey & side_key    = map_itr->first;
00344     SideVector    & side_vector = map_itr->second;
00345 
00346     // Only generated keys for sides in which this process
00347     // owns the first element in the side vector
00348     const EntitySideComponent & first_side = *(side_vector.begin());
00349     if ( first_side.entity->owner_rank() == mesh.parallel_rank()) {
00350 
00351       //to be used as the key in the reverse boundary map
00352       //a side can be identified in two ways
00353       //  1. vector of nodes in a correct permutation and a side-topology
00354       //  2. element-id side-ordinal for the side
00355       // the reverse boundary map is used to go from (2) to (1).
00356       const ElementIdSide elem_id_side( first_side);
00357 
00358       Entity & side = *(requested_sides[current_side]);
00359       EntityId side_id = side.identifier();
00360 
00361 
00362       PartVector add_parts ;
00363       {
00364         Part * topo_part = &fem::FEMMetaData::get(mesh).get_cell_topology_root_part(side_key.first);
00365         add_parts.push_back( topo_part);
00366         if (skin_part) {
00367           add_parts.push_back(skin_part);
00368         }
00369       }
00370       mesh.change_entity_parts(side, add_parts);
00371 
00372 
00373       //declare the side->node relations
00374       const EntityVector & nodes = side_key.second;
00375 
00376       for (size_t i = 0; i<nodes.size(); ++i) {
00377         Entity & node = *nodes[i];
00378         mesh.declare_relation( side, node, i);
00379       }
00380 
00381       //declare the elem->side relations
00382       for (SideVector::iterator side_itr = side_vector.begin();
00383           side_itr != side_vector.end();
00384           ++side_itr)
00385       {
00386         Entity & elem = *(side_itr->entity);
00387         //only declare relations for owned elements
00388         if (elem.owner_rank() == mesh.parallel_rank()) {
00389           const RelationIdentifier elem_side_ordinal = side_itr->side_ordinal;
00390           mesh.declare_relation( elem, side, elem_side_ordinal);
00391         }
00392         else {
00393           //add this side to the communication set
00394           side_comm_helper_set.insert( SideCommHelper( elem.owner_rank(), elem_id_side, side_id));
00395         }
00396       }
00397       ++current_side;
00398     }
00399   }
00400 
00401   CommAll comm( mesh.parallel() );
00402 
00403   //pack send buffers
00404   for (int allocation_pass=0; allocation_pass<2; ++allocation_pass) {
00405     if (allocation_pass==1) {
00406      comm.allocate_buffers( mesh.parallel_size() /4, 0);
00407     }
00408 
00409     for (std::set<SideCommHelper>::const_iterator i=side_comm_helper_set.begin();
00410          i != side_comm_helper_set.end();
00411          ++i)
00412     {
00413       const SideCommHelper & side_helper = *i;
00414       comm.send_buffer(side_helper.proc_to)
00415         .pack<EntityId>(side_helper.creating_elem_id_side.elem_id)
00416         .pack<RelationIdentifier>(side_helper.creating_elem_id_side.side_ordinal)
00417         .pack<EntityId>(side_helper.generated_side_id);
00418     }
00419   }
00420 
00421   comm.communicate();
00422 
00423   for ( unsigned ip = 0 ; ip < mesh.parallel_size() ; ++ip ) {
00424     CommBuffer & buf = comm.recv_buffer( ip );
00425     while ( buf.remaining() ) {
00426       ElementIdSide creating_elem_id_side;
00427       EntityId      generated_side_id;
00428 
00429       buf.unpack<EntityId>(creating_elem_id_side.elem_id)
00430          .unpack<RelationIdentifier>(creating_elem_id_side.side_ordinal)
00431          .unpack<EntityId>(generated_side_id);
00432 
00433       //find the side in the reverse boundary map
00434       const SideKey & side_key = reverse_side_map[creating_elem_id_side];
00435 
00436       //get the SideVector for the corresponding side_key
00437       SideVector & side_vector = side_map[side_key];
00438 
00439       PartVector add_parts ;
00440       {
00441         Part * topo_part = &fem::FEMMetaData::get(mesh).get_cell_topology_root_part(side_key.first);
00442         add_parts.push_back( topo_part);
00443         if (skin_part) {
00444           add_parts.push_back(skin_part);
00445         }
00446       }
00447       Entity & side = mesh.declare_entity(element_rank-1,generated_side_id,add_parts);
00448 
00449       //declare the side->node relations
00450       const EntityVector & nodes = side_key.second;
00451 
00452       for (size_t i = 0; i<nodes.size(); ++i) {
00453         Entity & node = *nodes[i];
00454         mesh.declare_relation( side, node, i);
00455       }
00456 
00457       //declare the elem->side relations
00458       for (SideVector::iterator side_itr = side_vector.begin();
00459            side_itr != side_vector.end();
00460            ++side_itr)
00461       {
00462         Entity & elem = *(side_itr->entity);
00463         //only declare relations for owned elements
00464         if (elem.owner_rank() == mesh.parallel_rank()) {
00465           const RelationIdentifier elem_side_ordinal = side_itr->side_ordinal;
00466           mesh.declare_relation( elem, side, elem_side_ordinal);
00467         }
00468       }
00469     }
00470   }
00471   mesh.modification_end();
00472 }
00473 
00474 }
00475 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends