Sierra Toolkit  Version of the Day
CreateAdjacentEntities.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/GetBuckets.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/FEMMetaData.hpp>
15 #include <stk_mesh/fem/FEMHelpers.hpp>
16 #include <stk_mesh/fem/CellTopology.hpp>
17 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
18 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
19 
20 #include <stk_util/parallel/ParallelComm.hpp>
21 
22 namespace stk_classic {
23 namespace mesh {
24 
25 namespace {
26 
27 
28 struct EntitySubcellComponent {
29  public:
30  EntitySubcellComponent()
31  : entity(NULL)
32  , subcell_rank(0)
33  , subcell_id(0)
34  {}
35 
36  EntitySubcellComponent(
37  Entity * arg_entity,
38  EntityRank arg_subcell_rank,
39  unsigned arg_subcell_id
40  )
41  : entity(arg_entity)
42  , subcell_rank(arg_subcell_rank)
43  , subcell_id(arg_subcell_id)
44  {}
45 
46  Entity * entity;
47  EntityRank subcell_rank;
48  unsigned subcell_id;
49 };
50 
51 
52 
53 void get_entities_with_given_subcell(
54  const CellTopologyData * subcell_topology,
55  const EntityRank subcell_rank,
56  const EntityVector & subcell_nodes,
57  const EntityRank entities_rank,
58  std::vector< EntitySubcellComponent> & entities_with_subcell
59  )
60 {
61  // Get all entities that have relations to all the subcell nodes
62  EntityVector entities;
63  get_entities_through_relations(subcell_nodes,
64  entities_rank,
65  entities);
66 
67  // For all such entities, add id info for the subcell if the subcell
68  // nodes compose a valid subcell of the entity
69  for (EntityVector::const_iterator eitr = entities.begin();
70  eitr != entities.end(); ++eitr) {
71  int local_subcell_num = fem::get_entity_subcell_id(
72  **eitr,
73  subcell_rank,
74  subcell_topology,
75  subcell_nodes);
76  if ( local_subcell_num != -1) {
77  entities_with_subcell.push_back(EntitySubcellComponent(*eitr, subcell_rank, local_subcell_num));
78  }
79  }
80 }
81 
82 
83 
84 // Check if 3d element topology topo is degenerate
85 bool is_degenerate( const fem::CellTopology & topo)
86 {
87  return topo.getSideCount() < 3;
88 }
89 
90 // Check if entity has a specific relation to an entity of subcell_rank
91 bool relation_exist( const Entity & entity, EntityRank subcell_rank, RelationIdentifier subcell_id )
92 {
93  bool found = false;
94  PairIterRelation relations = entity.relations(subcell_rank);
95 
96  for (; !relations.empty(); ++relations) {
97  if (relations->identifier() == subcell_id) {
98  found = true;
99  break;
100  }
101  }
102 
103  return found;
104 }
105 
106 
107 
108 void internal_count_entities_to_create( BulkData & mesh, std::vector<size_t> & entities_to_request) {
109 
110  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
111  const EntityRank element_rank = fem_meta.element_rank();
112  const EntityRank side_rank = fem_meta.side_rank();
113  const EntityRank edge_rank = fem_meta.edge_rank();
114 
115  Selector select_owned = fem::FEMMetaData::get(mesh).locally_owned_part();
116 
117 
118  BucketVector element_buckets;
119 
120  get_buckets( select_owned, mesh.buckets(element_rank), element_buckets);
121 
122 
123  for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
124  for (BucketVector::iterator bitr = element_buckets.begin();
125  bitr != element_buckets.end();
126  ++bitr)
127  {
128  Bucket & b = **bitr;
129  const fem::CellTopology topo = fem::get_cell_topology(b);
130 
131  ThrowErrorMsgIf( is_degenerate(topo),
132  "stk_classic::mesh::create_adjacent_entities(...) does not yet support degenerate topologies (i.e. shells and beams)");
133 
134 
135  if ( !is_degenerate(topo) ) { // don't loop over shell elements
136 
137  for (size_t i = 0; i<b.size(); ++i) {
138 
139  Entity & elem = b[i];
140 
141  const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
142 
143  for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
144 
145  if ( ! relation_exist( elem, subcell_rank, subcell_id) ) { //
146 
147 
148  EntityVector subcell_nodes;
149 
150  const CellTopologyData * subcell_topology =
152  elem,
153  subcell_rank,
154  subcell_id,
155  subcell_nodes
156  );
157 
158  std::vector<EntitySubcellComponent> adjacent_elements;
159 
160  get_entities_with_given_subcell(
161  subcell_topology,
162  subcell_rank,
163  subcell_nodes,
164  element_rank,
165  adjacent_elements
166  );
167 
168  std::reverse( subcell_nodes.begin(), subcell_nodes.end());
169 
170  get_entities_with_given_subcell(
171  subcell_topology,
172  subcell_rank,
173  subcell_nodes,
174  element_rank,
175  adjacent_elements
176  );
177 
178  bool current_elem_has_lowest_id = true;
179  //does this process own the element with the lowest id?
180 
181  for (std::vector<EntitySubcellComponent>::iterator adjacent_itr = adjacent_elements.begin();
182  adjacent_itr != adjacent_elements.end();
183  ++adjacent_itr)
184  {
185  if (adjacent_itr->entity->identifier() < elem.identifier()) {
186  current_elem_has_lowest_id = false;
187  break;
188  }
189  }
190 
191  // This process owns the lowest element so
192  // needs to generate a request to create
193  // the subcell
194  if (current_elem_has_lowest_id) {
195  entities_to_request[subcell_rank]++;
196  }
197  }
198  }
199  }
200  }
201  }
202  }
203 }
204 
205 void request_entities(
206  BulkData & mesh,
207  std::vector<size_t> & entities_to_request,
208  std::vector< EntityVector > & requested_entities)
209 {
210  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
211  const size_t num_ranks = fem_meta.entity_rank_count();
212 
213  requested_entities.clear();
214  requested_entities.resize(num_ranks);
215 
216  EntityVector requested_entities_flat_vector;
217  mesh.generate_new_entities(entities_to_request, requested_entities_flat_vector);
218 
219  EntityVector::iterator b_itr = requested_entities_flat_vector.begin();
220 
221  for (size_t i=0; i<num_ranks; ++i) {
222  EntityVector & temp = requested_entities[i];
223  temp.insert(temp.begin(), b_itr, b_itr + entities_to_request[i]);
224  b_itr += entities_to_request[i];
225  }
226 
227  ThrowRequire(b_itr == requested_entities_flat_vector.end());
228 
229 }
230 
231 void internal_create_adjacent_entities( BulkData & mesh, const PartVector & arg_add_parts, std::vector<size_t> & entities_to_request) {
232 
233  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
234  const EntityRank element_rank = fem_meta.element_rank();
235  const EntityRank side_rank = fem_meta.side_rank();
236  const EntityRank edge_rank = fem_meta.edge_rank();
237 
238  const size_t num_ranks = fem_meta.entity_rank_count();
239 
240  Selector select_owned = fem_meta.locally_owned_part();
241 
242  BucketVector element_buckets;
243 
244  get_buckets( select_owned, mesh.buckets(element_rank), element_buckets);
245 
246 
247  mesh.modification_begin();
248 
249 
250  std::vector< EntityVector > requested_entities;
251 
252  request_entities(
253  mesh,
254  entities_to_request,
255  requested_entities
256  );
257 
258  std::vector<size_t> entities_used(num_ranks, 0);
259 
260  for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
261  for (BucketVector::iterator bitr = element_buckets.begin();
262  bitr != element_buckets.end();
263  ++bitr)
264  {
265  Bucket & b = **bitr;
266  const fem::CellTopology topo = fem::get_cell_topology(b);
267 
268  if ( !is_degenerate(topo) ) { // don't loop over shell elements
269 
270  for (size_t i = 0; i<b.size(); ++i) {
271 
272  Entity & elem = b[i];
273 
274  const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
275 
276  for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
277 
278  if ( ! relation_exist( elem, subcell_rank, subcell_id) ) { //
279 
280 
281  EntityVector subcell_nodes;
282 
283  const CellTopologyData * subcell_topology =
285  elem,
286  subcell_rank,
287  subcell_id,
288  subcell_nodes
289  );
290 
291  std::vector<EntitySubcellComponent> adjacent_elements;
292 
293  get_entities_with_given_subcell(
294  subcell_topology,
295  subcell_rank,
296  subcell_nodes,
297  element_rank,
298  adjacent_elements
299  );
300 
301  std::reverse( subcell_nodes.begin(), subcell_nodes.end());
302 
303  get_entities_with_given_subcell(
304  subcell_topology,
305  subcell_rank,
306  subcell_nodes,
307  element_rank,
308  adjacent_elements
309  );
310 
311  bool current_elem_has_lowest_id = true;
312  //does this process own the element with the lowest id?
313 
314  for (std::vector<EntitySubcellComponent>::iterator adjacent_itr = adjacent_elements.begin();
315  adjacent_itr != adjacent_elements.end();
316  ++adjacent_itr)
317  {
318  if (adjacent_itr->entity->identifier() < elem.identifier()) {
319  current_elem_has_lowest_id = false;
320  break;
321  }
322  }
323 
324  // This process owns the lowest element so
325  // needs to generate a request to create
326  // the subcell
327  if (current_elem_has_lowest_id) {
328  Entity & subcell = * requested_entities[subcell_rank][entities_used[subcell_rank]++];
329 
330 
331  //declare the node relations for this subcell
332  for (size_t n = 0; n<subcell_nodes.size(); ++n) {
333  Entity & node = *subcell_nodes[n];
334  mesh.declare_relation( subcell, node, n);
335  }
336 
337  mesh.declare_relation( elem, subcell, subcell_id);
338 
339 
340  PartVector empty_remove_parts;
341 
342  PartVector add_parts = arg_add_parts;
343  add_parts.push_back( & fem_meta.get_cell_topology_root_part(topo.getCellTopologyData(subcell_rank,subcell_id)));
344 
345  mesh.change_entity_parts(subcell, add_parts, empty_remove_parts);
346 
347  }
348  }
349  }
350  }
351  }
352  }
353  }
354 
355  mesh.modification_end();
356 }
357 
358 void complete_connectivity( BulkData & mesh ) {
359 
360  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
361  const EntityRank element_rank = fem_meta.element_rank();
362  const EntityRank side_rank = fem_meta.side_rank();
363  const EntityRank edge_rank = fem_meta.edge_rank();
364 
365  Selector select_owned_or_shared = fem_meta.locally_owned_part() | fem_meta.globally_shared_part();
366 
367  BucketVector element_buckets;
368 
369  // Add the relationship to the correct entities. So far, we've added the
370  // edges/side to the sharing element w/ lowest id, but all other elements
371  // that contain that edge/side still need to have the relationship set up.
372  // We do that below...
373 
374  for ( EntityRank subcell_rank = side_rank; subcell_rank >= edge_rank; --subcell_rank) {
375 
376  mesh.modification_begin();
377  for (EntityRank entity_rank = element_rank; entity_rank > subcell_rank; --entity_rank) {
378 
379 
380  BucketVector entity_buckets;
381 
382 
383  get_buckets(select_owned_or_shared, mesh.buckets(entity_rank),entity_buckets);
384 
385  for (BucketVector::iterator bitr = entity_buckets.begin();
386  bitr != entity_buckets.end();
387  ++bitr)
388  {
389  Bucket & b = **bitr;
390  const fem::CellTopology topo = fem::get_cell_topology(b);
391 
392  ThrowErrorMsgIf( is_degenerate(topo),
393  "stk_classic::mesh::create_adjacent_entities(...) does not yet support degenerate topologies (i.e. shells and beams)");
394 
395  {
396  for (size_t i = 0; i<b.size(); ++i) {
397  Entity & entity = b[i];
398 
399  const unsigned subcell_count = topo.getSubcellCount(subcell_rank);
400 
401  for (size_t subcell_id = 0; subcell_id < subcell_count; ++subcell_id ) {
402 
403  if ( !relation_exist(entity, subcell_rank, subcell_id) ) {
404 
405  EntityVector subcell_nodes;
406 
407  const CellTopologyData * subcell_topology =
409  entity,
410  subcell_rank,
411  subcell_id,
412  subcell_nodes
413  );
414 
415  std::vector<EntitySubcellComponent> adjacent_entities;
416 
417  // add polarity information to newly created relations
418  // polarity information is required to correctly attached
419  // degenerate elements to the correct faces and edges
420 
421  get_entities_with_given_subcell(
422  subcell_topology,
423  subcell_rank,
424  subcell_nodes,
425  subcell_rank,
426  adjacent_entities
427  );
428 
429  std::reverse( subcell_nodes.begin(), subcell_nodes.end());
430 
431  get_entities_with_given_subcell(
432  subcell_topology,
433  subcell_rank,
434  subcell_nodes,
435  subcell_rank,
436  adjacent_entities
437  );
438 
439 
440  if ( !adjacent_entities.empty()) {
441 
442  mesh.declare_relation( entity, *adjacent_entities[0].entity, subcell_id);
443  }
444  }
445  }
446  }
447  }
448  }
449  }
450  mesh.modification_end();
451  }
452 
453 }
454 
455 } // un-named namespace
456 
457 void create_adjacent_entities( BulkData & mesh, PartVector & arg_add_parts)
458 {
459  ThrowErrorMsgIf(mesh.synchronized_state() == BulkData::MODIFIABLE,
460  "stk_classic::mesh::skin_mesh is not SYNCHRONIZED");
461 
462  // to handle degenerate topologies we anticipate the following order of operations
463  //
464  // complete_connectivity
465  // count degenerate entities to create
466  // create degenerate entities
467  // complete_connectivity
468  // count non degenerate entities to create
469  // create non degenerate entities
470  // complete_connectivity
471  //
472  // to complete the connectivity (with degenerate elements) we require that
473  // polarity information to be stored on each relation
474 
475 
476  complete_connectivity(mesh);
477 
478 
479  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(mesh);
480  const size_t num_ranks = fem_meta.entity_rank_count();
481  std::vector<size_t> entities_to_request(num_ranks, 0);
482 
483  internal_count_entities_to_create( mesh, entities_to_request);
484 
485  internal_create_adjacent_entities( mesh, arg_add_parts, entities_to_request);
486 
487 
488  complete_connectivity(mesh);
489 
490 
491 }
492 
493 }
494 }
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
int get_entity_subcell_id(const Entity &entity, const EntityRank subcell_rank, const CellTopologyData *side_topology, const EntityVector &side_nodes)
Given an entity and collection of nodes, return the local id of the subcell that contains those nodes...
const CellTopologyData * get_subcell_nodes(const Entity &entity, EntityRank subcell_rank, unsigned subcell_identifier, EntityVector &subcell_nodes)
Definition: FEMHelpers.cpp:239
PairIterRelation relations() const
All Entity relations for which this entity is a member. The relations are ordered from lowest entity-...
Definition: Entity.hpp:161
Sierra Toolkit.
void get_entities_through_relations(const std::vector< Entity *> &entities, std::vector< Entity *> &entities_related)
Query which mesh entities have a relation to all of the input mesh entities.
Definition: Relation.cpp:156
AllSelectedBucketsRange get_buckets(const Selector &selector, const BulkData &mesh)
Definition: GetBuckets.cpp:26
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.
EntityRank entity_rank(const EntityKey &key)
Given an entity key, return an entity type (rank).