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> 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> 19 #include <stk_util/parallel/ParallelComm.hpp> 26 typedef std::pair< const CellTopologyData *, EntityVector > SideKey;
27 typedef std::vector< EntitySideComponent > SideVector;
28 typedef std::map< SideKey, SideVector> BoundaryMap;
32 class EntitySideComponentLess {
34 bool operator () (
const EntitySideComponent & lhs,
const EntitySideComponent &rhs)
36 const EntityId lhs_elem_id = lhs.entity->identifier();
37 const EntityId rhs_elem_id = rhs.entity->identifier();
39 return (lhs_elem_id != rhs_elem_id) ?
40 (lhs_elem_id < rhs_elem_id) :
41 (lhs.side_ordinal < rhs.side_ordinal);
50 RelationIdentifier side_ordinal;
53 elem_id(0), side_ordinal(0) {}
55 ElementIdSide( EntityId
id, RelationIdentifier ordinal) :
56 elem_id(id), side_ordinal(ordinal) {}
58 ElementIdSide(
const EntitySideComponent & esc) :
59 elem_id(esc.entity->identifier()), side_ordinal(esc.side_ordinal) {}
61 ElementIdSide & operator = (
const ElementIdSide & rhs) {
62 elem_id = rhs.elem_id;
63 side_ordinal = rhs.side_ordinal;
68 bool operator < (
const ElementIdSide & rhs)
const {
69 const ElementIdSide & lhs = *
this;
71 return (lhs.elem_id != rhs.elem_id) ?
72 (lhs.elem_id < rhs.elem_id) :
73 (lhs.side_ordinal < rhs.side_ordinal);
79 typedef std::map< ElementIdSide, SideKey> ReverseBoundaryMap;
82 class SideCommHelper {
85 ElementIdSide creating_elem_id_side;
86 EntityId generated_side_id;
89 proc_to(0), creating_elem_id_side(), generated_side_id(0) {};
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) {}
94 SideCommHelper(
const SideCommHelper & sch) :
96 creating_elem_id_side(sch.creating_elem_id_side),
97 generated_side_id(sch.generated_side_id)
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;
108 bool operator < (
const SideCommHelper & rhs)
const {
109 const SideCommHelper & lhs = *
this;
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);
118 void ensure_consistent_order(EntityVector & side_entities)
120 ThrowRequire( !side_entities.empty() );
122 EntityId lowest_id = side_entities.front()->identifier();
123 unsigned idx_of_lowest_id = 0;
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;
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());
144 void add_owned_sides_to_map(
145 const BulkData & mesh,
146 const EntityRank element_rank,
147 const EntitySideVector & boundary,
148 BoundaryMap & side_map)
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);
157 if ( inside_entity.owner_rank() == mesh.parallel_rank() &&
158 outside.entity == NULL ) {
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);
166 if (existing_sides.first == existing_sides.second) {
180 ensure_consistent_order(side_key.second);
183 side_map[side_key].push_back(inside);
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)
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);
209 if ( inside_entity.owner_rank() != mesh.parallel_rank() &&
210 outside.entity == NULL ) {
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);
218 if (existing_sides.first == existing_sides.second) {
228 ensure_consistent_order(side_key.second);
231 if ( side_map.find(side_key) != side_map.end()) {
232 side_map[side_key].push_back(inside);
252 size_t determine_creating_processes(
253 const BulkData & mesh,
254 BoundaryMap & side_map,
255 ReverseBoundaryMap & reverse_side_map)
257 int num_sides_to_create = 0;
258 for (BoundaryMap::iterator i = side_map.begin();
262 const SideKey & side_key = i->first;
263 SideVector & side_vector = i->second;
266 std::sort( side_vector.begin(), side_vector.end(), EntitySideComponentLess() );
268 const EntitySideComponent & first_side = *side_vector.begin();
272 if (first_side.entity->owner_rank() == mesh.parallel_rank()) {
273 ++num_sides_to_create;
275 const ElementIdSide elem_id_side(first_side);
277 reverse_side_map[elem_id_side] = side_key;
280 return num_sides_to_create;
285 void skin_mesh( BulkData & mesh, EntityRank element_rank, Part * skin_part) {
286 ThrowErrorMsgIf( mesh.synchronized_state() == BulkData::MODIFIABLE,
287 "mesh is not SYNCHRONIZED" );
289 EntityVector owned_elements;
294 mesh.buckets(element_rank),
297 reskin_mesh(mesh, element_rank, owned_elements, skin_part);
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" );
304 EntityVector elements_closure;
307 find_closure( mesh, owned_elements, elements_closure );
310 EntitySideVector boundary;
311 boundary_analysis( mesh, elements_closure, element_rank, boundary);
313 BoundaryMap side_map;
315 add_owned_sides_to_map(mesh, element_rank, boundary, side_map);
317 add_non_owned_sides_to_map(mesh, element_rank, boundary, side_map);
319 ReverseBoundaryMap reverse_side_map;
321 size_t num_sides_to_create = determine_creating_processes(mesh, side_map, reverse_side_map);
324 mesh.modification_begin();
328 requests[element_rank -1] = num_sides_to_create;
331 EntityVector requested_sides;
332 mesh.generate_new_entities(requests, requested_sides);
335 std::set<SideCommHelper> side_comm_helper_set;
338 size_t current_side = 0;
339 for ( BoundaryMap::iterator map_itr = side_map.begin();
340 map_itr!= side_map.end();
343 const SideKey & side_key = map_itr->first;
344 SideVector & side_vector = map_itr->second;
348 const EntitySideComponent & first_side = *(side_vector.begin());
349 if ( first_side.entity->owner_rank() == mesh.parallel_rank()) {
356 const ElementIdSide elem_id_side( first_side);
358 Entity & side = *(requested_sides[current_side]);
359 EntityId side_id = side.identifier();
365 Part * topo_part = &fem_meta_data.get_cell_topology_root_part(side_key.first);
366 add_parts.push_back( topo_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) {
373 if (topo == fem_meta_data.get_cell_topology(**i)) {
374 if (std::string(side_key.first->name) ==
"Quadrilateral_4") {
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")) {
378 add_parts.push_back(*i);
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")) {
382 add_parts.push_back(*i);
386 add_parts.push_back(*i);
392 mesh.change_entity_parts(side, add_parts);
396 const EntityVector & nodes = side_key.second;
398 for (
size_t i = 0; i<nodes.size(); ++i) {
399 Entity & node = *nodes[i];
400 mesh.declare_relation( side, node, i);
404 for (SideVector::iterator side_itr = side_vector.begin();
405 side_itr != side_vector.end();
408 Entity & elem = *(side_itr->entity);
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);
416 side_comm_helper_set.insert( SideCommHelper( elem.owner_rank(), elem_id_side, side_id));
423 CommAll comm( mesh.parallel() );
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);
431 for (std::set<SideCommHelper>::const_iterator i=side_comm_helper_set.begin();
432 i != side_comm_helper_set.end();
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);
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;
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);
456 const SideKey & side_key = reverse_side_map[creating_elem_id_side];
459 SideVector & side_vector = side_map[side_key];
464 Part * topo_part = &fem_meta_data.get_cell_topology_root_part(side_key.first);
465 add_parts.push_back( topo_part);
467 add_parts.push_back(skin_part);
468 fem::CellTopology topo = fem_meta_data.get_cell_topology(*topo_part);
470 for (stk_classic::mesh::PartVector::const_iterator i=topo_parts.begin(); i!=topo_parts.end(); ++i) {
472 if (topo == fem_meta_data.get_cell_topology(**i)) {
473 if (std::string(side_key.first->name) ==
"Quadrilateral_4") {
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")) {
477 add_parts.push_back(*i);
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")) {
481 add_parts.push_back(*i);
485 add_parts.push_back(*i);
491 Entity & side = mesh.declare_entity(element_rank-1,generated_side_id,add_parts);
494 const EntityVector & nodes = side_key.second;
496 for (
size_t i = 0; i<nodes.size(); ++i) {
497 Entity & node = *nodes[i];
498 mesh.declare_relation( side, node, i);
502 for (SideVector::iterator side_itr = side_vector.begin();
503 side_itr != side_vector.end();
506 Entity & elem = *(side_itr->entity);
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);
515 mesh.modification_end();
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...
const CellTopologyData * get_subcell_nodes(const Entity &entity, EntityRank subcell_rank, unsigned subcell_identifier, EntityVector &subcell_nodes)
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.