Sierra Toolkit  Version of the Day
BoundaryAnalysis.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 #include <vector>
10 #include <set>
11 #include <algorithm>
12 
13 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
14 #include <stk_mesh/fem/FEMMetaData.hpp>
15 #include <stk_mesh/fem/FEMHelpers.hpp>
16 
17 #include <stk_mesh/base/BulkData.hpp>
18 #include <stk_mesh/base/Entity.hpp>
19 #include <stk_mesh/base/Part.hpp>
20 
21 namespace stk_classic {
22 namespace mesh {
23 
24 namespace {
25 
26 const EntityRank NODE_RANK = fem::FEMMetaData::NODE_RANK;
27 
28 void filter_superimposed_entities(const Entity & entity, EntityVector & entities)
29 {
30  // Get the node entities for the nodes that make up the entity, we'll
31  // use this to check for superimposed entities
32  PairIterRelation irel = entity.relations(NODE_RANK);
33  EntityVector entity_nodes;
34  entity_nodes.reserve(irel.size());
35  for ( ; !irel.empty(); ++irel ) {
36  entity_nodes.push_back(irel->entity());
37  }
38  std::sort(entity_nodes.begin(), entity_nodes.end());
39 
40  // Make sure to remove the all superimposed entities from the list
41  unsigned num_nodes_in_orig_entity = entity_nodes.size();
42  EntityVector current_nodes;
43  current_nodes.resize(num_nodes_in_orig_entity);
44  EntityVector::iterator itr = entities.begin();
45  while ( itr != entities.end() ) {
46  Entity * current_entity = *itr;
47  PairIterRelation relations = current_entity->relations(NODE_RANK);
48 
49  if (current_entity == &entity) {
50  // Superimposed with self by definition
51  itr = entities.erase(itr);
52  }
53  else if (relations.size() != num_nodes_in_orig_entity) {
54  // current_entity has a different number of nodes than entity, they
55  // cannot be superimposed
56  ++itr;
57  }
58  else {
59  for (unsigned i = 0; relations.first != relations.second;
60  ++relations.first, ++i ) {
61  current_nodes[i] = relations.first->entity();
62  }
63  std::sort(current_nodes.begin(), current_nodes.end());
64 
65  bool entities_are_superimposed = entity_nodes == current_nodes;
66  if (entities_are_superimposed) {
67  itr = entities.erase(itr);
68  }
69  else {
70  ++itr;
71  }
72  }
73  }
74 }
75 
88 void get_adjacent_entities( const Entity & entity ,
89  EntityRank subcell_rank ,
90  unsigned subcell_identifier ,
91  std::vector< EntitySideComponent> & adjacent_entities)
92 {
93  adjacent_entities.clear();
94 
95  // Get nodes that make up the subcell we're looking at
96  EntityVector subcell_nodes;
97  const CellTopologyData * subcell_topology = fem::get_subcell_nodes(entity,
98  subcell_rank,
99  subcell_identifier,
100  subcell_nodes);
101 
102 
103  // Given the nodes related to the subcell, find all entities
104  // with the same rank that have a relation to all of these nodes
105  EntityVector potentially_adjacent_entities;
106 
107  get_entities_through_relations(subcell_nodes,
108  entity.entity_rank(),
109  potentially_adjacent_entities);
110 
111  // We don't want to include entities that are superimposed with
112  // the input entity
113  filter_superimposed_entities(entity, potentially_adjacent_entities);
114 
115  // Add the local ids, from the POV of the adj entitiy, to the return value.
116  // Reverse the nodes so that the adjacent entity has them in the positive
117  // orientation
118  std::reverse(subcell_nodes.begin(),subcell_nodes.end());
119 
120  for (EntityVector::const_iterator eitr = potentially_adjacent_entities.begin();
121  eitr != potentially_adjacent_entities.end(); ++eitr) {
122  int local_subcell_num = fem::get_entity_subcell_id(**eitr,
123  subcell_rank,
124  subcell_topology,
125  subcell_nodes);
126  if ( local_subcell_num != -1) {
127  adjacent_entities.push_back(EntitySideComponent(*eitr, local_subcell_num));
128  }
129  }
130 }
131 
132 } // unnamed namespace
133 
134 void boundary_analysis(const BulkData& bulk_data,
135  const EntityVector & entities_closure,
136  EntityRank closure_rank,
137  EntitySideVector& boundary)
138 {
139  const Selector locally_used = fem::FEMMetaData::get(bulk_data).locally_owned_part()
141 
142  // find an iterator that points to the last item in the closure that is of a
143  // lower-order than the closure_rank
144  EntityVector::const_iterator itr =
145  std::lower_bound(entities_closure.begin(),
146  entities_closure.end(),
147  EntityKey(closure_rank, 0),
148  EntityLess());
149 
150  // iterate over all the entities in the closure up to the iterator we computed above
151  for ( ; itr != entities_closure.end() && (*itr)->entity_rank() == closure_rank; ++itr) {
152  // some temporaries for clarity
153  std::vector<EntitySideComponent > adjacent_entities;
154  Entity& curr_entity = **itr;
155  const CellTopologyData* celltopology = fem::get_cell_topology(curr_entity).getCellTopologyData();
156  if (celltopology == NULL) {
157  continue;
158  }
159 
160  unsigned subcell_rank = closure_rank - 1;
161  PairIterRelation relations = curr_entity.relations(NODE_RANK);
162 
163  // iterate over the subcells of the current entity
164  for (unsigned nitr = 0; nitr < celltopology->subcell_count[subcell_rank]; ++nitr) {
165  // find the entities (same rank as subcell) adjacent to this subcell
166  unsigned subcell_identifier = nitr;
167  get_adjacent_entities(curr_entity,
168  closure_rank - 1,
169  subcell_identifier,
170  adjacent_entities);
171 
172  // try to figure out if we want to keep ((curr_entity, subcell_identifier),
173  // (adjacent_entity, [k]))
174  // if so, push into boundary
175 
176  // it is a keeper if adjacent entities[k] is not in the entities closure
177  // AND if either curr_entity OR adjacent entities[k] is a member of the
178  // bulk_data.locally_used
179 
180  if (adjacent_entities.empty()) {
181  EntitySide keeper;
182  keeper.inside.entity = &curr_entity;
183  keeper.inside.side_ordinal = subcell_identifier;
184  keeper.outside.entity = NULL;
185  keeper.outside.side_ordinal = 0;
186  boundary.push_back(keeper);
187  continue;
188  }
189 
190  // iterate over adjacent entities (our neighbors)
191  for (std::vector<EntitySideComponent >::const_iterator
192  adj_itr = adjacent_entities.begin();
193  adj_itr != adjacent_entities.end(); ++adj_itr) {
194  // grab a reference to this neighbor for clarity
195  const Entity& neighbor = *(adj_itr->entity);
196 
197  // see if this neighbor is in the closure, if so, not a keeper
198  bool neighbor_is_in_closure =
199  std::binary_search(entities_closure.begin(),
200  entities_closure.end(),
201  neighbor,
202  EntityLess());
203 
204  if (neighbor_is_in_closure) {
205  continue;
206  }
207 
208  // if neighbor or curr_entity is locally-used, add it to keeper
209  if ( locally_used( neighbor.bucket()) || locally_used( curr_entity.bucket() ) ) {
210  EntitySide keeper;
211  keeper.inside.entity = &curr_entity;
212  keeper.inside.side_ordinal = subcell_identifier;
213  keeper.outside = *adj_itr;
214 
215  boundary.push_back(keeper);
216  }
217  }
218  }
219  }
220 }
221 
222 
223 }
224 }
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
Part & globally_shared_part() const
Subset for the problem domain that is shared with another process. Ghost entities are not members of ...
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
EntityRank entity_rank() const
The rank of this entity.
Definition: Entity.hpp:128
static FEMMetaData & get(const MetaData &meta)
Getter for FEMMetaData off of a MetaData object.