1 #include <stk_mesh/fem/FEMHelpers.hpp> 2 #include <stk_mesh/fem/FEMMetaData.hpp> 4 #include <Shards_CellTopologyTraits.hpp> 6 #include <stk_mesh/base/Types.hpp> 7 #include <stk_mesh/base/BulkData.hpp> 9 #include <stk_util/parallel/ParallelReduce.hpp> 20 void verify_declare_element_side(
21 const BulkData & mesh,
23 const unsigned local_side_id
26 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
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 ;
32 ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
33 "For elem " << print_entity_key(elem) <<
34 ", Bulkdata for 'elem' and mesh are different");
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 );
40 ThrowErrorMsgIf( side_top == NULL,
41 "For elem " << print_entity_key(elem) <<
", local_side_id " << local_side_id <<
", " <<
42 "No element topology found");
45 void verify_declare_element_edge(
46 const BulkData & mesh,
48 const unsigned local_edge_id
51 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
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 ;
57 ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
58 "For elem " << print_entity_key(elem) <<
59 ", Bulkdata for 'elem' and mesh are different");
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 );
65 ThrowErrorMsgIf( edge_top == NULL,
66 "For elem " << print_entity_key(elem) <<
", local_edge_id " << local_edge_id <<
", " <<
67 "No element topology found");
74 const EntityId elem_id ,
75 const EntityId node_id[] )
78 const CellTopologyData *
const top = fem_meta.
get_cell_topology( part ).getCellTopologyData();
80 ThrowErrorMsgIf(top == NULL,
81 "Part " << part.
name() <<
" does not have a local topology");
90 const EntityRank node_rank = fem_meta.
node_rank();
92 for (
unsigned i = 0 ; i < top->node_count ; ++i ) {
107 const unsigned local_side_id ,
110 BulkData & mesh = BulkData::get(side);
112 verify_declare_element_side(mesh, elem, local_side_id);
114 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
116 ThrowErrorMsgIf( elem_top == NULL,
117 "Element[" << elem.
identifier() <<
"] has no defined topology" );
119 const CellTopologyData *
const side_top = elem_top->side[ local_side_id ].topology;
121 ThrowErrorMsgIf( side_top == NULL,
122 "Element[" << elem.
identifier() <<
"], local_side_id = " <<
123 local_side_id <<
", side has no defined topology" );
125 const unsigned *
const side_node_map = elem_top->side[ local_side_id ].node ;
129 if ( part ) { add_parts.push_back( part ); }
137 for (
unsigned i = 0 ; i < side_top->node_count ; ++i ) {
138 Entity & node = * rel[ side_node_map[i] ].entity();
148 const unsigned local_edge_id ,
151 BulkData & mesh = BulkData::get(edge);
155 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
157 ThrowErrorMsgIf( elem_top == NULL,
158 "Element[" << elem.
identifier() <<
"] has no defined topology" );
160 const CellTopologyData *
const edge_top = elem_top->edge[ local_edge_id ].topology;
162 ThrowErrorMsgIf( edge_top == NULL,
163 "Element[" << elem.
identifier() <<
"], local_edge_id = " <<
164 local_edge_id <<
", edge has no defined topology" );
166 const unsigned *
const edge_node_map = elem_top->edge[ local_edge_id ].node ;
170 if ( part ) { add_parts.push_back( part ); }
178 for (
unsigned i = 0 ; i < edge_top->node_count ; ++i ) {
179 Entity & node = * rel[ edge_node_map[i] ].entity();
188 const stk_classic::mesh::EntityId global_side_id ,
190 const unsigned local_side_id ,
193 verify_declare_element_side(mesh, elem, local_side_id);
195 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
197 ThrowErrorMsgIf( elem_top == NULL,
198 "Element[" << elem.
identifier() <<
"] has no defined topology");
200 const CellTopologyData *
const side_top = elem_top->side[ local_side_id ].topology;
202 ThrowErrorMsgIf( side_top == NULL,
203 "Element[" << elem.
identifier() <<
"], local_side_id = " <<
204 local_side_id <<
", side has no defined topology" );
213 const stk_classic::mesh::EntityId global_edge_id ,
215 const unsigned local_edge_id ,
218 verify_declare_element_edge(mesh, elem, local_edge_id);
220 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
222 ThrowErrorMsgIf( elem_top == NULL,
223 "Element[" << elem.
identifier() <<
"] has no defined topology");
226 const CellTopologyData *
const edge_top = elem_top->edge[ local_edge_id ].topology;
228 ThrowErrorMsgIf( edge_top == NULL,
229 "Element[" << elem.
identifier() <<
"], local_edge_id = " <<
230 local_edge_id <<
", edge has no defined topology" );
240 EntityRank subcell_rank ,
241 unsigned subcell_identifier ,
242 EntityVector & subcell_nodes)
244 subcell_nodes.clear();
247 const CellTopologyData* celltopology = get_cell_topology(entity).getCellTopologyData();
252 if (celltopology == NULL) {
257 const bool bad_rank = subcell_rank >= celltopology->dimension;
258 ThrowInvalidArgMsgIf( bad_rank,
"subcell_rank is >= celltopology dimension\n");
261 const bool bad_id = subcell_identifier >= celltopology->subcell_count[subcell_rank];
262 ThrowInvalidArgMsgIf( bad_id,
"subcell_id is >= subcell_count\n");
266 const CellTopologyData * subcell_topology =
267 celltopology->subcell[subcell_rank][subcell_identifier].topology;
269 const int num_nodes_in_subcell = subcell_topology->node_count;
272 const unsigned* subcell_node_local_ids =
273 celltopology->subcell[subcell_rank][subcell_identifier].node;
276 const EntityRank node_rank = fem_meta.
node_rank();
279 subcell_nodes.reserve(num_nodes_in_subcell);
281 for (
int i = 0; i < num_nodes_in_subcell; ++i ) {
282 subcell_nodes.push_back( node_relations[subcell_node_local_ids[i]].entity() );
285 return subcell_topology;
290 const EntityRank subcell_rank,
291 const CellTopologyData * subcell_topology,
292 const std::vector<Entity*>& subcell_nodes )
294 const int INVALID_SIDE = -1;
296 unsigned num_nodes = subcell_topology->node_count;
298 if (num_nodes != subcell_nodes.size()) {
303 const CellTopologyData* entity_topology = get_cell_topology(entity).getCellTopologyData();
304 if (entity_topology == NULL) {
310 const EntityRank node_rank = fem_meta.node_rank();
311 PairIterRelation relations = entity.relations(node_rank);
313 const int num_permutations = subcell_topology->permutation_count;
316 for (
unsigned local_subcell_ordinal = 0;
317 local_subcell_ordinal < entity_topology->subcell_count[subcell_rank];
318 ++local_subcell_ordinal) {
321 const CellTopologyData* curr_subcell_topology =
322 entity_topology->subcell[subcell_rank][local_subcell_ordinal].topology;
325 if (subcell_topology == curr_subcell_topology) {
327 const unsigned*
const subcell_node_map = entity_topology->subcell[subcell_rank][local_subcell_ordinal].node;
333 for (
int p = 0; p < num_permutations; ++p) {
335 if (curr_subcell_topology->permutation[p].polarity ==
336 CELL_PERMUTATION_POLARITY_POSITIVE) {
338 const unsigned *
const perm_node =
339 curr_subcell_topology->permutation[p].node ;
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()) {
352 return local_subcell_ordinal ;
363 std::vector<size_t> & counts ,
366 const size_t zero = 0 ;
372 const size_t comm_count = entity_rank_count + 1 ;
374 std::vector<size_t> local( comm_count , zero );
375 std::vector<size_t> global( comm_count , zero );
380 for (
unsigned i = 0 ; i < entity_rank_count ; ++i ) {
381 const std::vector<Bucket*> & ks = M.
buckets( i );
383 std::vector<Bucket*>::const_iterator ik ;
385 for ( ik = ks.begin() ; ik != ks.end() ; ++ik ) {
387 local[i] += (*ik)->size();
392 local[ entity_rank_count ] = local_flag ;
396 counts.assign( global.begin() , global.begin() + entity_rank_count );
398 return 0 < global[ entity_rank_count ] ;
402 const Entity & side ,
int local_side_id )
407 const CellTopologyData *
const elem_top = get_cell_topology( elem ).getCellTopologyData();
409 const unsigned side_count = ! elem_top ? 0 : (
410 is_side ? elem_top->side_count
411 : elem_top->edge_count );
413 ThrowErrorMsgIf( elem_top == NULL,
414 "For Element[" << elem.
identifier() <<
"], element has no defined topology");
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");
422 const CellTopologyData *
const side_top =
423 is_side ? elem_top->side[ local_side_id ].topology
424 : elem_top->edge[ local_side_id ].topology ;
426 const unsigned *
const side_map =
427 is_side ? elem_top->side[ local_side_id ].node
428 : elem_top->edge[ local_side_id ].node ;
433 const unsigned n = side_top->node_count;
435 for (
unsigned i = 0 ; !good && i < n ; ++i ) {
437 for (
unsigned j = 0; good && j < n ; ++j ) {
438 good = side_nodes[(j+i)%n].entity() == elem_nodes[ side_map[j] ].entity();
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.
bool has_superset(const Bucket &bucket, const unsigned &ordinal)
Is this bucket a subset of the given part by partID.
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...
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.
void all_reduce_sum(ParallelMachine comm, const double *local, double *global, unsigned count)
Parallel summation to all processors.
const std::vector< Bucket * > & buckets(EntityRank rank) const
Query all buckets of a given entity rank.
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
An application-defined subset of a problem domain.
void change_entity_parts(Entity &entity, const PartVector &add_parts, const PartVector &remove_parts=PartVector())
Change the parallel-locally-owned entity's part membership by adding and/or removing parts...
ParallelMachine parallel() const
The parallel machine.
bool comm_mesh_counts(BulkData &M, std::vector< size_t > &counts, bool local_flag)
Global counts for a mesh's entities.
const CellTopologyData * get_subcell_nodes(const Entity &entity, EntityRank subcell_rank, unsigned subcell_identifier, EntityVector &subcell_nodes)
PairIterRelation relations() const
All Entity relations for which this entity is a member. The relations are ordered from lowest entity-...
const std::string & name() const
Application-defined text name of this part.
Manager for an integrated collection of entities, entity relations, and buckets of field data...
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
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.
EntityRank entity_rank() const
The rank of this entity.
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.
EntityId identifier() const
Identifier for this entity which is globally unique for a given entity type.
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
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.