44 #include "Panzer_Workset_Builder.hpp"
45 #include "Teuchos_Assert.hpp"
47 #include <stk_mesh/base/Selector.hpp>
48 #include <stk_mesh/base/GetEntities.hpp>
50 namespace panzer_stk {
54 const std::string & eBlock,
57 using namespace workset_utils;
59 std::vector<std::string> element_blocks;
61 std::vector<std::size_t> local_cell_ids;
62 Kokkos::DynRankView<double,PHX::Device> cell_node_coordinates;
64 getIdsAndNodes(mesh, eBlock, local_cell_ids, cell_node_coordinates);
76 const std::string & sideset,
77 const std::string & eBlock,
80 using namespace workset_utils;
83 std::vector<stk::mesh::Entity> sideEntities;
95 std::vector<std::string> sideSets;
99 ss << e.what() <<
"\nChoose one of:\n";
100 for(std::size_t i=0;i<sideSets.size();i++)
101 ss <<
"\"" << sideSets[i] <<
"\"\n";
106 std::stringstream ss;
107 std::vector<std::string> elementBlocks;
111 ss << e.what() <<
"\nChoose one of:\n";
112 for(std::size_t i=0;i<elementBlocks.size();i++)
113 ss <<
"\"" << elementBlocks[i] <<
"\"\n";
117 catch(std::logic_error & e) {
118 std::stringstream ss;
119 ss << e.what() <<
"\nUnrecognized logic error.\n";
124 std::vector<stk::mesh::Entity> elements;
125 std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > local_cell_ids;
128 std::vector<std::size_t> local_side_ids;
130 sideEntities,local_side_ids,elements);
133 for(std::size_t elm=0;elm<elements.size();++elm) {
134 stk::mesh::Entity element = elements[elm];
136 local_cell_ids[std::make_pair(subcell_dim,local_side_ids[elm])].push_back(mesh.
elementLocalId(element));
140 std::vector<std::size_t> local_subcell_ids, subcell_dim;
142 sideEntities,subcell_dim,local_subcell_ids,elements);
145 for(std::size_t elm=0;elm<elements.size();++elm) {
146 stk::mesh::Entity element = elements[elm];
148 local_cell_ids[std::make_pair(subcell_dim[elm],local_subcell_ids[elm])].push_back(mesh.
elementLocalId(element));
156 if(elements.size()!=0) {
163 for(std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator itr=local_cell_ids.begin();
164 itr!=local_cell_ids.end();++itr) {
166 if(itr->second.size()==0)
169 Kokkos::DynRankView<double,PHX::Device> nodes;
176 for(std::size_t w=0;w<current->size();w++) {
177 (*current)[w].subcell_dim = itr->first.first;
178 (*current)[w].subcell_index = itr->first.second;
182 worksets->insert(worksets->end(),current->begin(),current->end());
189 return Teuchos::rcp(
new std::vector<panzer::Workset>());
195 const std::string & blockid_a,
197 const std::string & blockid_b,
198 const std::string & sideset)
200 using namespace workset_utils;
203 std::vector<stk::mesh::Entity> sideEntities;
213 stk::mesh::Part * sidePart = mesh.
getSideset(sideset);
215 "Unknown side set \"" << sideset <<
"\"");
217 stk::mesh::Selector side = *sidePart;
224 std::stringstream ss;
225 std::vector<std::string> elementBlocks;
229 ss << e.what() <<
"\nChoose one of:\n";
230 for(std::size_t i=0;i<elementBlocks.size();i++)
231 ss <<
"\"" << elementBlocks[i] <<
"\"\n";
235 catch(std::logic_error & e) {
236 std::stringstream ss;
237 ss << e.what() <<
"\nUnrecognized logic error.\n";
242 std::vector<stk::mesh::Entity> elements_a, elements_b;
243 std::vector<std::size_t> local_cell_ids_a, local_cell_ids_b;
244 std::vector<std::size_t> local_side_ids_a, local_side_ids_b;
248 local_side_ids_a,elements_a,
249 local_side_ids_b,elements_b);
252 "For a DG type boundary, the number of elements on the \"left\" and \"right\" is not the same.");
258 if(elements_a.size()==0)
259 return Teuchos::rcp(
new std::map<unsigned,panzer::Workset>);
264 for(std::size_t elm=0;elm<elements_a.size();++elm) {
265 stk::mesh::Entity element_a = elements_a[elm];
266 stk::mesh::Entity element_b = elements_b[elm];
272 Kokkos::DynRankView<double,PHX::Device> node_coordinates_a, node_coordinates_b;
277 return buildBCWorkset(needs_a,blockid_a, local_cell_ids_a, local_side_ids_a, node_coordinates_a,
278 needs_b,blockid_b, local_cell_ids_b, local_side_ids_b, node_coordinates_b);
284 const std::string & eblockID,
285 const std::string & sidesetID)
287 using namespace workset_utils;
290 std::vector<stk::mesh::Entity> sideEntities;
295 mesh.
getMySides(sidesetID,eblockID,sideEntities);
298 std::stringstream ss;
299 std::vector<std::string> sideSets;
303 ss << e.what() <<
"\nChoose one of:\n";
304 for(std::size_t i=0;i<sideSets.size();i++)
305 ss <<
"\"" << sideSets[i] <<
"\"\n";
310 std::stringstream ss;
311 std::vector<std::string> elementBlocks;
315 ss << e.what() <<
"\nChoose one of:\n";
316 for(std::size_t i=0;i<elementBlocks.size();i++)
317 ss <<
"\"" << elementBlocks[i] <<
"\"\n";
321 catch(std::logic_error & e) {
322 std::stringstream ss;
323 ss << e.what() <<
"\nUnrecognized logic error.\n";
328 std::vector<stk::mesh::Entity> elements;
329 std::vector<std::size_t> local_cell_ids;
330 std::vector<std::size_t> local_side_ids;
332 sideEntities,local_side_ids,elements);
335 for(std::size_t elm=0;elm<elements.size();++elm) {
336 stk::mesh::Entity element = elements[elm];
345 if(elements.size()!=0) {
349 Kokkos::DynRankView<double,PHX::Device> nodes;
355 return Teuchos::null;
358 namespace workset_utils {
361 const std::string & blockId,
362 const std::vector<stk::mesh::Entity> & entities,
363 std::vector<std::size_t> & localEntityIds,
364 std::vector<stk::mesh::Entity> & elements)
369 stk::mesh::BulkData& bulkData = *mesh.
getBulkData();
373 std::vector<stk::mesh::Entity>::const_iterator entityItr;
374 for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
375 stk::mesh::Entity entity = *entityItr;
377 const size_t num_rels = bulkData.num_elements(entity);
378 stk::mesh::Entity
const* relations = bulkData.begin_elements(entity);
379 stk::mesh::ConnectivityOrdinal
const* ordinals = bulkData.begin_element_ordinals(entity);
380 for(std::size_t e=0; e<num_rels; ++e) {
381 stk::mesh::Entity element = relations[e];
382 std::size_t entityId = ordinals[e];
385 stk::mesh::Bucket
const& bucket = bulkData.bucket(element);
386 bool inBlock = bucket.member(*blockPart);
387 bool onProc = bucket.member(*ownedPart);
388 if(inBlock && onProc) {
390 elements.push_back(element);
391 localEntityIds.push_back(entityId);
398 const std::string & blockId,
399 const std::vector<stk::mesh::Entity> & entities,
400 std::vector<std::size_t> & localEntityIds,
401 std::vector<stk::mesh::Entity> & elements,
402 std::vector<std::size_t> & missingElementIndices)
406 stk::mesh::Part * universalPart = &mesh.
getMetaData()->universal_part();
407 stk::mesh::BulkData& bulkData = *mesh.
getBulkData();
411 std::size_t entityIndex =-1;
412 std::vector<stk::mesh::Entity>::const_iterator entityItr;
413 for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
414 stk::mesh::Entity entity = *entityItr;
417 const size_t num_rels = bulkData.num_elements(entity);
418 stk::mesh::Entity
const* element_rels = bulkData.begin_elements(entity);
419 stk::mesh::ConnectivityOrdinal
const* ordinals = bulkData.begin_element_ordinals(entity);
420 for(std::size_t e=0; e<num_rels; ++e) {
421 stk::mesh::Entity element = element_rels[e];
422 std::size_t entityId = ordinals[e];
425 stk::mesh::Bucket
const& bucket = bulkData.bucket(element);
426 bool inBlock = bucket.member(*blockPart);
427 bool onProc = bucket.member(*universalPart);
428 if(inBlock && onProc) {
430 elements.push_back(element);
431 localEntityIds.push_back(entityId);
432 }
else if(!inBlock && (num_rels == 1)) {
435 missingElementIndices.push_back(entityIndex);
442 const std::string & blockId,
443 const std::vector<stk::mesh::Entity> & sides,
444 std::vector<std::size_t> & localSubcellDim,
445 std::vector<std::size_t> & localSubcellIds,
446 std::vector<stk::mesh::Entity> & elements)
453 std::vector<std::vector<stk::mesh::Entity> > subcells;
455 subcells.push_back(sides);
460 for(std::size_t d=0;d<subcells.size();d++) {
461 std::vector<std::size_t> subcellIds;
462 std::vector<stk::mesh::Entity> subcellElements;
471 localSubcellDim.insert(localSubcellDim.end(),subcellElements.size(),d);
472 localSubcellIds.insert(localSubcellIds.end(),subcellIds.begin(),subcellIds.end());
473 elements.insert(elements.end(),subcellElements.begin(),subcellElements.end());
478 const std::string & blockId,
479 const std::vector<stk::mesh::Entity> & sides,
480 std::vector<std::size_t> & localSideIds,
481 std::vector<stk::mesh::Entity> & elements)
487 const std::string & blockId_a,
488 const std::string & blockId_b,
489 const std::vector<stk::mesh::Entity> & sides,
490 std::vector<std::size_t> & localSideIds_a,
491 std::vector<stk::mesh::Entity> & elements_a,
492 std::vector<std::size_t> & localSideIds_b,
493 std::vector<stk::mesh::Entity> & elements_b)
499 stk::mesh::Part * universalPart = &mesh.
getMetaData()->universal_part();
500 stk::mesh::BulkData& bulkData = *mesh.
getBulkData();
504 std::vector<stk::mesh::Entity>::const_iterator sidesItr;
505 for(sidesItr=sides.begin();sidesItr!=sides.end();++sidesItr) {
506 stk::mesh::Entity side = *sidesItr;
509 stk::mesh::Entity element_a = stk::mesh::Entity(), element_b = stk::mesh::Entity();
510 std::size_t entityId_a=0, entityId_b=0;
512 const size_t num_rels = bulkData.num_elements(side);
513 stk::mesh::Entity
const* element_rels = bulkData.begin_elements(side);
514 stk::mesh::ConnectivityOrdinal
const* ordinals = bulkData.begin_element_ordinals(side);
515 for(std::size_t e=0; e<num_rels; ++e) {
516 stk::mesh::Entity element = element_rels[e];
517 std::size_t entityId = ordinals[e];
520 stk::mesh::Bucket
const& bucket = bulkData.bucket(element);
521 bool inBlock_a = bucket.member(*blockPart_a);
522 bool inBlock_b = bucket.member(*blockPart_b);
523 bool onProc = bucket.member(*ownedPart);
524 bool unProc = bucket.member(*universalPart);
526 if(inBlock_a && onProc) {
529 entityId_a = entityId;
531 if(inBlock_b && unProc) {
534 entityId_b = entityId;
538 if(element_a!=stk::mesh::Entity() && element_b!=stk::mesh::Entity()) {
539 elements_a.push_back(element_a);
540 localSideIds_a.push_back(entityId_a);
543 elements_b.push_back(element_b);
544 localSideIds_b.push_back(entityId_b);
550 const std::string & blockId,
551 const std::vector<stk::mesh::Entity> & nodes,
552 std::vector<std::size_t> & localNodeIds,
553 std::vector<stk::mesh::Entity> & elements)
559 const std::vector<stk::mesh::Entity> & entities,
560 std::vector<std::vector<stk::mesh::Entity> > & subcells)
563 if(entities.size()==0) {
568 stk::mesh::BulkData& bulkData = *mesh.
getBulkData();
569 stk::mesh::EntityRank master_rank = bulkData.entity_rank(entities[0]);
571 std::vector<std::set<stk::mesh::Entity> > subcells_set(master_rank);
575 std::vector<stk::mesh::Entity>::const_iterator entityItr;
576 for(entityItr=entities.begin();entityItr!=entities.end();++entityItr) {
577 stk::mesh::Entity entity = *entityItr;
582 for(
int i=0; i<static_cast<int>(master_rank); i++) {
583 stk::mesh::EntityRank
const to_rank =
static_cast<stk::mesh::EntityRank
>(i);
584 const size_t num_rels = bulkData.num_connectivity(entity, to_rank);
585 stk::mesh::Entity
const* relations = bulkData.begin(entity, to_rank);
589 for(std::size_t e=0; e<num_rels; ++e) {
590 stk::mesh::Entity subcell = relations[e];
592 subcells_set[i].insert(subcell);
599 subcells.resize(subcells_set.size());
600 for(std::size_t i=0;i<subcells_set.size();i++)
601 subcells[i].assign(subcells_set[i].begin(),subcells_set[i].end());
void getSidesetNames(std::vector< std::string > &name) const
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
void getElementBlockNames(std::vector< std::string > &names) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
void getSubcellEntities(const panzer_stk::STK_Interface &mesh, const std::vector< stk::mesh::Entity > &entities, std::vector< std::vector< stk::mesh::Entity > > &subcells)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void getElementNodes(const std::vector< std::size_t > &localIds, ArrayT &nodes) const
Teuchos::RCP< std::vector< Workset > > buildWorksets(const WorksetNeeds &needs, const std::string &elementBlock, const std::vector< std::size_t > &local_cell_ids, const ArrayT &node_coordinates)
Teuchos::RCP< std::vector< panzer::Workset > > buildWorksets(const panzer_stk::STK_Interface &mesh, const std::string &eBlock, const panzer::WorksetNeeds &needs)
void getIdsAndNodes(const panzer_stk::STK_Interface &mesh, std::string blockId, std::vector< std::size_t > &localIds, ArrayT &nodes)
void getSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements)
Teuchos::RCP< std::map< unsigned, panzer::Workset > > buildBCWorksets(const panzer_stk::STK_Interface &mesh, const panzer::WorksetNeeds &needs_a, const std::string &blockid_a, const panzer::WorksetNeeds &needs_b, const std::string &blockid_b, const std::string &sideset)
stk::mesh::Part * getSideset(const std::string &name) const
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void getSideElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &sides, std::vector< std::size_t > &localSideIds, std::vector< stk::mesh::Entity > &elements)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
stk::mesh::EntityRank getSideRank() const
void getUniversalSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements, std::vector< std::size_t > &missingElementIndices)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
Teuchos::RCP< std::map< unsigned, Workset > > buildBCWorkset(const WorksetNeeds &needs, const std::string &elementBlock, const std::vector< std::size_t > &local_cell_ids, const std::vector< std::size_t > &local_side_ids, const ArrayT &node_coordinates, const bool populate_value_arrays=true)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
#define TEUCHOS_ASSERT(assertion_test)
void getSideElementCascade(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &sides, std::vector< std::size_t > &localSubcellDim, std::vector< std::size_t > &localSubcellIds, std::vector< stk::mesh::Entity > &elements)
void getNodeElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &nodes, std::vector< std::size_t > &localNodeIds, std::vector< stk::mesh::Entity > &elements)
stk::mesh::Part * getOwnedPart() const
Get a pointer to the locally owned part.
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
#define TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(throw_exception_test, Exception, msg)