11 #include <PanzerAdaptersSTK_config.hpp>
15 #include <Teuchos_RCPStdSharedPtrConversions.hpp>
19 #include <stk_mesh/base/FieldBase.hpp>
20 #include <stk_mesh/base/Comm.hpp>
21 #include <stk_mesh/base/Selector.hpp>
22 #include <stk_mesh/base/GetEntities.hpp>
23 #include <stk_mesh/base/MeshBuilder.hpp>
24 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
31 #include <stk_util/parallel/ParallelReduce.hpp>
32 #include <stk_util/parallel/CommSparse.hpp>
34 #include <stk_tools/mesh_tools/FixNodeSharingViaSearch.hpp>
36 #ifdef PANZER_HAVE_IOSS
37 #include <Ionit_Initializer.h>
38 #include <stk_io/IossBridge.hpp>
39 #include <stk_io/WriteMesh.hpp>
42 #ifdef PANZER_HAVE_PERCEPT
43 #include <percept/PerceptMesh.hpp>
44 #include <adapt/UniformRefinerPattern.hpp>
45 #include <adapt/UniformRefiner.hpp>
55 namespace panzer_stk {
59 : gid_(gid), nodes_(nodes) {}
78 : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
85 : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
92 : dimension_(dim), initialized_(false), currentLocalId_(0), useFieldCoordinates_(false)
94 std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
95 entity_rank_names.push_back(
"FAMILY_TREE");
108 stk::mesh::Part * sideset =
metaData_->get_part(name);
110 sideset = &
metaData_->declare_part_with_topology(name,
111 stk::mesh::get_topology(shards::CellTopology(ctData),
dimension_));
112 sidesets_.insert(std::make_pair(name,sideset));
120 stk::mesh::Part * nodeset =
metaData_->get_part(name);
121 if(nodeset==
nullptr) {
122 const CellTopologyData * ctData = shards::getCellTopologyData<shards::Node>();
123 nodeset = &
metaData_->declare_part_with_topology(name,
124 stk::mesh::get_topology(shards::CellTopology(ctData),
dimension_));
126 nodesets_.insert(std::make_pair(name,nodeset));
132 "Unknown element block \"" << blockId <<
"\"");
133 std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
139 field = &
metaData_->declare_field<
double>(stk::topology::NODE_RANK, fieldName);
142 stk::mesh::put_field_on_mesh(*field,
metaData_->universal_part(),
nullptr);
151 "Unknown element block \"" << blockId <<
"\"");
152 std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
158 field = &
metaData_->declare_field<
double>(stk::topology::ELEMENT_RANK, fieldName);
162 stk::mesh::put_field_on_mesh(*field,
metaData_->universal_part(),
nullptr);
171 "Unknown element block \"" << blockId <<
"\"");
172 std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
178 field = &
metaData_->declare_field<
double>(stk::topology::EDGE_RANK, fieldName);
183 stk::mesh::put_field_on_mesh(*field,
metaData_->universal_part(),
nullptr);
192 "Unknown element block \"" << blockId <<
"\"");
193 std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
199 field = &
metaData_->declare_field<
double>(stk::topology::FACE_RANK, fieldName);
204 stk::mesh::put_field_on_mesh(*field,
metaData_->universal_part(),
nullptr);
211 const std::vector<std::string> & coordNames,
212 const std::string & dispPrefix)
218 "Unknown element block \"" << blockId <<
"\"");
222 "STK_Interface::addMeshCoordFields: Can't set more than one set of coordinate "
223 "fields for element block \""+blockId+
"\".");
236 std::pair<std::string,std::string> key = std::make_pair(coordNames[i],blockId);
237 std::string dispName = dispPrefix+coordNames[i];
239 dispFields[i] = dispName;
247 field = &
metaData_->declare_field<
double>(stk::topology::NODE_RANK, dispName);
260 const bool buildRefinementSupport)
283 #ifdef PANZER_HAVE_IOSS
290 std::map<std::string, stk::mesh::Part*>::iterator itr;
292 if(!stk::io::is_part_io_part(*itr->second))
293 stk::io::put_io_part_attribute(*itr->second);
298 std::map<std::string, stk::mesh::Part*>::iterator itr;
300 if(!stk::io::is_part_edge_block_io_part(*itr->second))
301 stk::io::put_edge_block_io_part_attribute(*itr->second);
306 std::map<std::string, stk::mesh::Part*>::iterator itr;
308 if(!stk::io::is_part_face_block_io_part(*itr->second))
309 stk::io::put_face_block_io_part_attribute(*itr->second);
314 std::map<std::string, stk::mesh::Part*>::iterator itr;
316 if(!stk::io::is_part_io_part(*itr->second))
317 stk::io::put_io_part_attribute(*itr->second);
322 std::map<std::string, stk::mesh::Part*>::iterator itr;
324 if(!stk::io::is_part_io_part(*itr->second))
325 stk::io::put_io_part_attribute(*itr->second);
333 stk::io::set_field_output_type(*
coordinatesField_, stk::io::FieldOutputType::VECTOR_3D);
334 stk::io::set_field_role(*
edgesField_, Ioss::Field::MESH);
335 stk::io::set_field_output_type(*
edgesField_, stk::io::FieldOutputType::VECTOR_3D);
337 stk::io::set_field_role(*
facesField_, Ioss::Field::MESH);
338 stk::io::set_field_output_type(*
facesField_, stk::io::FieldOutputType::VECTOR_3D);
346 if (buildRefinementSupport) {
347 #ifdef PANZER_HAVE_PERCEPT
352 breakPattern_ =
Teuchos::rcp(
new percept::URP_Heterogeneous_3D(*refinedMesh_));
355 "ERROR: Uniform refinement requested. This requires the Percept package to be enabled in Trilinos!");
370 std::set<SolutionFieldType*> uniqueFields;
371 std::map<std::pair<std::string,std::string>,
SolutionFieldType*>::const_iterator fieldIter;
372 for(fieldIter=nameToField.begin();fieldIter!=nameToField.end();++fieldIter)
373 uniqueFields.insert(fieldIter->second);
376 std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
377 for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
378 stk::mesh::put_field_on_mesh(*(*uniqueFieldIter),
metaData_->universal_part(),
nullptr);
381 #ifdef PANZER_HAVE_IOSS
384 std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
385 for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
386 stk::io::set_field_role(*(*uniqueFieldIter), Ioss::Field::TRANSIENT);
397 std::unique_ptr<stk::mesh::BulkData> bulkUPtr = stk::mesh::MeshBuilder(*
mpiComm_->getRawMpiComm()).create(Teuchos::get_shared_ptr(
metaData_));
404 "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"beginModification\"");
412 "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"endModification\"");
421 if (find_and_set_shared_nodes_in_stk) {
422 const double radius = 10.0 * std::numeric_limits<double>::epsilon();
423 stk::tools::fix_node_sharing_via_search(*
bulkData_,radius);
435 "STK_Interface::addNode: STK_Interface must be modifiable to add a node");
437 "STK_Interface::addNode: number of coordinates in vector must mation dimension");
439 "STK_Interface::addNode: STK has STUPID restriction of no zero GIDs, pick something else");
446 for(std::size_t i=0;i<coord.size();++i)
447 fieldCoords[i] = coord[i];
452 std::vector<stk::mesh::Part*> sidesetV;
453 sidesetV.push_back(sideset);
455 bulkData_->change_entity_parts(entity,sidesetV);
460 std::vector<stk::mesh::Part*> nodesetV;
461 nodesetV.push_back(nodeset);
463 bulkData_->change_entity_parts(entity,nodesetV);
468 std::vector<stk::mesh::Part*> edgeblockV;
469 edgeblockV.push_back(edgeblock);
471 bulkData_->change_entity_parts(entity,edgeblockV);
475 if (entities.size() > 0) {
476 std::vector<stk::mesh::Part*> edgeblockV;
477 edgeblockV.push_back(edgeblock);
479 bulkData_->change_entity_parts(entities,edgeblockV);
485 std::vector<stk::mesh::Part*> faceblockV;
486 faceblockV.push_back(faceblock);
488 bulkData_->change_entity_parts(entity,faceblockV);
492 if (entities.size() > 0) {
493 std::vector<stk::mesh::Part*> faceblockV;
494 faceblockV.push_back(faceblock);
496 bulkData_->change_entity_parts(entities,faceblockV);
502 std::vector<stk::mesh::Part*> blockVec;
503 blockVec.push_back(block);
507 stk::mesh::Entity element =
bulkData_->declare_entity(elementRank,ed->getGID(),blockVec);
510 const std::vector<stk::mesh::EntityId> & nodes = ed->getNodes();
511 for(std::size_t i=0;i<nodes.size();++i) {
513 stk::mesh::Entity node =
bulkData_->get_entity(nodeRank,nodes[i]);
515 bulkData_->declare_relation(element,node,i);
519 procId[0] = Teuchos::as<ProcIdData>(
procRank_);
527 std::vector<stk::mesh::Entity> localElmts;
529 std::vector<stk::mesh::Entity>::const_iterator itr;
530 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
531 stk::mesh::Entity element = (*itr);
532 stk::mesh::EntityId gid =
bulkData_->identifier(element);
533 std::vector<stk::mesh::EntityId> subcellIds;
536 for(std::size_t i=0;i<subcellIds.size();++i) {
537 stk::mesh::Entity edge =
bulkData_->get_entity(edgeRank,subcellIds[i]);
538 stk::mesh::Entity
const* relations =
bulkData_->begin_nodes(edge);
544 double * edgeCoords = stk::mesh::field_data(*
edgesField_,edge);
546 edgeCoords[j] = (node_coord_1[j]+node_coord_2[j])/2.0;
555 std::vector<stk::mesh::Entity> localElmts;
557 std::vector<stk::mesh::Entity>::const_iterator itr;
558 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
559 stk::mesh::Entity element = (*itr);
561 std::vector<stk::mesh::EntityId> subcellIds;
564 for(std::size_t i=0;i<subcellIds.size();++i) {
565 stk::mesh::Entity face =
bulkData_->get_entity(faceRank,subcellIds[i]);
566 stk::mesh::Entity
const* relations =
bulkData_->begin_nodes(face);
567 const size_t num_relations =
bulkData_->num_nodes(face);
570 double * faceCoords = stk::mesh::field_data(*
facesField_,face);
573 for(std::size_t k=0;k<num_relations;++k)
575 faceCoords[j] /= double(num_relations);
583 const size_t num_rels =
bulkData_->num_connectivity(src, tgt_rank);
584 stk::mesh::Entity
const* relations =
bulkData_->begin(src, tgt_rank);
585 stk::mesh::ConnectivityOrdinal
const* ordinals =
bulkData_->begin_ordinals(src, tgt_rank);
586 for (
size_t i = 0; i < num_rels; ++i) {
587 if (ordinals[i] == static_cast<stk::mesh::ConnectivityOrdinal>(rel_id)) {
592 return stk::mesh::Entity();
618 const bool append_after_restart_time,
619 const double restart_time)
621 std::vector<Ioss::Property> ioss_properties;
625 append_after_restart_time,
632 const std::vector<Ioss::Property>& ioss_properties,
634 const bool append_after_restart_time,
635 const double restart_time)
637 using std::runtime_error;
638 using stk::io::StkMeshIoBroker;
639 using stk::mesh::FieldVector;
640 using stk::ParallelMachine;
642 PANZER_FUNC_TIME_MONITOR(
"STK_Interface::setupExodusFile(filename)");
643 #ifdef PANZER_HAVE_IOSS
646 ParallelMachine comm = *
mpiComm_->getRawMpiComm();
647 meshData_ =
rcp(
new StkMeshIoBroker(comm));
648 meshData_->set_bulk_data(Teuchos::get_shared_ptr(
bulkData_));
649 Ioss::PropertyManager props;
650 props.add(Ioss::Property(
"LOWER_CASE_VARIABLE_NAMES",
"FALSE"));
651 for (
const auto& p : ioss_properties ) {
655 if (append_after_restart_time) {
656 meshIndex_ = meshData_->create_output_mesh(filename, stk::io::APPEND_RESULTS,
657 props, restart_time);
660 meshIndex_ = meshData_->create_output_mesh(filename, stk::io::APPEND_RESULTS, props);
663 meshIndex_ = meshData_->create_output_mesh(filename, stk::io::WRITE_RESULTS, props);
664 const FieldVector& fields =
metaData_->get_fields();
665 for (
size_t i(0); i < fields.size(); ++i) {
671 auto role = stk::io::get_field_role(*fields[i]);
672 if (role !=
nullptr) {
673 if (*role != Ioss::Field::MESH)
674 meshData_->add_field(meshIndex_, *fields[i]);
676 meshData_->add_field(meshIndex_, *fields[i]);
682 meshData_->add_info_records(meshIndex_, deduped_info_records);
701 using Teuchos::rcpFromRef;
702 PANZER_FUNC_TIME_MONITOR(
"STK_Interface::writeToExodus(timestep)");
703 #ifdef PANZER_HAVE_IOSS
704 if (not meshData_.is_null())
707 globalToExodus(GlobalVariable::ADD);
709 meshData_->write_defined_output_fields(meshIndex_);
710 globalToExodus(GlobalVariable::WRITE);
711 meshData_->end_output_step(meshIndex_);
716 out.setOutputToRootOnly(0);
717 out <<
"WARNING: Exodus I/O has been disabled or not setup properly; "
718 <<
"not writing to Exodus." << endl;
733 const GlobalVariable& flag)
735 using std::invalid_argument;
742 for (
auto i = globalData_.begin(); i != globalData_.end(); ++i)
744 const string& name = globalData_.name(i);
747 if (globalData_.isType<
int>(name))
749 const auto& value = globalData_.get<
int>(name);
750 if (flag == GlobalVariable::ADD)
754 meshData_->add_global(meshIndex_, name, value,
755 stk::util::ParameterType::INTEGER);
763 meshData_->write_global(meshIndex_, name, value,
764 stk::util::ParameterType::INTEGER);
768 else if (globalData_.isType<
double>(name))
770 const auto& value = globalData_.get<
double>(name);
771 if (flag == GlobalVariable::ADD)
775 meshData_->add_global(meshIndex_, name, value,
776 stk::util::ParameterType::DOUBLE);
784 meshData_->write_global(meshIndex_, name, value,
785 stk::util::ParameterType::DOUBLE);
789 else if (globalData_.isType<Array<int>>(name))
791 const auto& value = globalData_.get<Array<int>>(name).toVector();
792 if (flag == GlobalVariable::ADD)
796 meshData_->add_global(meshIndex_, name, value,
797 stk::util::ParameterType::INTEGERVECTOR);
805 meshData_->write_global(meshIndex_, name, value,
806 stk::util::ParameterType::INTEGERVECTOR);
810 else if (globalData_.isType<Array<double>>(name))
812 const auto& value = globalData_.get<Array<double>>(name).toVector();
813 if (flag == GlobalVariable::ADD)
817 meshData_->add_global(meshIndex_, name, value,
818 stk::util::ParameterType::DOUBLEVECTOR);
826 meshData_->write_global(meshIndex_, name, value,
827 stk::util::ParameterType::DOUBLEVECTOR);
833 "STK_Interface::globalToExodus(): The global variable to be added " \
834 "to the Exodus output file is of an invalid type. Valid types are " \
835 "int and double, along with std::vectors of those types.")
847 const std::string& key,
850 globalData_.set(key, value);
861 const std::string& key,
864 globalData_.set(key, value);
875 const std::string& key,
876 const std::vector<int>& value)
879 globalData_.set(key, Array<int>(value));
890 const std::string& key,
891 const std::vector<double>& value)
894 globalData_.set(key, Array<double>(value));
899 #ifdef PANZER_HAVE_IOSS
911 stk::mesh::Entity node =
bulkData_->get_entity(nodeRank,nodeId);
912 const size_t numElements =
bulkData_->num_elements(node);
913 stk::mesh::Entity
const* relations =
bulkData_->begin_elements(node);
916 elements.insert(elements.end(), relations, relations + numElements);
920 std::vector<int> & relIds)
const
923 const size_t numElements =
bulkData_->num_elements(node);
924 stk::mesh::Entity
const* relations =
bulkData_->begin_elements(node);
925 stk::mesh::ConnectivityOrdinal
const* rel_ids =
bulkData_->begin_element_ordinals(node);
928 for (
size_t i = 0; i < numElements; ++i) {
929 stk::mesh::Entity element = relations[i];
933 elements.push_back(element);
934 relIds.push_back(rel_ids[i]);
940 std::vector<int> & relIds,
unsigned int matchType)
const
942 stk::mesh::EntityRank rank;
945 else if(matchType == 1)
947 else if(matchType == 2)
952 stk::mesh::Entity node =
bulkData_->get_entity(rank,nodeId);
959 std::vector<stk::mesh::Entity> current;
962 std::sort(current.begin(),current.end());
965 for(std::size_t n=1;n<nodeIds.size();++n) {
967 std::vector<stk::mesh::Entity> nextNode;
969 std::sort(nextNode.begin(),nextNode.end());
972 std::vector<stk::mesh::Entity> intersection(std::min(nextNode.size(),current.size()));
973 std::vector<stk::mesh::Entity>::const_iterator endItr
974 = std::set_intersection(current.begin(),current.end(),
975 nextNode.begin(),nextNode.end(),
976 intersection.begin());
977 std::size_t newLength = endItr-intersection.begin();
978 intersection.resize(newLength);
982 current = intersection;
991 stk::mesh::Entity
const* nodeRel =
getBulkData()->begin_nodes(element);
992 const size_t numNodes =
getBulkData()->num_nodes(element);
994 nodeIds.reserve(numNodes);
995 for(
size_t i = 0; i < numNodes; ++i) {
1010 const auto entityRankCount =
metaData_->entity_rank_count();
1011 const size_t commCount = 10;
1016 stk::ParallelMachine mach = *
mpiComm_->getRawMpiComm();
1018 std::vector<stk::mesh::EntityId> local(commCount,0);
1021 stk::mesh::Selector ownedPart =
metaData_->locally_owned_part();
1022 for(stk::mesh::EntityRank i=stk::topology::NODE_RANK;
1023 i < static_cast<stk::mesh::EntityRank>(entityRankCount); ++i) {
1024 std::vector<stk::mesh::Entity> entities;
1026 stk::mesh::get_selected_entities(ownedPart,
bulkData_->buckets(i),entities);
1029 std::vector<stk::mesh::Entity>::const_iterator itr;
1030 for(itr=entities.begin();itr!=entities.end();++itr) {
1031 stk::mesh::EntityId
id =
bulkData_->identifier(*itr);
1038 stk::all_reduce(mach,stk::ReduceMax<10>(&local[0]));
1039 maxEntityId_.assign(local.begin(),local.begin()+entityRankCount+1);
1045 "STK_Interface::getEntityCounts: Entity counts do not include rank: " << entityRank);
1053 "STK_Interface::getMaxEntityId: Max entity ids do not include rank: " << entityRank);
1060 stk::mesh::PartVector emptyPartVector;
1061 stk::mesh::create_adjacent_entities(*
bulkData_,emptyPartVector);
1083 std::vector<stk::mesh::EntityId> & subcellIds)
const
1086 stk::mesh::Entity cell =
bulkData_->get_entity(elementRank,elementId);
1089 "STK_Interface::getSubcellIndices: could not find element requested (GID = " << elementId <<
")");
1091 const size_t numSubcells =
bulkData_->num_connectivity(cell, static_cast<stk::mesh::EntityRank>(entityRank));
1092 stk::mesh::Entity
const* subcells =
bulkData_->begin(cell, static_cast<stk::mesh::EntityRank>(entityRank));
1094 subcellIds.resize(numSubcells,0);
1097 for(
size_t i = 0; i < numSubcells; ++i) {
1098 stk::mesh::Entity subcell = subcells[i];
1099 subcellIds[i] =
bulkData_->identifier(subcell);
1106 stk::mesh::Selector ownedPart =
metaData_->locally_owned_part();
1110 stk::mesh::get_selected_entities(ownedPart,
bulkData_->buckets(elementRank),elements);
1121 stk::mesh::Selector ownedBlock =
metaData_->locally_owned_part() & (*elementBlock);
1125 stk::mesh::get_selected_entities(ownedBlock,
bulkData_->buckets(elementRank),elements);
1131 stk::mesh::Selector neighborBlock = (!
metaData_->locally_owned_part());
1135 stk::mesh::get_selected_entities(neighborBlock,
bulkData_->buckets(elementRank),elements);
1145 stk::mesh::Selector neighborBlock = (!
metaData_->locally_owned_part()) & (*elementBlock);
1149 stk::mesh::get_selected_entities(neighborBlock,
bulkData_->buckets(elementRank),elements);
1155 stk::mesh::Selector ownedPart =
metaData_->locally_owned_part();
1163 stk::mesh::Part * edgeBlockPart =
getEdgeBlock(edgeBlockName);
1165 "Unknown edge block \"" << edgeBlockName <<
"\"");
1167 stk::mesh::Selector edge_block = *edgeBlockPart;
1168 stk::mesh::Selector owned_block =
metaData_->locally_owned_part() & edge_block;
1174 void STK_Interface::getMyEdges(
const std::string & edgeBlockName,
const std::string & blockName,std::vector<stk::mesh::Entity> & edges)
const
1176 stk::mesh::Part * edgeBlockPart =
getEdgeBlock(edgeBlockName);
1179 "Unknown edge block \"" << edgeBlockName <<
"\"");
1181 "Unknown element block \"" << blockName <<
"\"");
1183 stk::mesh::Selector edge_block = *edgeBlockPart;
1184 stk::mesh::Selector element_block = *elmtPart;
1185 stk::mesh::Selector owned_block =
metaData_->locally_owned_part() & element_block & edge_block;
1193 stk::mesh::Part * edgeBlockPart =
getEdgeBlock(edgeBlockName);
1195 "Unknown edge block \"" << edgeBlockName <<
"\"");
1197 stk::mesh::Selector edge_block = *edgeBlockPart;
1205 stk::mesh::Part * edgeBlockPart =
getEdgeBlock(edgeBlockName);
1208 "Unknown edge block \"" << edgeBlockName <<
"\"");
1210 "Unknown element block \"" << blockName <<
"\"");
1212 stk::mesh::Selector edge_block = *edgeBlockPart;
1213 stk::mesh::Selector element_block = *elmtPart;
1214 stk::mesh::Selector element_edge_block = element_block & edge_block;
1223 stk::mesh::Selector ownedPart =
metaData_->locally_owned_part();
1227 stk::mesh::get_selected_entities(ownedPart,
bulkData_->buckets(faceRank),faces);
1232 stk::mesh::Part * faceBlockPart =
getFaceBlock(faceBlockName);
1234 "Unknown face block \"" << faceBlockName <<
"\"");
1236 stk::mesh::Selector face_block = *faceBlockPart;
1237 stk::mesh::Selector owned_block =
metaData_->locally_owned_part() & face_block;
1243 void STK_Interface::getMyFaces(
const std::string & faceBlockName,
const std::string & blockName,std::vector<stk::mesh::Entity> & faces)
const
1245 stk::mesh::Part * faceBlockPart =
getFaceBlock(faceBlockName);
1248 "Unknown face block \"" << faceBlockName <<
"\"");
1250 "Unknown element block \"" << blockName <<
"\"");
1252 stk::mesh::Selector face_block = *faceBlockPart;
1253 stk::mesh::Selector element_block = *elmtPart;
1254 stk::mesh::Selector owned_block =
metaData_->locally_owned_part() & element_block & face_block;
1262 stk::mesh::Part * faceBlockPart =
getFaceBlock(faceBlockName);
1264 "Unknown face block \"" << faceBlockName <<
"\"");
1266 stk::mesh::Selector face_block = *faceBlockPart;
1274 stk::mesh::Part * faceBlockPart =
getFaceBlock(faceBlockName);
1277 "Unknown face block \"" << faceBlockName <<
"\"");
1279 "Unknown element block \"" << blockName <<
"\"");
1281 stk::mesh::Selector face_block = *faceBlockPart;
1282 stk::mesh::Selector element_block = *elmtPart;
1283 stk::mesh::Selector element_face_block = element_block & face_block;
1291 stk::mesh::Part * sidePart =
getSideset(sideName);
1293 "Unknown side set \"" << sideName <<
"\"");
1295 stk::mesh::Selector side = *sidePart;
1296 stk::mesh::Selector ownedBlock =
metaData_->locally_owned_part() & side;
1304 stk::mesh::Part * sidePart =
getSideset(sideName);
1307 "Unknown side set \"" << sideName <<
"\"");
1309 "Unknown element block \"" << blockName <<
"\"");
1311 stk::mesh::Selector side = *sidePart;
1312 stk::mesh::Selector block = *elmtPart;
1313 stk::mesh::Selector ownedBlock =
metaData_->locally_owned_part() & block & side;
1321 stk::mesh::Part * sidePart =
getSideset(sideName);
1323 "Unknown side set \"" << sideName <<
"\"");
1325 stk::mesh::Selector side = *sidePart;
1333 stk::mesh::Part * sidePart =
getSideset(sideName);
1336 "Unknown side set \"" << sideName <<
"\"");
1338 "Unknown element block \"" << blockName <<
"\"");
1340 stk::mesh::Selector side = *sidePart;
1341 stk::mesh::Selector block = *elmtPart;
1342 stk::mesh::Selector sideBlock = block & side;
1350 stk::mesh::Part * nodePart =
getNodeset(nodesetName);
1353 "Unknown node set \"" << nodesetName <<
"\"");
1355 "Unknown element block \"" << blockName <<
"\"");
1357 stk::mesh::Selector nodeset = *nodePart;
1358 stk::mesh::Selector block = *elmtPart;
1359 stk::mesh::Selector ownedBlock =
metaData_->locally_owned_part() & block & nodeset;
1372 std::map<std::string, stk::mesh::Part*>::const_iterator blkItr;
1374 names.push_back(blkItr->first);
1384 std::map<std::string, stk::mesh::Part*>::const_iterator sideItr;
1386 names.push_back(sideItr->first);
1394 std::map<std::string, stk::mesh::Part*>::const_iterator nodeItr;
1396 names.push_back(nodeItr->first);
1404 std::map<std::string, stk::mesh::Part*>::const_iterator edgeBlockItr;
1406 names.push_back(edgeBlockItr->first);
1414 std::map<std::string, stk::mesh::Part*>::const_iterator faceBlockItr;
1416 names.push_back(faceBlockItr->first);
1432 std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr =
localIDHash_.find(gid);
1444 std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr =
localEdgeIDHash_.find(gid);
1458 std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr =
localEdgeIDHash_.find(gid);
1470 std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr =
localFaceIDHash_.find(gid);
1484 std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr =
localFaceIDHash_.find(gid);
1493 if(
bulkData_->bucket(elmt).member(*(eb_pair.second)))
1494 return eb_pair.first;
1499 const std::string & blockId)
const
1502 std::map<std::pair<std::string,std::string>,
SolutionFieldType*>::const_iterator
1507 "Solution field name \"" << fieldName <<
"\" in block ID \"" << blockId <<
"\" was not found");
1509 return iter->second;
1513 const std::string & blockId)
const
1516 std::map<std::pair<std::string,std::string>,
SolutionFieldType*>::const_iterator
1521 "Cell field named \"" << fieldName <<
"\" in block ID \"" << blockId <<
"\" was not found");
1523 return iter->second;
1527 const std::string & blockId)
const
1530 std::map<std::pair<std::string,std::string>,
SolutionFieldType*>::const_iterator
1535 "Edge field named \"" << fieldName <<
"\" in block ID \"" << blockId <<
"\" was not found");
1537 return iter->second;
1541 const std::string & blockId)
const
1544 std::map<std::pair<std::string,std::string>,
SolutionFieldType*>::const_iterator
1549 "Face field named \"" << fieldName <<
"\" in block ID \"" << blockId <<
"\" was not found");
1551 return iter->second;
1571 stk::mesh::Part * block =
metaData_->get_part(name);
1573 block = &
metaData_->declare_part_with_topology(name, stk::mesh::get_topology(shards::CellTopology(ctData),
dimension_));
1599 const std::string & edgeBlockName,
1600 const stk::topology & topology)
1604 stk::mesh::Part * edge_block =
metaData_->get_part(edgeBlockName);
1606 edge_block = &
metaData_->declare_part_with_topology(edgeBlockName, topology);
1612 stk::mesh::Part * elem_block =
metaData_->get_part(elemBlockName);
1613 metaData_->declare_part_subset(*elem_block, *edge_block);
1616 edgeBlocks_.insert(std::make_pair(edgeBlockName,edge_block));
1620 const std::string & edgeBlockName,
1621 const CellTopologyData * ctData)
1627 stk::mesh::get_topology(shards::CellTopology(ctData),
dimension_));
1644 const std::string & faceBlockName,
1645 const stk::topology & topology)
1649 stk::mesh::Part * face_block =
metaData_->get_part(faceBlockName);
1651 face_block = &
metaData_->declare_part_with_topology(faceBlockName, topology);
1657 stk::mesh::Part * elem_block =
metaData_->get_part(elemBlockName);
1658 metaData_->declare_part_subset(*elem_block, *face_block);
1661 faceBlocks_.insert(std::make_pair(faceBlockName,face_block));
1665 const std::string & faceBlockName,
1666 const CellTopologyData * ctData)
1672 stk::mesh::get_topology(shards::CellTopology(ctData),
dimension_));
1706 std::vector<stk::mesh::Entity> elements;
1709 for(std::size_t index=0;index<elements.size();++index) {
1710 stk::mesh::Entity element = elements[index];
1714 procId[0] = Teuchos::as<ProcIdData>(
procRank_);
1727 for(std::size_t index=0;index<elements.size();++index) {
1728 stk::mesh::Entity element = elements[index];
1732 procId[0] = Teuchos::as<ProcIdData>(
procRank_);
1744 std::vector<std::string> names;
1747 for(std::size_t b=0;b<names.size();b++) {
1749 std::map<std::string,double>::const_iterator bw_itr =
blockWeights_.find(names[b]);
1750 double blockWeight = (bw_itr!=
blockWeights_.end()) ? bw_itr->second : 1.0;
1752 std::vector<stk::mesh::Entity> elements;
1755 for(std::size_t index=0;index<elements.size();++index) {
1757 double * loadBal = stk::mesh::field_data(*
loadBalField_,elements[index]);
1758 loadBal[0] = blockWeight;
1770 std::vector<stk::mesh::Entity> edges;
1773 for(std::size_t index=0;index<edges.size();++index) {
1774 stk::mesh::Entity edge = edges[index];
1790 std::vector<stk::mesh::Entity> faces;
1793 for(std::size_t index=0;index<faces.size();++index) {
1794 stk::mesh::Entity face = faces[index];
1805 const std::string & fieldName,
1808 std::map<std::string,std::vector<std::string> >::const_iterator blkItr =
meshCoordFields_.find(eBlock);
1814 for(axis=0;axis<Teuchos::as<int>(blkItr->second.size());axis++) {
1815 if(blkItr->second[axis]==fieldName)
1819 if(axis>=Teuchos::as<int>(blkItr->second.size()))
1830 const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & matchers =
getPeriodicBCVector();
1832 std::vector<std::vector<std::string> > matchedSides(3);
1835 for(std::size_t m=0;m<matchers.size();m++){
1837 if(matchers[m]->getType() ==
"coord")
1839 else if(matchers[m]->getType() ==
"edge")
1841 else if(matchers[m]->getType() ==
"face")
1845 #ifdef PANZER_HAVE_STKSEARCH
1847 if (useBBoxSearch) {
1848 vec = matchers[m]->getMatchedPair(*
this,matchedSides[type],vec);
1850 vec = matchers[m]->getMatchedPair(*
this,vec);
1854 "panzer::STK_Interface::getPeriodicNodePairing(): Requested bounding box search, but "
1855 "did not compile with STK_SEARCH enabled.");
1856 vec = matchers[m]->getMatchedPair(*
this,vec);
1858 type_vec->insert(type_vec->begin(),vec->size()-type_vec->size(),type);
1859 matchedSides[type].push_back(matchers[m]->getLeftSidesetName());
1862 return std::make_pair(vec,type_vec);
1867 std::map<std::string, stk::mesh::Part*>::const_iterator blkItr =
elementBlocks_.find(blockId);
1874 std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1880 os <<
"STK Mesh data:\n";
1883 os <<
" Entity counts (Nodes, Edges, Cells) = ( "
1888 os <<
" Entity counts (Nodes, Edges, Faces, Cells) = ( "
1894 os <<
" Entity counts (Nodes, Cells) = ( "
1898 os <<
" Element blocks = ";
1899 for(std::size_t i=0;i<blockNames.size();i++)
1900 os <<
"\"" << blockNames[i] <<
"\" ";
1902 os <<
" Sidesets = ";
1903 for(std::size_t i=0;i<sidesetNames.size();i++)
1904 os <<
"\"" << sidesetNames[i] <<
"\" ";
1906 os <<
" Nodesets = ";
1907 for(std::size_t i=0;i<nodesetNames.size();i++)
1908 os <<
"\"" << nodesetNames[i] <<
"\" ";
1912 const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1914 if(bcVector.size()!=0) {
1915 os <<
" Periodic BCs:\n";
1916 for(std::size_t i=0;i<bcVector.size();i++)
1917 os <<
" " << bcVector[i]->getString() <<
"\n";
1924 std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1930 os <<
"STK Meta data:\n";
1931 os <<
" Element blocks = ";
1932 for(std::size_t i=0;i<blockNames.size();i++)
1933 os <<
"\"" << blockNames[i] <<
"\" ";
1935 os <<
" Sidesets = ";
1936 for(std::size_t i=0;i<sidesetNames.size();i++)
1937 os <<
"\"" << sidesetNames[i] <<
"\" ";
1939 os <<
" Nodesets = ";
1940 for(std::size_t i=0;i<nodesetNames.size();i++)
1941 os <<
"\"" << nodesetNames[i] <<
"\" ";
1945 const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1947 if(bcVector.size()!=0) {
1948 os <<
" Periodic BCs:\n";
1949 for(std::size_t i=0;i<bcVector.size();i++)
1950 os <<
" " << bcVector[i]->getString() <<
"\n";
1956 const stk::mesh::FieldVector & fv =
metaData_->get_fields();
1957 for(std::size_t i=0;i<fv.size();i++)
1958 os <<
"\"" << fv[i]->name() <<
"\" ";
1964 std::map<std::string, Teuchos::RCP<shards::CellTopology> >::const_iterator itr;
1968 std::stringstream ss;
1971 "STK_Interface::getCellTopology: No such element block \"" +eBlock +
"\" available.\n\n"
1972 <<
"STK Meta Data follows: \n" << ss.str());
1981 const int err = MPI_Comm_dup (parallelMach, &newComm);
1983 "panzer::STK_Interface: MPI_Comm_dup failed with error \""
1984 << Teuchos::mpiErrorCodeToString (err) <<
"\".");
1996 stk::mesh::Selector selector(
getMetaData()->universal_part());
1997 stk::mesh::Selector owned_selector(
getMetaData()->locally_owned_part());
2005 if(params.begin()!=params.end())
2006 graph.
sublist(stk::rebalance::Zoltan::default_parameters_name()) = params;
2025 #ifdef PANZER_HAVE_PERCEPT
2027 "ERROR: Number of levels for uniform refinement must be greater than 0");
2031 refinedMesh_->setCoordinatesField();
2033 percept::UniformRefiner breaker(*refinedMesh_,*breakPattern_);
2034 breaker.setRemoveOldElements(deleteParentElements);
2036 for (
int i=0; i < numberOfLevels; ++i)
2041 "ERROR: Uniform refinement requested. This requires the Percept package to be enabled in Trilinos!");
void initializeFromMetaData()
Teuchos::RCP< ElementDescriptor > buildElementDescriptor(stk::mesh::EntityId elmtId, std::vector< stk::mesh::EntityId > &nodes)
Teuchos::RCP< Teuchos::MpiComm< int > > getSafeCommunicator(stk::ParallelMachine parallelMach) const
void getSidesetNames(std::vector< std::string > &name) const
RCP< const T > getConst() const
std::pair< Teuchos::RCP< std::vector< std::pair< std::size_t, std::size_t > > >, Teuchos::RCP< std::vector< unsigned int > > > getPeriodicNodePairing() const
const bool & useBoundingBoxSearch() const
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedEdgeVector_
void addNodeset(const std::string &name)
Teuchos::RCP< stk::mesh::BulkData > bulkData_
void getFaceBlockNames(std::vector< std::string > &names) const
void print(std::ostream &os) const
void getMyEdges(std::vector< stk::mesh::Entity > &edges) const
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getElementsOrderedByLID() const
std::vector< stk::mesh::Part * > facesPartVec_
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToSolution_
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void addGlobalToExodus(const std::string &key, const int &value)
Add an int global variable to the information to be written to the Exodus output file.
void getElementBlockNames(std::vector< std::string > &names) const
void getSubcellIndices(unsigned entityRank, stk::mesh::EntityId elementId, std::vector< stk::mesh::EntityId > &subcellIds) const
void getElementsSharingNodes(const std::vector< stk::mesh::EntityId > nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing multiple nodes
stk::mesh::Field< double > SolutionFieldType
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
std::vector< stk::mesh::Part * > edgesPartVec_
bool isFaceLocal(stk::mesh::Entity face) const
std::map< std::string, double > blockWeights_
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getEdgesOrderedByLID() const
VectorFieldType * coordinatesField_
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::size_t getEntityCounts(unsigned entityRank) const
get the global counts for the entity of specified rank
void printMetaData(std::ostream &os) const
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
std::map< std::string, stk::mesh::Part * > faceBlocks_
stk::mesh::Part * nodesPart_
void addSolutionField(const std::string &fieldName, const std::string &blockId)
std::map< std::string, Teuchos::RCP< shards::CellTopology > > elementBlockCT_
void addEntityToFaceBlock(stk::mesh::Entity entity, stk::mesh::Part *faceblock)
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
std::size_t currentLocalId_
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToFaceField_
std::map< std::string, stk::mesh::Part * > sidesets_
std::vector< std::size_t > entityCounts_
bool validBlockId(const std::string &blockId) const
stk::mesh::Part * getNodeset(const std::string &name) const
void getEdgeBlockNames(std::vector< std::string > &names) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
void addInformationRecords(const std::vector< std::string > &info_records)
void endModification(const bool find_and_set_shared_nodes_in_stk=true)
void getAllFaces(const std::string &faceBlockName, std::vector< stk::mesh::Entity > &faces) const
stk::mesh::Field< double > * getSolutionField(const std::string &fieldName, const std::string &blockId) const
SolutionFieldType * loadBalField_
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const
stk::mesh::Part * getSideset(const std::string &name) const
static const std::string nodesString
void instantiateBulkData(stk::ParallelMachine parallelMach)
std::set< std::string > informationRecords_
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
stk::mesh::Field< double > * getEdgeField(const std::string &fieldName, const std::string &blockId) const
std::size_t elementLocalId(stk::mesh::Entity elmt) const
unsigned getDimension() const
get the dimension
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
std::map< std::string, std::vector< std::string > > meshDispFields_
void setupExodusFile(const std::string &filename, const bool append=false, const bool append_after_restart_time=false, const double restart_time=0.0)
Set up an output Exodus file for writing results.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
VectorFieldType * edgesField_
stk::mesh::EntityRank getSideRank() const
void getAllEdges(const std::string &edgeBlockName, std::vector< stk::mesh::Entity > &edges) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
static const std::string coordsString
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
std::map< std::string, stk::mesh::Part * > edgeBlocks_
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
std::map< std::string, stk::mesh::Part * > elementBlocks_
std::vector< stk::mesh::Part * > nodesPartVec_
static const std::string edgeBlockString
static const std::string edgesString
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
std::string containingBlockId(stk::mesh::Entity elmt) const
stk::mesh::EntityRank getFaceRank() const
std::size_t edgeLocalId(stk::mesh::Entity elmt) const
void getOwnedElementsSharingNode(stk::mesh::Entity node, std::vector< stk::mesh::Entity > &elements, std::vector< int > &relIds) const
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void buildSubcells()
force the mesh to build subcells: edges and faces
Teuchos::RCP< Teuchos::MpiComm< int > > mpiComm_
Teuchos::RCP< stk::mesh::MetaData > metaData_
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
virtual ~ElementDescriptor()
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedElementVector_
bool isInitialized() const
Has initialize been called on this mesh object?
void writeToExodus(const std::string &filename, const bool append=false)
Write this mesh and associated fields to the given output file.
void getMyNodes(const std::string &sideName, const std::string &blockName, std::vector< stk::mesh::Entity > &nodes) const
bool nonnull(const boost::shared_ptr< T > &p)
void addEdgeField(const std::string &fieldName, const std::string &blockId)
void addMeshCoordFields(const std::string &blockId, const std::vector< std::string > &coordField, const std::string &dispPrefix)
void applyElementLoadBalanceWeights()
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
bool isMeshCoordField(const std::string &eBlock, const std::string &fieldName, int &axis) const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
static const std::string faceBlockString
std::unordered_map< stk::mesh::EntityId, std::size_t > localFaceIDHash_
void addEntityToEdgeBlock(stk::mesh::Entity entity, stk::mesh::Part *edgeblock)
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToEdgeField_
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::unordered_map< stk::mesh::EntityId, std::size_t > localEdgeIDHash_
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
void rebalance(const Teuchos::ParameterList ¶ms)
std::map< std::string, std::vector< std::string > > meshCoordFields_
stk::mesh::EntityRank getEdgeRank() const
static const std::string facesString
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
void initializeFieldsInSTK(const std::map< std::pair< std::string, std::string >, SolutionFieldType * > &nameToField, bool setupIO)
const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector() const
void refineMesh(const int numberOfLevels, const bool deleteParentElements)
stk::mesh::Field< double > * getFaceField(const std::string &fieldName, const std::string &blockId) const
void buildLocalElementIDs()
std::size_t faceLocalId(stk::mesh::Entity elmt) const
bool isModifiable() const
void addFaceBlock(const std::string &elemBlockName, const std::string &faceBlockName, const stk::topology &topology)
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedFaceVector_
std::unordered_map< stk::mesh::EntityId, std::size_t > localIDHash_
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
std::vector< stk::mesh::EntityId > maxEntityId_
void getElementsSharingNode(stk::mesh::EntityId nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing a single node
void addFaceField(const std::string &fieldName, const std::string &blockId)
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
#define TEUCHOS_ASSERT(assertion_test)
stk::mesh::EntityRank getNodeRank() const
void getMyFaces(std::vector< stk::mesh::Entity > &faces) const
void addCellField(const std::string &fieldName, const std::string &blockId)
void getNodeIdsForElement(stk::mesh::Entity element, std::vector< stk::mesh::EntityId > &nodeIds) const
get a list of node ids for nodes connected to an element
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getFacesOrderedByLID() const
stk::mesh::Part * edgesPart_
stk::mesh::EntityRank getElementRank() const
stk::mesh::Part * facesPart_
ProcIdFieldType * processorIdField_
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
const VectorFieldType & getCoordinatesField() const
bool isEdgeLocal(stk::mesh::Entity edge) const
void getNeighborElements(std::vector< stk::mesh::Entity > &elements) const
void getNodesetNames(std::vector< std::string > &name) const
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToCellField_
std::map< std::string, stk::mesh::Part * > nodesets_
VectorFieldType * facesField_