Sierra Toolkit  Version of the Day
SkinMesh.cpp
1 #include <map>
2 #include <set>
3 #include <algorithm>
4 
5 #include <stk_mesh/base/Types.hpp>
6 #include <stk_mesh/base/BulkData.hpp>
7 #include <stk_mesh/base/BulkModification.hpp>
8 #include <stk_mesh/base/Entity.hpp>
9 #include <stk_mesh/base/GetEntities.hpp>
10 #include <stk_mesh/base/MetaData.hpp>
11 #include <stk_mesh/base/Selector.hpp>
12 #include <stk_mesh/base/Relation.hpp>
13 
14 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
15 #include <stk_mesh/fem/FEMHelpers.hpp>
16 #include <stk_mesh/fem/FEMMetaData.hpp>
17 #include <stk_mesh/fem/SkinMesh.hpp>
18 
19 #include <stk_util/parallel/ParallelComm.hpp>
20 
21 namespace stk_classic {
22 namespace mesh {
23 
24 namespace {
25 
26 typedef std::pair< const CellTopologyData *, EntityVector > SideKey;
27 typedef std::vector< EntitySideComponent > SideVector;
28 typedef std::map< SideKey, SideVector> BoundaryMap;
29 
30 //Comparator class to sort EntitySideComponent
31 //first by entity identifier and then by side_ordinal
32 class EntitySideComponentLess {
33  public:
34  bool operator () (const EntitySideComponent & lhs, const EntitySideComponent &rhs)
35  {
36  const EntityId lhs_elem_id = lhs.entity->identifier();
37  const EntityId rhs_elem_id = rhs.entity->identifier();
38 
39  return (lhs_elem_id != rhs_elem_id) ?
40  (lhs_elem_id < rhs_elem_id) :
41  (lhs.side_ordinal < rhs.side_ordinal);
42  }
43 };
44 
45 //Convience class to help with communication.
46 //sort first by element identifier and then by side_ordinal
47 class ElementIdSide {
48  public:
49  EntityId elem_id;
50  RelationIdentifier side_ordinal;
51 
52  ElementIdSide() :
53  elem_id(0), side_ordinal(0) {}
54 
55  ElementIdSide( EntityId id, RelationIdentifier ordinal) :
56  elem_id(id), side_ordinal(ordinal) {}
57 
58  ElementIdSide( const EntitySideComponent & esc) :
59  elem_id(esc.entity->identifier()), side_ordinal(esc.side_ordinal) {}
60 
61  ElementIdSide & operator = ( const ElementIdSide & rhs) {
62  elem_id = rhs.elem_id;
63  side_ordinal = rhs.side_ordinal;
64  return *this;
65  }
66 
67  //compare first by elem_id and then by side_ordinal
68  bool operator < ( const ElementIdSide & rhs) const {
69  const ElementIdSide & lhs = *this;
70 
71  return (lhs.elem_id != rhs.elem_id) ?
72  (lhs.elem_id < rhs.elem_id) :
73  (lhs.side_ordinal < rhs.side_ordinal);
74  }
75 };
76 
77 //a reverse map used in unpacking to determine what side needs
78 //to be created
79 typedef std::map< ElementIdSide, SideKey> ReverseBoundaryMap;
80 
81 //a convience class to help with packing the comm buffer
82 class SideCommHelper {
83  public:
84  unsigned proc_to;
85  ElementIdSide creating_elem_id_side; //used as a look up in the reverse boundary map
86  EntityId generated_side_id;
87 
88  SideCommHelper() :
89  proc_to(0), creating_elem_id_side(), generated_side_id(0) {};
90 
91  SideCommHelper( unsigned p, const ElementIdSide & creating_eid, const EntityId side_id) :
92  proc_to(p), creating_elem_id_side(creating_eid), generated_side_id(side_id) {}
93 
94  SideCommHelper( const SideCommHelper & sch) :
95  proc_to(sch.proc_to),
96  creating_elem_id_side(sch.creating_elem_id_side),
97  generated_side_id(sch.generated_side_id)
98  {}
99 
100  SideCommHelper & operator = (const SideCommHelper & rhs) {
101  proc_to = rhs.proc_to;
102  creating_elem_id_side = rhs.creating_elem_id_side;
103  generated_side_id = rhs.generated_side_id;
104  return *this;
105  }
106 
107  //compare first by proc_to, then by creating elem_id
108  bool operator < ( const SideCommHelper & rhs) const {
109  const SideCommHelper & lhs = *this;
110 
111  return (lhs.proc_to != rhs.proc_to) ?
112  (lhs.proc_to < rhs.proc_to) :
113  (lhs.creating_elem_id_side < rhs.creating_elem_id_side);
114  }
115 };
116 
117 // Use permutation that starts with lowest entity id
118 void ensure_consistent_order(EntityVector & side_entities)
119 {
120  ThrowRequire( !side_entities.empty() );
121 
122  EntityId lowest_id = side_entities.front()->identifier();
123  unsigned idx_of_lowest_id = 0;
124 
125  for (unsigned idx = 1; idx < side_entities.size(); ++idx) {
126  EntityId curr_id = side_entities[idx]->identifier();
127  if (curr_id < lowest_id) {
128  idx_of_lowest_id = idx;
129  lowest_id = curr_id;
130  }
131  }
132 
133  if (idx_of_lowest_id != 0) {
134  std::rotate(side_entities.begin(),
135  side_entities.begin() + idx_of_lowest_id,
136  side_entities.end());
137  }
138 }
139 
140 // populate the side_map with 'owned' sides that need to be created
141 //
142 // a side needs to be created if the outside is NULL and the element
143 // does not have a current relation to a side for the side_ordinal
144 void add_owned_sides_to_map(
145  const BulkData & mesh,
146  const EntityRank element_rank,
147  const EntitySideVector & boundary,
148  BoundaryMap & side_map)
149 {
150  for (stk_classic::mesh::EntitySideVector::const_iterator itr = boundary.begin();
151  itr != boundary.end(); ++itr) {
152  const EntitySideComponent & inside = itr->inside;
153  const EntitySideComponent & outside = itr->outside;
154  const RelationIdentifier side_ordinal = inside.side_ordinal;
155  const Entity& inside_entity = *(inside.entity);
156 
157  if ( inside_entity.owner_rank() == mesh.parallel_rank() &&
158  outside.entity == NULL ) {
159  // search through existing sides
160  PairIterRelation existing_sides = inside_entity.relations(element_rank -1);
161  for (; existing_sides.first != existing_sides.second &&
162  existing_sides.first->identifier() != side_ordinal ;
163  ++existing_sides.first);
164 
165  // a relation the side was not found
166  if (existing_sides.first == existing_sides.second) {
167  //create the side_key
168  //side_key.first := CellTopologyData * of the side topology
169  //side_key.second := EntityVector * of the side_nodes with the nodes in the correct
170  // permutation for the side starting with
171  // the node with the smallest identifier
172  SideKey side_key;
173 
174  side_key.first = fem::get_subcell_nodes(
175  inside_entity,
176  element_rank - 1, // subcell rank
177  side_ordinal, // subcell identifier
178  side_key.second // subcell nodes
179  );
180  ensure_consistent_order(side_key.second);
181 
182  //add this side to the side_map
183  side_map[side_key].push_back(inside);
184  }
185  }
186  }
187 }
188 
189 // populate the side_map with 'non-owned' sides that need to be created
190 // who's side_key is already present in the side_map. This process
191 // may need to communicate with the process which own these elements
192 //
193 // a side needs to be created if the outside is NULL and the element
194 // does not have a current relation to a side for the side_ordinal
195 void add_non_owned_sides_to_map(
196  const BulkData & mesh,
197  const EntityRank element_rank,
198  const EntitySideVector & boundary,
199  BoundaryMap & side_map)
200 {
201  for (stk_classic::mesh::EntitySideVector::const_iterator itr = boundary.begin();
202  itr != boundary.end(); ++itr) {
203  const EntitySideComponent & inside = itr->inside;
204  const EntitySideComponent & outside = itr->outside;
205  const RelationIdentifier side_ordinal = inside.side_ordinal;
206  const Entity& inside_entity = *(inside.entity);
207 
208  // If this process does NOT own the inside and the outside entity does not exist
209  if ( inside_entity.owner_rank() != mesh.parallel_rank() &&
210  outside.entity == NULL ) {
211  // search through existing sides
212  PairIterRelation existing_sides = inside_entity.relations(element_rank -1);
213  for (; existing_sides.first != existing_sides.second &&
214  existing_sides.first->identifier() != side_ordinal ;
215  ++existing_sides.first);
216 
217  // a relation to the side was not found
218  if (existing_sides.first == existing_sides.second) {
219  // Get the nodes for the inside entity
220  SideKey side_key;
221 
222  side_key.first = fem::get_subcell_nodes(
223  inside_entity,
224  element_rank - 1, // subcell rank
225  side_ordinal, // subcell identifier
226  side_key.second // subcell nodes
227  );
228  ensure_consistent_order(side_key.second);
229 
230  //only add the side if the side_key currently exist in the map
231  if ( side_map.find(side_key) != side_map.end()) {
232  side_map[side_key].push_back(inside);
233  }
234  }
235  }
236  }
237 }
238 
239 //sort the SideVector for each side_key.
240 //the process who owns the element that appears
241 //first is responsible for generating the identifier
242 //for the side and communicating the identifier to the
243 //remaining process
244 //
245 //the ElementIdSide of the first side
246 //is enough to uniquely identify a side.
247 //the reverse_map is used to find the side_key
248 //from this ElementIdSide.
249 //
250 //return the number of sides this process
251 //is responsible for creating
252 size_t determine_creating_processes(
253  const BulkData & mesh,
254  BoundaryMap & side_map,
255  ReverseBoundaryMap & reverse_side_map)
256 {
257  int num_sides_to_create = 0;
258  for (BoundaryMap::iterator i = side_map.begin();
259  i != side_map.end();
260  ++i)
261  {
262  const SideKey & side_key = i->first;
263  SideVector & side_vector = i->second;
264 
265  //sort the side vectors base on entity identifier and side ordinal
266  std::sort( side_vector.begin(), side_vector.end(), EntitySideComponentLess() );
267 
268  const EntitySideComponent & first_side = *side_vector.begin();
269 
270  //does this process create the first side
271  //if so it needs to create a side_id
272  if (first_side.entity->owner_rank() == mesh.parallel_rank()) {
273  ++num_sides_to_create;
274  }
275  const ElementIdSide elem_id_side(first_side);
276 
277  reverse_side_map[elem_id_side] = side_key;
278  }
279 
280  return num_sides_to_create;
281 }
282 
283 } //end un-named namespace
284 
285 void skin_mesh( BulkData & mesh, EntityRank element_rank, Part * skin_part) {
286  ThrowErrorMsgIf( mesh.synchronized_state() == BulkData::MODIFIABLE,
287  "mesh is not SYNCHRONIZED" );
288 
289  EntityVector owned_elements;
290 
291  // select owned
292  Selector owned = fem::FEMMetaData::get(mesh).locally_owned_part();
293  get_selected_entities( owned,
294  mesh.buckets(element_rank),
295  owned_elements);
296 
297  reskin_mesh(mesh, element_rank, owned_elements, skin_part);
298 }
299 
300 void reskin_mesh( BulkData & mesh, EntityRank element_rank, EntityVector & owned_elements, Part * skin_part) {
301  ThrowErrorMsgIf( mesh.synchronized_state() == BulkData::MODIFIABLE,
302  "mesh is not SYNCHRONIZED" );
303 
304  EntityVector elements_closure;
305 
306  // compute owned closure
307  find_closure( mesh, owned_elements, elements_closure );
308 
309  // compute boundary
310  EntitySideVector boundary;
311  boundary_analysis( mesh, elements_closure, element_rank, boundary);
312 
313  BoundaryMap side_map;
314 
315  add_owned_sides_to_map(mesh, element_rank, boundary, side_map);
316 
317  add_non_owned_sides_to_map(mesh, element_rank, boundary, side_map);
318 
319  ReverseBoundaryMap reverse_side_map;
320 
321  size_t num_sides_to_create = determine_creating_processes(mesh, side_map, reverse_side_map);
322 
323  //begin modification
324  mesh.modification_begin();
325 
326  // formulate request ids for the new sides
327  std::vector<size_t> requests(fem::FEMMetaData::get(mesh).entity_rank_count(), 0);
328  requests[element_rank -1] = num_sides_to_create;
329 
330  // create the new sides
331  EntityVector requested_sides;
332  mesh.generate_new_entities(requests, requested_sides);
333 
334  //set to aid comm packing
335  std::set<SideCommHelper> side_comm_helper_set;
336 
337 
338  size_t current_side = 0;
339  for ( BoundaryMap::iterator map_itr = side_map.begin();
340  map_itr!= side_map.end();
341  ++map_itr)
342  {
343  const SideKey & side_key = map_itr->first;
344  SideVector & side_vector = map_itr->second;
345 
346  // Only generated keys for sides in which this process
347  // owns the first element in the side vector
348  const EntitySideComponent & first_side = *(side_vector.begin());
349  if ( first_side.entity->owner_rank() == mesh.parallel_rank()) {
350 
351  //to be used as the key in the reverse boundary map
352  //a side can be identified in two ways
353  // 1. vector of nodes in a correct permutation and a side-topology
354  // 2. element-id side-ordinal for the side
355  // the reverse boundary map is used to go from (2) to (1).
356  const ElementIdSide elem_id_side( first_side);
357 
358  Entity & side = *(requested_sides[current_side]);
359  EntityId side_id = side.identifier();
360 
361 
362  PartVector add_parts ;
363  {
364  fem::FEMMetaData &fem_meta_data = fem::FEMMetaData::get(mesh);
365  Part * topo_part = &fem_meta_data.get_cell_topology_root_part(side_key.first);
366  add_parts.push_back( topo_part);
367  if (skin_part) {
368  add_parts.push_back(skin_part);
369  fem::CellTopology topo = fem_meta_data.get_cell_topology(*topo_part);
370  const PartVector topo_parts = skin_part->subsets();
371  for (stk_classic::mesh::PartVector::const_iterator i=topo_parts.begin(); i!=topo_parts.end(); ++i) {
372  // Decode side to element topology. Specific to tri and quad on linear hex and wedges
373  if (topo == fem_meta_data.get_cell_topology(**i)) {
374  if (std::string(side_key.first->name) == "Quadrilateral_4") { // A enum would be nice
375  // Quad could be the face of a hex or wedge.
376  if (std::string("Hexahedron_8") == stk_classic::mesh::fem::get_cell_topology(*side_vector.front().entity).getName() &&
377  std::string::npos != (*i)->name().find("hex8")) { // Magic string
378  add_parts.push_back(*i);
379  }
380  else if (std::string("Wedge_6") == stk_classic::mesh::fem::get_cell_topology(*side_vector.front().entity).getName() &&
381  std::string::npos != (*i)->name().find("wedge6")) { // Magic string
382  add_parts.push_back(*i);
383  }
384  }
385  else {
386  add_parts.push_back(*i);
387  }
388  }
389  }
390  }
391  }
392  mesh.change_entity_parts(side, add_parts);
393 
394 
395  //declare the side->node relations
396  const EntityVector & nodes = side_key.second;
397 
398  for (size_t i = 0; i<nodes.size(); ++i) {
399  Entity & node = *nodes[i];
400  mesh.declare_relation( side, node, i);
401  }
402 
403  //declare the elem->side relations
404  for (SideVector::iterator side_itr = side_vector.begin();
405  side_itr != side_vector.end();
406  ++side_itr)
407  {
408  Entity & elem = *(side_itr->entity);
409  //only declare relations for owned elements
410  if (elem.owner_rank() == mesh.parallel_rank()) {
411  const RelationIdentifier elem_side_ordinal = side_itr->side_ordinal;
412  mesh.declare_relation( elem, side, elem_side_ordinal);
413  }
414  else {
415  //add this side to the communication set
416  side_comm_helper_set.insert( SideCommHelper( elem.owner_rank(), elem_id_side, side_id));
417  }
418  }
419  ++current_side;
420  }
421  }
422 
423  CommAll comm( mesh.parallel() );
424 
425  //pack send buffers
426  for (int allocation_pass=0; allocation_pass<2; ++allocation_pass) {
427  if (allocation_pass==1) {
428  comm.allocate_buffers( mesh.parallel_size() /4, 0);
429  }
430 
431  for (std::set<SideCommHelper>::const_iterator i=side_comm_helper_set.begin();
432  i != side_comm_helper_set.end();
433  ++i)
434  {
435  const SideCommHelper & side_helper = *i;
436  comm.send_buffer(side_helper.proc_to)
437  .pack<EntityId>(side_helper.creating_elem_id_side.elem_id)
438  .pack<RelationIdentifier>(side_helper.creating_elem_id_side.side_ordinal)
439  .pack<EntityId>(side_helper.generated_side_id);
440  }
441  }
442 
443  comm.communicate();
444 
445  for ( unsigned ip = 0 ; ip < mesh.parallel_size() ; ++ip ) {
446  CommBuffer & buf = comm.recv_buffer( ip );
447  while ( buf.remaining() ) {
448  ElementIdSide creating_elem_id_side;
449  EntityId generated_side_id;
450 
451  buf.unpack<EntityId>(creating_elem_id_side.elem_id)
452  .unpack<RelationIdentifier>(creating_elem_id_side.side_ordinal)
453  .unpack<EntityId>(generated_side_id);
454 
455  //find the side in the reverse boundary map
456  const SideKey & side_key = reverse_side_map[creating_elem_id_side];
457 
458  //get the SideVector for the corresponding side_key
459  SideVector & side_vector = side_map[side_key];
460 
461  PartVector add_parts ;
462  {
463  fem::FEMMetaData &fem_meta_data = fem::FEMMetaData::get(mesh);
464  Part * topo_part = &fem_meta_data.get_cell_topology_root_part(side_key.first);
465  add_parts.push_back( topo_part);
466  if (skin_part) {
467  add_parts.push_back(skin_part);
468  fem::CellTopology topo = fem_meta_data.get_cell_topology(*topo_part);
469  PartVector topo_parts = skin_part->subsets();
470  for (stk_classic::mesh::PartVector::const_iterator i=topo_parts.begin(); i!=topo_parts.end(); ++i) {
471  // Decode side to element topology. Specific to tri and quad on linear hex and wedges
472  if (topo == fem_meta_data.get_cell_topology(**i)) {
473  if (std::string(side_key.first->name) == "Quadrilateral_4") { // A enum would be nice
474  // Quad could be the face of a hex or wedge.
475  if (std::string("Hexahedron_8") == stk_classic::mesh::fem::get_cell_topology(*side_vector.front().entity).getName() &&
476  std::string::npos != (*i)->name().find("hex8")) { // Magic string
477  add_parts.push_back(*i);
478  }
479  else if (std::string("Wedge_6") == stk_classic::mesh::fem::get_cell_topology(*side_vector.front().entity).getName() &&
480  std::string::npos != (*i)->name().find("wedge6")) { // Magic string
481  add_parts.push_back(*i);
482  }
483  }
484  else {
485  add_parts.push_back(*i);
486  }
487  }
488  }
489  }
490  }
491  Entity & side = mesh.declare_entity(element_rank-1,generated_side_id,add_parts);
492 
493  //declare the side->node relations
494  const EntityVector & nodes = side_key.second;
495 
496  for (size_t i = 0; i<nodes.size(); ++i) {
497  Entity & node = *nodes[i];
498  mesh.declare_relation( side, node, i);
499  }
500 
501  //declare the elem->side relations
502  for (SideVector::iterator side_itr = side_vector.begin();
503  side_itr != side_vector.end();
504  ++side_itr)
505  {
506  Entity & elem = *(side_itr->entity);
507  //only declare relations for owned elements
508  if (elem.owner_rank() == mesh.parallel_rank()) {
509  const RelationIdentifier elem_side_ordinal = side_itr->side_ordinal;
510  mesh.declare_relation( elem, side, elem_side_ordinal);
511  }
512  }
513  }
514  }
515  mesh.modification_end();
516 }
517 
518 }
519 }
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
void get_selected_entities(const Selector &selector, const std::vector< Bucket * > &input_buckets, std::vector< Entity * > &entities)
Get entities in selected buckets (selected by the given selector instance), and sorted by ID...
Definition: GetEntities.cpp:77
const CellTopologyData * get_subcell_nodes(const Entity &entity, EntityRank subcell_rank, unsigned subcell_identifier, EntityVector &subcell_nodes)
Definition: FEMHelpers.cpp:239
Sierra Toolkit.
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
static FEMMetaData & get(const MetaData &meta)
Getter for FEMMetaData off of a MetaData object.