Sierra Toolkit  Version of the Day
FEMHelpers.cpp
1 #include <stk_mesh/fem/FEMHelpers.hpp>
2 #include <stk_mesh/fem/FEMMetaData.hpp>
3 
4 #include <Shards_CellTopologyTraits.hpp>
5 
6 #include <stk_mesh/base/Types.hpp>
7 #include <stk_mesh/base/BulkData.hpp>
8 
9 #include <stk_util/parallel/ParallelReduce.hpp>
10 
11 #include <sstream>
12 #include <stdexcept>
13 
14 namespace stk_classic {
15 namespace mesh {
16 namespace fem {
17 
18 namespace {
19 
20 void verify_declare_element_side(
21  const BulkData & mesh,
22  const Entity & elem,
23  const unsigned local_side_id
24  )
25 {
26  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
27 
28  const CellTopologyData * const side_top =
29  ( elem_top && local_side_id < elem_top->side_count )
30  ? elem_top->side[ local_side_id ].topology : NULL ;
31 
32  ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
33  "For elem " << print_entity_key(elem) <<
34  ", Bulkdata for 'elem' and mesh are different");
35 
36  ThrowErrorMsgIf( elem_top && local_side_id >= elem_top->side_count,
37  "For elem " << print_entity_key(elem) << ", local_side_id " << local_side_id << ", " <<
38  "local_side_id exceeds " << elem_top->name << ".side_count = " << elem_top->side_count );
39 
40  ThrowErrorMsgIf( side_top == NULL,
41  "For elem " << print_entity_key(elem) << ", local_side_id " << local_side_id << ", " <<
42  "No element topology found");
43 }
44 
45 void verify_declare_element_edge(
46  const BulkData & mesh,
47  const Entity & elem,
48  const unsigned local_edge_id
49  )
50 {
51  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
52 
53  const CellTopologyData * const edge_top =
54  ( elem_top && local_edge_id < elem_top->edge_count )
55  ? elem_top->edge[ local_edge_id ].topology : NULL ;
56 
57  ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
58  "For elem " << print_entity_key(elem) <<
59  ", Bulkdata for 'elem' and mesh are different");
60 
61  ThrowErrorMsgIf( elem_top && local_edge_id >= elem_top->edge_count,
62  "For elem " << print_entity_key(elem) << ", local_edge_id " << local_edge_id << ", " <<
63  "local_edge_id exceeds " << elem_top->name << ".edge_count = " << elem_top->edge_count );
64 
65  ThrowErrorMsgIf( edge_top == NULL,
66  "For elem " << print_entity_key(elem) << ", local_edge_id " << local_edge_id << ", " <<
67  "No element topology found");
68 }
69 
70 } // unnamed namespace
71 
73  Part & part ,
74  const EntityId elem_id ,
75  const EntityId node_id[] )
76 {
77  FEMMetaData & fem_meta = FEMMetaData::get(mesh);
78  const CellTopologyData * const top = fem_meta.get_cell_topology( part ).getCellTopologyData();
79 
80  ThrowErrorMsgIf(top == NULL,
81  "Part " << part.name() << " does not have a local topology");
82 
83  PartVector empty ;
84  PartVector add( 1 ); add[0] = & part ;
85 
86  const EntityRank entity_rank = fem_meta.element_rank();
87 
88  Entity & elem = mesh.declare_entity( entity_rank, elem_id, add );
89 
90  const EntityRank node_rank = fem_meta.node_rank();
91 
92  for ( unsigned i = 0 ; i < top->node_count ; ++i ) {
93  //declare node if it doesn't already exist
94  Entity * node = mesh.get_entity( node_rank , node_id[i]);
95  if ( NULL == node) {
96  node = & mesh.declare_entity( node_rank , node_id[i], empty );
97  }
98 
99  mesh.declare_relation( elem , *node , i );
100  }
101  return elem ;
102 }
103 
105  Entity & elem ,
106  Entity & side,
107  const unsigned local_side_id ,
108  Part * part )
109 {
110  BulkData & mesh = BulkData::get(side);
111 
112  verify_declare_element_side(mesh, elem, local_side_id);
113 
114  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
115 
116  ThrowErrorMsgIf( elem_top == NULL,
117  "Element[" << elem.identifier() << "] has no defined topology" );
118 
119  const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
120 
121  ThrowErrorMsgIf( side_top == NULL,
122  "Element[" << elem.identifier() << "], local_side_id = " <<
123  local_side_id << ", side has no defined topology" );
124 
125  const unsigned * const side_node_map = elem_top->side[ local_side_id ].node ;
126 
127  PartVector add_parts ;
128 
129  if ( part ) { add_parts.push_back( part ); }
130 
131  mesh.change_entity_parts(side, add_parts);
132 
133  mesh.declare_relation( elem , side , local_side_id );
134 
135  PairIterRelation rel = elem.relations( FEMMetaData::NODE_RANK );
136 
137  for ( unsigned i = 0 ; i < side_top->node_count ; ++i ) {
138  Entity & node = * rel[ side_node_map[i] ].entity();
139  mesh.declare_relation( side , node , i );
140  }
141 
142  return side ;
143 }
144 
146  Entity & elem ,
147  Entity & edge,
148  const unsigned local_edge_id ,
149  Part * part )
150 {
151  BulkData & mesh = BulkData::get(edge);
152 
153  // verify_declare_element_edge(mesh, elem, local_edge_id);
154 
155  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
156 
157  ThrowErrorMsgIf( elem_top == NULL,
158  "Element[" << elem.identifier() << "] has no defined topology" );
159 
160  const CellTopologyData * const edge_top = elem_top->edge[ local_edge_id ].topology;
161 
162  ThrowErrorMsgIf( edge_top == NULL,
163  "Element[" << elem.identifier() << "], local_edge_id = " <<
164  local_edge_id << ", edge has no defined topology" );
165 
166  const unsigned * const edge_node_map = elem_top->edge[ local_edge_id ].node ;
167 
168  PartVector add_parts ;
169 
170  if ( part ) { add_parts.push_back( part ); }
171 
172  mesh.change_entity_parts(edge, add_parts);
173 
174  mesh.declare_relation( elem , edge , local_edge_id );
175 
176  PairIterRelation rel = elem.relations( FEMMetaData::NODE_RANK );
177 
178  for ( unsigned i = 0 ; i < edge_top->node_count ; ++i ) {
179  Entity & node = * rel[ edge_node_map[i] ].entity();
180  mesh.declare_relation( edge , node , i );
181  }
182 
183  return edge ;
184 }
185 
187  BulkData & mesh ,
188  const stk_classic::mesh::EntityId global_side_id ,
189  Entity & elem ,
190  const unsigned local_side_id ,
191  Part * part )
192 {
193  verify_declare_element_side(mesh, elem, local_side_id);
194 
195  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
196 
197  ThrowErrorMsgIf( elem_top == NULL,
198  "Element[" << elem.identifier() << "] has no defined topology");
199 
200  const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
201 
202  ThrowErrorMsgIf( side_top == NULL,
203  "Element[" << elem.identifier() << "], local_side_id = " <<
204  local_side_id << ", side has no defined topology" );
205 
206  PartVector empty_parts ;
207  Entity & side = mesh.declare_entity( side_top->dimension , global_side_id, empty_parts );
208  return declare_element_side( elem, side, local_side_id, part);
209 }
210 
212  BulkData & mesh ,
213  const stk_classic::mesh::EntityId global_edge_id ,
214  Entity & elem ,
215  const unsigned local_edge_id ,
216  Part * part )
217 {
218  verify_declare_element_edge(mesh, elem, local_edge_id);
219 
220  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
221 
222  ThrowErrorMsgIf( elem_top == NULL,
223  "Element[" << elem.identifier() << "] has no defined topology");
224 
225 
226  const CellTopologyData * const edge_top = elem_top->edge[ local_edge_id ].topology;
227 
228  ThrowErrorMsgIf( edge_top == NULL,
229  "Element[" << elem.identifier() << "], local_edge_id = " <<
230  local_edge_id << ", edge has no defined topology" );
231 
232  PartVector empty_parts ;
233  Entity & edge = mesh.declare_entity( edge_top->dimension , global_edge_id, empty_parts );
234  return declare_element_edge( elem, edge, local_edge_id, part);
235 }
236 
237 
238 
239 const CellTopologyData * get_subcell_nodes(const Entity & entity ,
240  EntityRank subcell_rank ,
241  unsigned subcell_identifier ,
242  EntityVector & subcell_nodes)
243 {
244  subcell_nodes.clear();
245 
246  // get cell topology
247  const CellTopologyData* celltopology = get_cell_topology(entity).getCellTopologyData();
248 
249  //error checking
250  {
251  //no celltopology defined
252  if (celltopology == NULL) {
253  return NULL;
254  }
255 
256  // valid ranks fall within the dimension of the cell topology
257  const bool bad_rank = subcell_rank >= celltopology->dimension;
258  ThrowInvalidArgMsgIf( bad_rank, "subcell_rank is >= celltopology dimension\n");
259 
260  // subcell_identifier must be less than the subcell count
261  const bool bad_id = subcell_identifier >= celltopology->subcell_count[subcell_rank];
262  ThrowInvalidArgMsgIf( bad_id, "subcell_id is >= subcell_count\n");
263  }
264 
265  // Get the cell topology of the subcell
266  const CellTopologyData * subcell_topology =
267  celltopology->subcell[subcell_rank][subcell_identifier].topology;
268 
269  const int num_nodes_in_subcell = subcell_topology->node_count;
270 
271  // For the subcell, get it's local nodes ids
272  const unsigned* subcell_node_local_ids =
273  celltopology->subcell[subcell_rank][subcell_identifier].node;
274 
275  FEMMetaData & fem_meta = FEMMetaData::get(entity);
276  const EntityRank node_rank = fem_meta.node_rank();
277  PairIterRelation node_relations = entity.relations(node_rank);
278 
279  subcell_nodes.reserve(num_nodes_in_subcell);
280 
281  for (int i = 0; i < num_nodes_in_subcell; ++i ) {
282  subcell_nodes.push_back( node_relations[subcell_node_local_ids[i]].entity() );
283  }
284 
285  return subcell_topology;
286 }
287 
288 
289 int get_entity_subcell_id( const Entity & entity ,
290  const EntityRank subcell_rank,
291  const CellTopologyData * subcell_topology,
292  const std::vector<Entity*>& subcell_nodes )
293 {
294  const int INVALID_SIDE = -1;
295 
296  unsigned num_nodes = subcell_topology->node_count;
297 
298  if (num_nodes != subcell_nodes.size()) {
299  return INVALID_SIDE;
300  }
301 
302  // get topology of elem
303  const CellTopologyData* entity_topology = get_cell_topology(entity).getCellTopologyData();
304  if (entity_topology == NULL) {
305  return INVALID_SIDE;
306  }
307 
308  // get nodal relations for entity
309  FEMMetaData & fem_meta = FEMMetaData::get(entity);
310  const EntityRank node_rank = fem_meta.node_rank();
311  PairIterRelation relations = entity.relations(node_rank);
312 
313  const int num_permutations = subcell_topology->permutation_count;
314 
315  // Iterate over the subcells of entity...
316  for (unsigned local_subcell_ordinal = 0;
317  local_subcell_ordinal < entity_topology->subcell_count[subcell_rank];
318  ++local_subcell_ordinal) {
319 
320  // get topological data for this subcell
321  const CellTopologyData* curr_subcell_topology =
322  entity_topology->subcell[subcell_rank][local_subcell_ordinal].topology;
323 
324  // If topologies are not the same, there is no way the subcells are the same
325  if (subcell_topology == curr_subcell_topology) {
326 
327  const unsigned* const subcell_node_map = entity_topology->subcell[subcell_rank][local_subcell_ordinal].node;
328 
329  // Taking all positive permutations into account, check if this subcell
330  // has the same nodes as the subcell_nodes argument. Note that this
331  // implementation preserves the node-order so that we can take
332  // entity-orientation into account.
333  for (int p = 0; p < num_permutations; ++p) {
334 
335  if (curr_subcell_topology->permutation[p].polarity ==
336  CELL_PERMUTATION_POLARITY_POSITIVE) {
337 
338  const unsigned * const perm_node =
339  curr_subcell_topology->permutation[p].node ;
340 
341  bool all_match = true;
342  for (unsigned j = 0 ; j < num_nodes; ++j ) {
343  if (subcell_nodes[j] !=
344  relations[subcell_node_map[perm_node[j]]].entity()) {
345  all_match = false;
346  break;
347  }
348  }
349 
350  // all nodes were the same, we have a match
351  if ( all_match ) {
352  return local_subcell_ordinal ;
353  }
354  }
355  }
356  }
357  }
358 
359  return INVALID_SIDE;
360 }
361 
363  std::vector<size_t> & counts ,
364  bool local_flag )
365 {
366  const size_t zero = 0 ;
367 
368  // Count locally owned entities
369 
370  const FEMMetaData & S = FEMMetaData::get(M);
371  const unsigned entity_rank_count = S.entity_rank_count();
372  const size_t comm_count = entity_rank_count + 1 ;
373 
374  std::vector<size_t> local( comm_count , zero );
375  std::vector<size_t> global( comm_count , zero );
376 
377  ParallelMachine comm = M.parallel();
378  Part & owns = S.locally_owned_part();
379 
380  for ( unsigned i = 0 ; i < entity_rank_count ; ++i ) {
381  const std::vector<Bucket*> & ks = M.buckets( i );
382 
383  std::vector<Bucket*>::const_iterator ik ;
384 
385  for ( ik = ks.begin() ; ik != ks.end() ; ++ik ) {
386  if ( has_superset( **ik , owns ) ) {
387  local[i] += (*ik)->size();
388  }
389  }
390  }
391 
392  local[ entity_rank_count ] = local_flag ;
393 
394  stk_classic::all_reduce_sum( comm , & local[0] , & global[0] , comm_count );
395 
396  counts.assign( global.begin() , global.begin() + entity_rank_count );
397 
398  return 0 < global[ entity_rank_count ] ;
399 }
400 
401 bool element_side_polarity( const Entity & elem ,
402  const Entity & side , int local_side_id )
403 {
404  // 09/14/10: TODO: tscoffe: Will this work in 1D?
405  FEMMetaData &fem_meta = FEMMetaData::get(elem);
406  const bool is_side = side.entity_rank() != fem_meta.edge_rank();
407  const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
408 
409  const unsigned side_count = ! elem_top ? 0 : (
410  is_side ? elem_top->side_count
411  : elem_top->edge_count );
412 
413  ThrowErrorMsgIf( elem_top == NULL,
414  "For Element[" << elem.identifier() << "], element has no defined topology");
415 
416  ThrowErrorMsgIf( local_side_id < 0 || static_cast<int>(side_count) <= local_side_id,
417  "For Element[" << elem.identifier() << "], " <<
418  "side: " << print_entity_key(side) << ", " <<
419  "local_side_id = " << local_side_id <<
420  " ; unsupported local_side_id");
421 
422  const CellTopologyData * const side_top =
423  is_side ? elem_top->side[ local_side_id ].topology
424  : elem_top->edge[ local_side_id ].topology ;
425 
426  const unsigned * const side_map =
427  is_side ? elem_top->side[ local_side_id ].node
428  : elem_top->edge[ local_side_id ].node ;
429 
430  const PairIterRelation elem_nodes = elem.relations( FEMMetaData::NODE_RANK );
431  const PairIterRelation side_nodes = side.relations( FEMMetaData::NODE_RANK );
432 
433  const unsigned n = side_top->node_count;
434  bool good = false ;
435  for ( unsigned i = 0 ; !good && i < n ; ++i ) {
436  good = true;
437  for ( unsigned j = 0; good && j < n ; ++j ) {
438  good = side_nodes[(j+i)%n].entity() == elem_nodes[ side_map[j] ].entity();
439  }
440  }
441  return good ;
442 }
443 
444 }
445 }
446 }
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 declare_relation(Entity &e_from, Entity &e_to, const RelationIdentifier local_id)
Declare a relation and its converse between entities in the same mesh.
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
bool has_superset(const Bucket &bucket, const unsigned &ordinal)
Is this bucket a subset of the given part by partID.
Definition: Bucket.cpp:127
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...
Entity & declare_element(BulkData &mesh, Part &part, const EntityId elem_id, const EntityId node_id[])
Declare an element member of a Part with a CellTopology and nodes conformal to that topology...
Definition: FEMHelpers.cpp:72
Entity & declare_element_edge(BulkData &mesh, const stk_classic::mesh::EntityId global_edge_id, Entity &elem, const unsigned local_edge_id, Part *part)
Create (or find) an element edge.
Definition: FEMHelpers.cpp:211
void all_reduce_sum(ParallelMachine comm, const double *local, double *global, unsigned count)
Parallel summation to all processors.
EntityRank element_rank() const
Returns the element rank which is always equal to spatial dimension.
const std::vector< Bucket * > & buckets(EntityRank rank) const
Query all buckets of a given entity rank.
Definition: BulkData.hpp:195
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
Definition: BulkData.hpp:211
An application-defined subset of a problem domain.
Definition: Part.hpp:49
void change_entity_parts(Entity &entity, const PartVector &add_parts, const PartVector &remove_parts=PartVector())
Change the parallel-locally-owned entity&#39;s part membership by adding and/or removing parts...
Definition: BulkData.hpp:249
ParallelMachine parallel() const
The parallel machine.
Definition: BulkData.hpp:79
bool comm_mesh_counts(BulkData &M, std::vector< size_t > &counts, bool local_flag)
Global counts for a mesh&#39;s entities.
Definition: FEMHelpers.cpp:362
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
const std::string & name() const
Application-defined text name of this part.
Definition: Part.hpp:67
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
Sierra Toolkit.
EntityRank edge_rank() const
Returns the edge rank which changes depending on spatial dimension.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
std::vector< std::string >::size_type entity_rank_count() const
Return the maximum entity rank.
Entity & declare_element_side(BulkData &mesh, const stk_classic::mesh::EntityId global_side_id, Entity &elem, const unsigned local_side_id, Part *part)
Create (or find) an element side.
Definition: FEMHelpers.cpp:186
EntityRank entity_rank() const
The rank of this entity.
Definition: Entity.hpp:128
Entity & declare_entity(EntityRank ent_rank, EntityId ent_id, const PartVector &parts)
Create or retrieve a locally owned entity of a given rank and id.
Definition: BulkData.cpp:215
EntityId identifier() const
Identifier for this entity which is globally unique for a given entity type.
Definition: Entity.hpp:133
EntityRank node_rank() const
Returns the node rank, which is always zero.
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).
bool element_side_polarity(const Entity &elem, const Entity &side, int local_side_id)
Determine the polarity of the local side, more efficient if the local_side_id is known.
Definition: FEMHelpers.cpp:401
fem::CellTopology get_cell_topology(const Part &part) const
Return the cell topology associated with the given part. The cell topology is set on a part through p...