13 #include <Teuchos_RCPStdSharedPtrConversions.hpp>
14 #include <PanzerAdaptersSTK_config.hpp>
19 #ifdef PANZER_HAVE_IOSS
21 #include <Ionit_Initializer.h>
22 #include <Ioss_ElementBlock.h>
23 #include <Ioss_EdgeBlock.h>
24 #include <Ioss_FaceBlock.h>
25 #include <Ioss_Region.h>
26 #include <stk_io/StkMeshIoBroker.hpp>
27 #include <stk_io/IossBridge.hpp>
28 #include <stk_mesh/base/FieldParallel.hpp>
30 #ifdef PANZER_HAVE_UMR
31 #include <Ioumr_DatabaseIO.hpp>
34 #include "Teuchos_StandardParameterEntryValidators.hpp"
36 namespace panzer_stk {
38 int getMeshDimension(
const std::string & meshStr,
39 stk::ParallelMachine parallelMach,
40 const std::string & typeStr)
42 stk::io::StkMeshIoBroker meshData(parallelMach);
43 meshData.use_simple_fields();
44 meshData.property_add(Ioss::Property(
"LOWER_CASE_VARIABLE_NAMES",
false));
45 meshData.add_mesh_database(meshStr, fileTypeToIOSSType(typeStr), stk::io::READ_MESH);
46 meshData.create_input_mesh();
47 return Teuchos::as<int>(meshData.meta_data_ptr()->spatial_dimension());
50 std::string fileTypeToIOSSType(
const std::string & fileType)
53 if (fileType==
"Exodus")
54 IOSSType =
"exodusii";
55 #ifdef PANZER_HAVE_UMR
56 else if (fileType==
"Exodus Refinement")
57 IOSSType =
"Refinement";
59 else if (fileType==
"Pamgen")
65 STK_ExodusReaderFactory::STK_ExodusReaderFactory()
66 : fileName_(
""), fileType_(
""), restartIndex_(0), userMeshScaling_(false), keepPerceptData_(false),
67 keepPerceptParentElements_(false), rebalancing_(
"default"),
68 meshScaleFactor_(0.0), levelsOfRefinement_(0),
69 createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_(
"")
72 STK_ExodusReaderFactory::STK_ExodusReaderFactory(
const std::string & fileName,
73 const int restartIndex)
74 : fileName_(fileName), fileType_(
"Exodus"), restartIndex_(restartIndex), userMeshScaling_(false),
75 keepPerceptData_(false), keepPerceptParentElements_(false), rebalancing_(
"default"),
76 meshScaleFactor_(0.0), levelsOfRefinement_(0), createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_(
"")
81 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::buildMesh()");
86 RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
92 mesh->initialize(parallelMach,
false,doPerceptRefinement());
94 completeMeshConstruction(*mesh,parallelMach);
105 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::buildUncomittedMesh()");
111 stk::io::StkMeshIoBroker* meshData =
new stk::io::StkMeshIoBroker(parallelMach);
112 meshData->use_simple_fields();
113 meshData->property_add(Ioss::Property(
"LOWER_CASE_VARIABLE_NAMES",
false));
116 std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
117 entity_rank_names.push_back(
"FAMILY_TREE");
118 meshData->set_rank_name_vector(entity_rank_names);
120 #ifdef PANZER_HAVE_UMR
122 Ioumr::IOFactory::factory();
124 meshData->property_add(Ioss::Property(
"GEOMETRY_FILE", geometryName_));
125 meshData->property_add(Ioss::Property(
"NUMBER_REFINEMENTS", levelsOfRefinement_));
128 meshData->add_mesh_database(fileName_, fileTypeToIOSSType(fileType_), stk::io::READ_MESH);
130 meshData->create_input_mesh();
131 RCP<stk::mesh::MetaData> metaData =
Teuchos::rcp(meshData->meta_data_ptr());
133 RCP<STK_Interface> mesh =
rcp(
new STK_Interface(metaData));
134 mesh->initializeFromMetaData();
135 mesh->instantiateBulkData(parallelMach);
136 meshData->set_bulk_data(Teuchos::get_shared_ptr(mesh->getBulkData()));
140 meshData->add_all_mesh_fields_as_input_fields();
144 mesh->getMetaData()->declare_attribute_with_delete(meshData);
147 registerElementBlocks(*mesh,*meshData);
148 registerSidesets(*mesh);
149 registerNodesets(*mesh);
151 if (createEdgeBlocks_) {
152 registerEdgeBlocks(*mesh,*meshData);
154 if (createFaceBlocks_ && mesh->getMetaData()->spatial_dimension() > 2) {
155 registerFaceBlocks(*mesh,*meshData);
158 buildMetaData(parallelMach, *mesh);
160 mesh->addPeriodicBCs(periodicBCVec_);
161 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
166 void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach)
const
168 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::completeMeshConstruction()");
174 if(not mesh.isInitialized()) {
175 mesh.initialize(parallelMach,
true,doPerceptRefinement());
179 stk::mesh::MetaData & metaData = *mesh.getMetaData();
180 stk::mesh::BulkData & bulkData = *mesh.getBulkData();
181 stk::io::StkMeshIoBroker * meshData =
182 const_cast<stk::io::StkMeshIoBroker *
>(metaData.get_attribute<stk::io::StkMeshIoBroker>());
190 meshData->populate_bulk_data();
192 if (doPerceptRefinement()) {
193 const bool deleteParentElements = !keepPerceptParentElements_;
194 mesh.refineMesh(levelsOfRefinement_,deleteParentElements);
199 if (userMeshScaling_)
201 stk::mesh::Field<double>* coord_field = metaData.get_field<
double>(stk::topology::NODE_RANK,
"coordinates");
203 stk::mesh::Selector select_all_local = metaData.locally_owned_part() | metaData.globally_shared_part();
204 stk::mesh::BucketVector
const& my_node_buckets = bulkData.get_buckets(stk::topology::NODE_RANK, select_all_local);
206 int mesh_dim = mesh.getDimension();
209 const double inv_msf = 1.0/meshScaleFactor_;
210 for (
size_t i=0; i < my_node_buckets.size(); ++i)
212 stk::mesh::Bucket& b = *(my_node_buckets[i]);
213 double* coordinate_data = field_data( *coord_field, b );
215 for (
size_t j=0; j < b.size(); ++j) {
216 for (
int k=0; k < mesh_dim; ++k) {
217 coordinate_data[mesh_dim*j + k] *= inv_msf;
225 int restartIndex = restartIndex_;
227 std::pair<int,double> lastTimeStep = meshData->get_input_ioss_region()->get_max_time();
228 restartIndex = 1+restartIndex+lastTimeStep.first;
232 meshData->read_defined_input_fields(restartIndex);
234 mesh.buildSubcells();
235 mesh.buildLocalElementIDs();
237 mesh.beginModification();
238 if (createEdgeBlocks_) {
239 mesh.buildLocalEdgeIDs();
242 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
243 mesh.buildLocalFaceIDs();
248 const bool find_and_set_shared_nodes_in_stk =
false;
249 mesh.endModification(find_and_set_shared_nodes_in_stk);
252 if (userMeshScaling_) {
253 stk::mesh::Field<double>* coord_field = metaData.get_field<
double>(stk::topology::NODE_RANK,
"coordinates");
254 std::vector< const stk::mesh::FieldBase *> fields;
255 fields.push_back(coord_field);
257 stk::mesh::communicate_field_data(bulkData, fields);
261 mesh.setInitialStateTime(meshData->get_input_ioss_region()->get_state_time(restartIndex));
263 mesh.setInitialStateTime(0.0);
268 if(rebalancing_ ==
"default")
270 this->rebalance(mesh);
271 else if(rebalancing_ !=
"none")
274 "ERROR: Rebalancing was not set to a valid choice");
283 "Error, the parameter {name=\"File Name\","
285 "\nis required in parameter (sub)list \""<< paramList->
name() <<
"\"."
293 paramList->
set<
int>(
"Restart Index", -1);
296 paramList->
set(
"File Type",
"Exodus");
298 if(!paramList->
isSublist(
"Periodic BCs"))
299 paramList->
sublist(
"Periodic BCs");
303 p_bcs.
set<
int>(
"Count", 0);
305 if(!paramList->
isParameter(
"Levels of Uniform Refinement"))
306 paramList->
set<
int>(
"Levels of Uniform Refinement", 0);
309 paramList->
set<
bool>(
"Keep Percept Data",
false);
311 if(!paramList->
isParameter(
"Keep Percept Parent Elements"))
312 paramList->
set<
bool>(
"Keep Percept Parent Elements",
false);
315 paramList->
set<std::string>(
"Rebalancing",
"default");
319 paramList->
set<
bool>(
"Create Edge Blocks",
false);
323 paramList->
set<
bool>(
"Create Face Blocks",
false);
326 paramList->
set(
"Geometry File Name",
"");
330 setMyParamList(paramList);
332 fileName_ = paramList->
get<std::string>(
"File Name");
334 geometryName_ = paramList->
get<std::string>(
"Geometry File Name");
336 restartIndex_ = paramList->
get<
int>(
"Restart Index");
338 fileType_ = paramList->
get<std::string>(
"File Type");
343 meshScaleFactor_ = paramList->
get<
double>(
"Scale Factor");
344 userMeshScaling_ =
true;
348 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->
sublist(
"Periodic BCs")),periodicBCVec_,useBBoxSearch_);
350 keepPerceptData_ = paramList->
get<
bool>(
"Keep Percept Data");
352 keepPerceptParentElements_ = paramList->
get<
bool>(
"Keep Percept Parent Elements");
354 rebalancing_ = paramList->
get<std::string>(
"Rebalancing");
356 levelsOfRefinement_ = paramList->
get<
int>(
"Levels of Uniform Refinement");
358 createEdgeBlocks_ = paramList->
get<
bool>(
"Create Edge Blocks");
359 createFaceBlocks_ = paramList->
get<
bool>(
"Create Face Blocks");
367 if(validParams==Teuchos::null) {
369 validParams->
set<std::string>(
"File Name",
"<file name not set>",
"Name of exodus file to be read",
371 validParams->
set<std::string>(
"Geometry File Name",
"<file name not set>",
"Name of geometry file for refinement",
374 validParams->
set<
int>(
"Restart Index",-1,
"Index of solution to read in",
377 validParams->
set<std::string>(
"File Type",
"Exodus",
378 "Choose input file type - either \"Exodus\", \"Exodus Refinement\" or \"Pamgen\"",
380 #ifdef PANZER_HAVE_UMR
385 validParams->
set<
double>(
"Scale Factor", 1.0,
"Scale factor to apply to mesh after read",
389 bcs.
set<
int>(
"Count",0);
391 validParams->
set(
"Levels of Uniform Refinement",0,
"Number of levels of inline uniform mesh refinement");
393 validParams->
set(
"Keep Percept Data",
false,
"Keep the Percept mesh after uniform refinement is applied");
395 validParams->
set(
"Keep Percept Parent Elements",
false,
"Keep the parent element information in the Percept data");
397 validParams->
set(
"Rebalancing",
"default",
"The type of rebalancing to be performed on the mesh after creation (default, none)");
400 validParams->
set(
"Create Edge Blocks",
false,
"Create or copy edge blocks in the mesh");
401 validParams->
set(
"Create Face Blocks",
false,
"Create or copy face blocks in the mesh");
407 void STK_ExodusReaderFactory::registerElementBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
411 RCP<stk::mesh::MetaData> femMetaData = mesh.getMetaData();
416 const Ioss::ElementBlockContainer & elem_blocks = meshData.get_input_ioss_region()->get_element_blocks();
417 for(Ioss::ElementBlockContainer::const_iterator itr=elem_blocks.begin();itr!=elem_blocks.end();++itr) {
418 Ioss::GroupingEntity * entity = *itr;
419 const std::string & name = entity->name();
421 const stk::mesh::Part * part = femMetaData->get_part(name);
422 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(femMetaData->get_topology(*part));
423 const CellTopologyData * ct = cellTopo.getCellTopologyData();
426 mesh.addElementBlock(part->name(),ct);
428 if (createEdgeBlocks_) {
429 createUniqueEdgeTopologyMap(mesh, part);
431 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
432 createUniqueFaceTopologyMap(mesh, part);
437 template <
typename SetType>
438 void buildSetNames(
const SetType & setData,std::vector<std::string> & names)
441 for(
typename SetType::const_iterator itr=setData.begin();itr!=setData.end();++itr) {
442 Ioss::GroupingEntity * entity = *itr;
443 names.push_back(entity->name());
447 void STK_ExodusReaderFactory::registerSidesets(STK_Interface & mesh)
const
451 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
452 const stk::mesh::PartVector & parts = metaData->get_parts();
454 stk::mesh::PartVector::const_iterator partItr;
455 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
456 const stk::mesh::Part * part = *partItr;
457 const stk::mesh::PartVector & subsets = part->subsets();
458 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
459 const CellTopologyData * ct = cellTopo.getCellTopologyData();
463 if(part->primary_entity_rank()==mesh.getSideRank() && ct==0 && subsets.size()>0) {
465 "STK_ExodusReaderFactory::registerSidesets error - part \"" << part->name() <<
466 "\" has more than one subset");
469 const stk::mesh::Part * ss_part = subsets[0];
470 shards::CellTopology ss_cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*ss_part));
471 const CellTopologyData * ss_ct = ss_cellTopo.getCellTopologyData();
475 mesh.addSideset(part->name(),ss_ct);
480 void STK_ExodusReaderFactory::registerNodesets(STK_Interface & mesh)
const
484 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
485 const stk::mesh::PartVector & parts = metaData->get_parts();
487 stk::mesh::PartVector::const_iterator partItr;
488 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
489 const stk::mesh::Part * part = *partItr;
490 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
491 const CellTopologyData * ct = cellTopo.getCellTopologyData();
495 if(part->primary_entity_rank()==mesh.getNodeRank() && ct==0) {
499 mesh.addNodeset(part->name());
504 void STK_ExodusReaderFactory::registerEdgeBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
508 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
516 const Ioss::EdgeBlockContainer & edge_blocks = meshData.get_input_ioss_region()->get_edge_blocks();
517 for(Ioss::EdgeBlockContainer::const_iterator ebc_iter=edge_blocks.begin();ebc_iter!=edge_blocks.end();++ebc_iter) {
518 Ioss::GroupingEntity * entity = *ebc_iter;
519 const stk::mesh::Part * edgeBlockPart = metaData->get_part(entity->name());
520 const stk::topology edgeBlockTopo = metaData->get_topology(*edgeBlockPart);
522 for (
auto ebuet_iter : elemBlockUniqueEdgeTopologies_) {
523 std::string elemBlockName = ebuet_iter.first;
524 std::vector<stk::topology> uniqueEdgeTopologies = ebuet_iter.second;
526 auto find_result = std::find(uniqueEdgeTopologies.begin(), uniqueEdgeTopologies.end(), edgeBlockTopo);
527 if (find_result != uniqueEdgeTopologies.end()) {
528 mesh.addEdgeBlock(elemBlockName, edgeBlockPart->name(), edgeBlockTopo);
534 void STK_ExodusReaderFactory::registerFaceBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
538 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
546 const Ioss::FaceBlockContainer & face_blocks = meshData.get_input_ioss_region()->get_face_blocks();
547 for(Ioss::FaceBlockContainer::const_iterator fbc_itr=face_blocks.begin();fbc_itr!=face_blocks.end();++fbc_itr) {
548 Ioss::GroupingEntity * entity = *fbc_itr;
549 const stk::mesh::Part * faceBlockPart = metaData->get_part(entity->name());
550 const stk::topology faceBlockTopo = metaData->get_topology(*faceBlockPart);
552 for (
auto ebuft_iter : elemBlockUniqueFaceTopologies_) {
553 std::string elemBlockName = ebuft_iter.first;
554 std::vector<stk::topology> uniqueFaceTopologies = ebuft_iter.second;
556 auto find_result = std::find(uniqueFaceTopologies.begin(), uniqueFaceTopologies.end(), faceBlockTopo);
557 if (find_result != uniqueFaceTopologies.end()) {
558 mesh.addFaceBlock(elemBlockName, faceBlockPart->name(), faceBlockTopo);
564 bool topo_less (stk::topology &i,stk::topology &j) {
return (i.value() < j.value()); }
565 bool topo_equal (stk::topology &i,stk::topology &j) {
return (i.value() == j.value()); }
567 void STK_ExodusReaderFactory::createUniqueEdgeTopologyMap(STK_Interface & mesh,
const stk::mesh::Part *elemBlockPart)
const
574 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
575 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
577 std::vector<stk::topology> edge_topologies;
578 for (
unsigned i=0;i<elemBlockTopo.num_edges();i++) {
579 edge_topologies.push_back(elemBlockTopo.edge_topology(i));
581 std::sort(edge_topologies.begin(), edge_topologies.end(), topo_less);
582 std::vector<stk::topology>::iterator new_end;
583 new_end = std::unique(edge_topologies.begin(), edge_topologies.end(), topo_equal);
584 edge_topologies.resize( std::distance(edge_topologies.begin(),new_end) );
586 elemBlockUniqueEdgeTopologies_[elemBlockPart->name()] = edge_topologies;
589 void STK_ExodusReaderFactory::createUniqueFaceTopologyMap(STK_Interface & mesh,
const stk::mesh::Part *elemBlockPart)
const
596 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
597 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
599 std::vector<stk::topology> face_topologies;
600 for (
unsigned i=0;i<elemBlockTopo.num_faces();i++) {
601 face_topologies.push_back(elemBlockTopo.face_topology(i));
603 std::sort(face_topologies.begin(), face_topologies.end(), topo_less);
604 std::vector<stk::topology>::iterator new_end;
605 new_end = std::unique(face_topologies.begin(), face_topologies.end(), topo_equal);
606 face_topologies.resize( std::distance(face_topologies.begin(),new_end) );
608 elemBlockUniqueFaceTopologies_[elemBlockPart->name()] = face_topologies;
613 void STK_ExodusReaderFactory::addEdgeBlocks(STK_Interface & mesh)
const
625 for (
auto iter : elemBlockUniqueEdgeTopologies_) {
626 std::string elemBlockName =
iter.first;
627 std::vector<stk::topology> uniqueEdgeTopologies =
iter.second;
629 for (
auto topo : uniqueEdgeTopologies ) {
630 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
631 const stk::mesh::Part & edgeTopoPart = metaData->get_topology_root_part(topo);
633 stk::mesh::Selector owned_block;
634 owned_block = *elemBlockPart;
635 owned_block &= edgeTopoPart;
638 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edge_block_name);
640 bulkData->change_entity_parts(owned_block, mesh.getEdgeRank(), stk::mesh::PartVector{edge_block});
647 void STK_ExodusReaderFactory::addFaceBlocks(STK_Interface & mesh)
const
659 for (
auto iter : elemBlockUniqueFaceTopologies_) {
660 std::string elemBlockName =
iter.first;
661 std::vector<stk::topology> uniqueFaceTopologies =
iter.second;
663 for (
auto topo : uniqueFaceTopologies ) {
664 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
665 const stk::mesh::Part & faceTopoPart = metaData->get_topology_root_part(topo);
667 stk::mesh::Selector owned_block;
668 owned_block = *elemBlockPart;
669 owned_block &= faceTopoPart;
672 stk::mesh::Part * face_block = mesh.getFaceBlock(face_block_name);
674 bulkData->change_entity_parts(owned_block, mesh.getFaceRank(), stk::mesh::PartVector{face_block});
679 void STK_ExodusReaderFactory::buildMetaData(stk::ParallelMachine , STK_Interface & mesh)
const
681 if (createEdgeBlocks_) {
685 for (
auto iter : elemBlockUniqueEdgeTopologies_) {
686 std::string elemBlockName =
iter.first;
687 std::vector<stk::topology> uniqueEdgeTopologies =
iter.second;
689 for (
auto topo : uniqueEdgeTopologies ) {
691 mesh.addEdgeBlock(elemBlockName, edge_block_name, topo);
695 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
699 for (
auto iter : elemBlockUniqueFaceTopologies_) {
700 std::string elemBlockName =
iter.first;
701 std::vector<stk::topology> uniqueFaceTopologies =
iter.second;
703 for (
auto topo : uniqueFaceTopologies ) {
705 mesh.addFaceBlock(elemBlockName, face_block_name, topo);
711 bool STK_ExodusReaderFactory::doPerceptRefinement()
const
713 return (fileType_!=
"Exodus Refinement") && (levelsOfRefinement_ > 0);
716 std::string STK_ExodusReaderFactory::mkBlockName(std::string base, std::string topo_name)
const
719 name = topo_name+
"_"+base;
720 std::transform(name.begin(), name.end(), name.begin(),
722 {
return char(std::tolower(c)); });
RCP< const T > getConst() const
const std::string & name() const
std::string currentParametersString() const
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static const std::string nodesString
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSublist(const std::string &name) const
static const std::string edgeBlockString
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
static const std::string faceBlockString
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(throw_exception_test, Exception, msg)