46 #include <Teuchos_RCPStdSharedPtrConversions.hpp>
47 #include <PanzerAdaptersSTK_config.hpp>
52 #ifdef PANZER_HAVE_IOSS
54 #include <Ionit_Initializer.h>
55 #include <Ioss_ElementBlock.h>
56 #include <Ioss_EdgeBlock.h>
57 #include <Ioss_FaceBlock.h>
58 #include <Ioss_Region.h>
59 #include <stk_mesh/base/GetBuckets.hpp>
60 #include <stk_io/StkMeshIoBroker.hpp>
61 #include <stk_io/IossBridge.hpp>
62 #include <stk_mesh/base/FieldParallel.hpp>
64 #ifdef PANZER_HAVE_UMR
65 #include <Ioumr_DatabaseIO.hpp>
68 #include "Teuchos_StandardParameterEntryValidators.hpp"
70 namespace panzer_stk {
72 int getMeshDimension(
const std::string & meshStr,
73 stk::ParallelMachine parallelMach,
74 const std::string & typeStr)
76 stk::io::StkMeshIoBroker meshData(parallelMach);
77 meshData.use_simple_fields();
78 meshData.property_add(Ioss::Property(
"LOWER_CASE_VARIABLE_NAMES",
false));
79 meshData.add_mesh_database(meshStr, fileTypeToIOSSType(typeStr), stk::io::READ_MESH);
80 meshData.create_input_mesh();
81 return Teuchos::as<int>(meshData.meta_data_ptr()->spatial_dimension());
84 std::string fileTypeToIOSSType(
const std::string & fileType)
87 if (fileType==
"Exodus")
88 IOSSType =
"exodusii";
89 #ifdef PANZER_HAVE_UMR
90 else if (fileType==
"Exodus Refinement")
91 IOSSType =
"Refinement";
93 else if (fileType==
"Pamgen")
99 STK_ExodusReaderFactory::STK_ExodusReaderFactory()
100 : fileName_(
""), fileType_(
""), restartIndex_(0), userMeshScaling_(false), keepPerceptData_(false),
101 keepPerceptParentElements_(false), rebalancing_(
"default"),
102 meshScaleFactor_(0.0), levelsOfRefinement_(0),
103 createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_(
"")
106 STK_ExodusReaderFactory::STK_ExodusReaderFactory(
const std::string & fileName,
107 const int restartIndex)
108 : fileName_(fileName), fileType_(
"Exodus"), restartIndex_(restartIndex), userMeshScaling_(false),
109 keepPerceptData_(false), keepPerceptParentElements_(false), rebalancing_(
"default"),
110 meshScaleFactor_(0.0), levelsOfRefinement_(0), createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_(
"")
115 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::buildMesh()");
120 RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
126 mesh->initialize(parallelMach,
false,doPerceptRefinement());
128 completeMeshConstruction(*mesh,parallelMach);
139 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::buildUncomittedMesh()");
145 stk::io::StkMeshIoBroker* meshData =
new stk::io::StkMeshIoBroker(parallelMach);
146 meshData->use_simple_fields();
147 meshData->property_add(Ioss::Property(
"LOWER_CASE_VARIABLE_NAMES",
false));
150 std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
151 entity_rank_names.push_back(
"FAMILY_TREE");
152 meshData->set_rank_name_vector(entity_rank_names);
154 #ifdef PANZER_HAVE_UMR
156 Ioumr::IOFactory::factory();
158 meshData->property_add(Ioss::Property(
"GEOMETRY_FILE", geometryName_));
159 meshData->property_add(Ioss::Property(
"NUMBER_REFINEMENTS", levelsOfRefinement_));
162 meshData->add_mesh_database(fileName_, fileTypeToIOSSType(fileType_), stk::io::READ_MESH);
164 meshData->create_input_mesh();
165 RCP<stk::mesh::MetaData> metaData =
Teuchos::rcp(meshData->meta_data_ptr());
167 RCP<STK_Interface> mesh =
rcp(
new STK_Interface(metaData));
168 mesh->initializeFromMetaData();
169 mesh->instantiateBulkData(parallelMach);
170 meshData->set_bulk_data(Teuchos::get_shared_ptr(mesh->getBulkData()));
174 meshData->add_all_mesh_fields_as_input_fields();
178 mesh->getMetaData()->declare_attribute_with_delete(meshData);
181 registerElementBlocks(*mesh,*meshData);
182 registerSidesets(*mesh);
183 registerNodesets(*mesh);
185 if (createEdgeBlocks_) {
186 registerEdgeBlocks(*mesh,*meshData);
188 if (createFaceBlocks_ && mesh->getMetaData()->spatial_dimension() > 2) {
189 registerFaceBlocks(*mesh,*meshData);
192 buildMetaData(parallelMach, *mesh);
194 mesh->addPeriodicBCs(periodicBCVec_);
195 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
200 void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach)
const
202 PANZER_FUNC_TIME_MONITOR(
"panzer::STK_ExodusReaderFactory::completeMeshConstruction()");
208 if(not mesh.isInitialized()) {
209 mesh.initialize(parallelMach,
true,doPerceptRefinement());
213 stk::mesh::MetaData & metaData = *mesh.getMetaData();
214 stk::mesh::BulkData & bulkData = *mesh.getBulkData();
215 stk::io::StkMeshIoBroker * meshData =
216 const_cast<stk::io::StkMeshIoBroker *
>(metaData.get_attribute<stk::io::StkMeshIoBroker>());
224 meshData->populate_bulk_data();
226 if (doPerceptRefinement()) {
227 const bool deleteParentElements = !keepPerceptParentElements_;
228 mesh.refineMesh(levelsOfRefinement_,deleteParentElements);
233 if (userMeshScaling_)
235 stk::mesh::Field<double>* coord_field = metaData.get_field<
double>(stk::topology::NODE_RANK,
"coordinates");
237 stk::mesh::Selector select_all_local = metaData.locally_owned_part() | metaData.globally_shared_part();
238 stk::mesh::BucketVector
const& my_node_buckets = bulkData.get_buckets(stk::topology::NODE_RANK, select_all_local);
240 int mesh_dim = mesh.getDimension();
243 const double inv_msf = 1.0/meshScaleFactor_;
244 for (
size_t i=0; i < my_node_buckets.size(); ++i)
246 stk::mesh::Bucket& b = *(my_node_buckets[i]);
247 double* coordinate_data = field_data( *coord_field, b );
249 for (
size_t j=0; j < b.size(); ++j) {
250 for (
int k=0; k < mesh_dim; ++k) {
251 coordinate_data[mesh_dim*j + k] *= inv_msf;
259 int restartIndex = restartIndex_;
261 std::pair<int,double> lastTimeStep = meshData->get_input_ioss_region()->get_max_time();
262 restartIndex = 1+restartIndex+lastTimeStep.first;
266 meshData->read_defined_input_fields(restartIndex);
268 mesh.buildSubcells();
269 mesh.buildLocalElementIDs();
271 mesh.beginModification();
272 if (createEdgeBlocks_) {
273 mesh.buildLocalEdgeIDs();
276 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
277 mesh.buildLocalFaceIDs();
280 mesh.endModification();
282 if (userMeshScaling_) {
283 stk::mesh::Field<double>* coord_field = metaData.get_field<
double>(stk::topology::NODE_RANK,
"coordinates");
284 std::vector< const stk::mesh::FieldBase *> fields;
285 fields.push_back(coord_field);
287 stk::mesh::communicate_field_data(bulkData, fields);
291 mesh.setInitialStateTime(meshData->get_input_ioss_region()->get_state_time(restartIndex));
293 mesh.setInitialStateTime(0.0);
298 if(rebalancing_ ==
"default")
300 this->rebalance(mesh);
301 else if(rebalancing_ !=
"none")
304 "ERROR: Rebalancing was not set to a valid choice");
313 "Error, the parameter {name=\"File Name\","
315 "\nis required in parameter (sub)list \""<< paramList->
name() <<
"\"."
323 paramList->
set<
int>(
"Restart Index", -1);
326 paramList->
set(
"File Type",
"Exodus");
328 if(!paramList->
isSublist(
"Periodic BCs"))
329 paramList->
sublist(
"Periodic BCs");
333 p_bcs.
set<
int>(
"Count", 0);
335 if(!paramList->
isParameter(
"Levels of Uniform Refinement"))
336 paramList->
set<
int>(
"Levels of Uniform Refinement", 0);
339 paramList->
set<
bool>(
"Keep Percept Data",
false);
341 if(!paramList->
isParameter(
"Keep Percept Parent Elements"))
342 paramList->
set<
bool>(
"Keep Percept Parent Elements",
false);
345 paramList->
set<std::string>(
"Rebalancing",
"default");
349 paramList->
set<
bool>(
"Create Edge Blocks",
false);
353 paramList->
set<
bool>(
"Create Face Blocks",
false);
356 paramList->
set(
"Geometry File Name",
"");
360 setMyParamList(paramList);
362 fileName_ = paramList->
get<std::string>(
"File Name");
364 geometryName_ = paramList->
get<std::string>(
"Geometry File Name");
366 restartIndex_ = paramList->
get<
int>(
"Restart Index");
368 fileType_ = paramList->
get<std::string>(
"File Type");
373 meshScaleFactor_ = paramList->
get<
double>(
"Scale Factor");
374 userMeshScaling_ =
true;
378 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->
sublist(
"Periodic BCs")),periodicBCVec_,useBBoxSearch_);
380 keepPerceptData_ = paramList->
get<
bool>(
"Keep Percept Data");
382 keepPerceptParentElements_ = paramList->
get<
bool>(
"Keep Percept Parent Elements");
384 rebalancing_ = paramList->
get<std::string>(
"Rebalancing");
386 levelsOfRefinement_ = paramList->
get<
int>(
"Levels of Uniform Refinement");
388 createEdgeBlocks_ = paramList->
get<
bool>(
"Create Edge Blocks");
389 createFaceBlocks_ = paramList->
get<
bool>(
"Create Face Blocks");
397 if(validParams==Teuchos::null) {
399 validParams->
set<std::string>(
"File Name",
"<file name not set>",
"Name of exodus file to be read",
401 validParams->
set<std::string>(
"Geometry File Name",
"<file name not set>",
"Name of geometry file for refinement",
404 validParams->
set<
int>(
"Restart Index",-1,
"Index of solution to read in",
407 validParams->
set<std::string>(
"File Type",
"Exodus",
408 "Choose input file type - either \"Exodus\", \"Exodus Refinement\" or \"Pamgen\"",
410 #ifdef PANZER_HAVE_UMR
415 validParams->
set<
double>(
"Scale Factor", 1.0,
"Scale factor to apply to mesh after read",
419 bcs.
set<
int>(
"Count",0);
421 validParams->
set(
"Levels of Uniform Refinement",0,
"Number of levels of inline uniform mesh refinement");
423 validParams->
set(
"Keep Percept Data",
false,
"Keep the Percept mesh after uniform refinement is applied");
425 validParams->
set(
"Keep Percept Parent Elements",
false,
"Keep the parent element information in the Percept data");
427 validParams->
set(
"Rebalancing",
"default",
"The type of rebalancing to be performed on the mesh after creation (default, none)");
430 validParams->
set(
"Create Edge Blocks",
false,
"Create or copy edge blocks in the mesh");
431 validParams->
set(
"Create Face Blocks",
false,
"Create or copy face blocks in the mesh");
437 void STK_ExodusReaderFactory::registerElementBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
441 RCP<stk::mesh::MetaData> femMetaData = mesh.getMetaData();
446 const Ioss::ElementBlockContainer & elem_blocks = meshData.get_input_ioss_region()->get_element_blocks();
447 for(Ioss::ElementBlockContainer::const_iterator itr=elem_blocks.begin();itr!=elem_blocks.end();++itr) {
448 Ioss::GroupingEntity * entity = *itr;
449 const std::string & name = entity->name();
451 const stk::mesh::Part * part = femMetaData->get_part(name);
452 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(femMetaData->get_topology(*part));
453 const CellTopologyData * ct = cellTopo.getCellTopologyData();
456 mesh.addElementBlock(part->name(),ct);
458 if (createEdgeBlocks_) {
459 createUniqueEdgeTopologyMap(mesh, part);
461 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
462 createUniqueFaceTopologyMap(mesh, part);
467 template <
typename SetType>
468 void buildSetNames(
const SetType & setData,std::vector<std::string> & names)
471 for(
typename SetType::const_iterator itr=setData.begin();itr!=setData.end();++itr) {
472 Ioss::GroupingEntity * entity = *itr;
473 names.push_back(entity->name());
477 void STK_ExodusReaderFactory::registerSidesets(STK_Interface & mesh)
const
481 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
482 const stk::mesh::PartVector & parts = metaData->get_parts();
484 stk::mesh::PartVector::const_iterator partItr;
485 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
486 const stk::mesh::Part * part = *partItr;
487 const stk::mesh::PartVector & subsets = part->subsets();
488 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
489 const CellTopologyData * ct = cellTopo.getCellTopologyData();
493 if(part->primary_entity_rank()==mesh.getSideRank() && ct==0 && subsets.size()>0) {
495 "STK_ExodusReaderFactory::registerSidesets error - part \"" << part->name() <<
496 "\" has more than one subset");
499 const stk::mesh::Part * ss_part = subsets[0];
500 shards::CellTopology ss_cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*ss_part));
501 const CellTopologyData * ss_ct = ss_cellTopo.getCellTopologyData();
505 mesh.addSideset(part->name(),ss_ct);
510 void STK_ExodusReaderFactory::registerNodesets(STK_Interface & mesh)
const
514 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
515 const stk::mesh::PartVector & parts = metaData->get_parts();
517 stk::mesh::PartVector::const_iterator partItr;
518 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
519 const stk::mesh::Part * part = *partItr;
520 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
521 const CellTopologyData * ct = cellTopo.getCellTopologyData();
525 if(part->primary_entity_rank()==mesh.getNodeRank() && ct==0) {
529 mesh.addNodeset(part->name());
534 void STK_ExodusReaderFactory::registerEdgeBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
538 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
546 const Ioss::EdgeBlockContainer & edge_blocks = meshData.get_input_ioss_region()->get_edge_blocks();
547 for(Ioss::EdgeBlockContainer::const_iterator ebc_iter=edge_blocks.begin();ebc_iter!=edge_blocks.end();++ebc_iter) {
548 Ioss::GroupingEntity * entity = *ebc_iter;
549 const stk::mesh::Part * edgeBlockPart = metaData->get_part(entity->name());
550 const stk::topology edgeBlockTopo = metaData->get_topology(*edgeBlockPart);
552 for (
auto ebuet_iter : elemBlockUniqueEdgeTopologies_) {
553 std::string elemBlockName = ebuet_iter.first;
554 std::vector<stk::topology> uniqueEdgeTopologies = ebuet_iter.second;
556 auto find_result = std::find(uniqueEdgeTopologies.begin(), uniqueEdgeTopologies.end(), edgeBlockTopo);
557 if (find_result != uniqueEdgeTopologies.end()) {
558 mesh.addEdgeBlock(elemBlockName, edgeBlockPart->name(), edgeBlockTopo);
564 void STK_ExodusReaderFactory::registerFaceBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData)
const
568 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
576 const Ioss::FaceBlockContainer & face_blocks = meshData.get_input_ioss_region()->get_face_blocks();
577 for(Ioss::FaceBlockContainer::const_iterator fbc_itr=face_blocks.begin();fbc_itr!=face_blocks.end();++fbc_itr) {
578 Ioss::GroupingEntity * entity = *fbc_itr;
579 const stk::mesh::Part * faceBlockPart = metaData->get_part(entity->name());
580 const stk::topology faceBlockTopo = metaData->get_topology(*faceBlockPart);
582 for (
auto ebuft_iter : elemBlockUniqueFaceTopologies_) {
583 std::string elemBlockName = ebuft_iter.first;
584 std::vector<stk::topology> uniqueFaceTopologies = ebuft_iter.second;
586 auto find_result = std::find(uniqueFaceTopologies.begin(), uniqueFaceTopologies.end(), faceBlockTopo);
587 if (find_result != uniqueFaceTopologies.end()) {
588 mesh.addFaceBlock(elemBlockName, faceBlockPart->name(), faceBlockTopo);
594 bool topo_less (stk::topology &i,stk::topology &j) {
return (i.value() < j.value()); }
595 bool topo_equal (stk::topology &i,stk::topology &j) {
return (i.value() == j.value()); }
597 void STK_ExodusReaderFactory::createUniqueEdgeTopologyMap(STK_Interface & mesh,
const stk::mesh::Part *elemBlockPart)
const
604 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
605 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
607 std::vector<stk::topology> edge_topologies;
608 for (
unsigned i=0;i<elemBlockTopo.num_edges();i++) {
609 edge_topologies.push_back(elemBlockTopo.edge_topology(i));
611 std::sort(edge_topologies.begin(), edge_topologies.end(), topo_less);
612 std::vector<stk::topology>::iterator new_end;
613 new_end = std::unique(edge_topologies.begin(), edge_topologies.end(), topo_equal);
614 edge_topologies.resize( std::distance(edge_topologies.begin(),new_end) );
616 elemBlockUniqueEdgeTopologies_[elemBlockPart->name()] = edge_topologies;
619 void STK_ExodusReaderFactory::createUniqueFaceTopologyMap(STK_Interface & mesh,
const stk::mesh::Part *elemBlockPart)
const
626 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
627 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
629 std::vector<stk::topology> face_topologies;
630 for (
unsigned i=0;i<elemBlockTopo.num_faces();i++) {
631 face_topologies.push_back(elemBlockTopo.face_topology(i));
633 std::sort(face_topologies.begin(), face_topologies.end(), topo_less);
634 std::vector<stk::topology>::iterator new_end;
635 new_end = std::unique(face_topologies.begin(), face_topologies.end(), topo_equal);
636 face_topologies.resize( std::distance(face_topologies.begin(),new_end) );
638 elemBlockUniqueFaceTopologies_[elemBlockPart->name()] = face_topologies;
643 void STK_ExodusReaderFactory::addEdgeBlocks(STK_Interface & mesh)
const
652 for (
auto iter : elemBlockUniqueEdgeTopologies_) {
653 std::string elemBlockName =
iter.first;
654 std::vector<stk::topology> uniqueEdgeTopologies =
iter.second;
656 for (
auto topo : uniqueEdgeTopologies ) {
657 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
658 const stk::mesh::Part & edgeTopoPart = metaData->get_topology_root_part(topo);
660 stk::mesh::Selector owned_block;
661 owned_block = *elemBlockPart;
662 owned_block &= edgeTopoPart;
663 owned_block &= metaData->locally_owned_part();
666 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edge_block_name);
668 std::vector<stk::mesh::Entity> all_edges_for_topo;
669 bulkData->get_entities(mesh.getEdgeRank(),owned_block,all_edges_for_topo);
670 mesh.addEntitiesToEdgeBlock(all_edges_for_topo, edge_block);
677 void STK_ExodusReaderFactory::addFaceBlocks(STK_Interface & mesh)
const
686 for (
auto iter : elemBlockUniqueFaceTopologies_) {
687 std::string elemBlockName =
iter.first;
688 std::vector<stk::topology> uniqueFaceTopologies =
iter.second;
690 for (
auto topo : uniqueFaceTopologies ) {
691 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
692 const stk::mesh::Part & faceTopoPart = metaData->get_topology_root_part(topo);
694 stk::mesh::Selector owned_block;
695 owned_block = *elemBlockPart;
696 owned_block &= faceTopoPart;
697 owned_block &= metaData->locally_owned_part();
700 stk::mesh::Part * face_block = mesh.getFaceBlock(face_block_name);
702 std::vector<stk::mesh::Entity> all_faces_for_topo;
703 bulkData->get_entities(mesh.getFaceRank(),owned_block,all_faces_for_topo);
704 mesh.addEntitiesToFaceBlock(all_faces_for_topo, face_block);
709 void STK_ExodusReaderFactory::buildMetaData(stk::ParallelMachine , STK_Interface & mesh)
const
711 if (createEdgeBlocks_) {
715 for (
auto iter : elemBlockUniqueEdgeTopologies_) {
716 std::string elemBlockName =
iter.first;
717 std::vector<stk::topology> uniqueEdgeTopologies =
iter.second;
719 for (
auto topo : uniqueEdgeTopologies ) {
721 mesh.addEdgeBlock(elemBlockName, edge_block_name, topo);
725 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
729 for (
auto iter : elemBlockUniqueFaceTopologies_) {
730 std::string elemBlockName =
iter.first;
731 std::vector<stk::topology> uniqueFaceTopologies =
iter.second;
733 for (
auto topo : uniqueFaceTopologies ) {
735 mesh.addFaceBlock(elemBlockName, face_block_name, topo);
741 bool STK_ExodusReaderFactory::doPerceptRefinement()
const
743 return (fileType_!=
"Exodus Refinement") && (levelsOfRefinement_ > 0);
746 std::string STK_ExodusReaderFactory::mkBlockName(std::string base, std::string topo_name)
const
749 name = topo_name+
"_"+base;
750 std::transform(name.begin(), name.end(), name.begin(),
752 {
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)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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)