Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_ExodusReaderFactory.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #include <Teuchos_TimeMonitor.hpp>
12 #include <Teuchos_RCP.hpp>
13 #include <Teuchos_RCPStdSharedPtrConversions.hpp>
14 #include <PanzerAdaptersSTK_config.hpp>
15 
17 #include "Panzer_STK_Interface.hpp"
18 
19 #ifdef PANZER_HAVE_IOSS
20 
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>
29 
30 #ifdef PANZER_HAVE_UMR
31 #include <Ioumr_DatabaseIO.hpp>
32 #endif
33 
34 #include "Teuchos_StandardParameterEntryValidators.hpp"
35 
36 namespace panzer_stk {
37 
38 int getMeshDimension(const std::string & meshStr,
39  stk::ParallelMachine parallelMach,
40  const std::string & typeStr)
41 {
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());
48 }
49 
50 std::string fileTypeToIOSSType(const std::string & fileType)
51 {
52  std::string IOSSType;
53  if (fileType=="Exodus")
54  IOSSType = "exodusii";
55 #ifdef PANZER_HAVE_UMR
56  else if (fileType=="Exodus Refinement")
57  IOSSType = "Refinement";
58 #endif
59  else if (fileType=="Pamgen")
60  IOSSType = "pamgen";
61 
62  return IOSSType;
63 }
64 
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_("")
70 { }
71 
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_("")
77 { }
78 
79 Teuchos::RCP<STK_Interface> STK_ExodusReaderFactory::buildMesh(stk::ParallelMachine parallelMach) const
80 {
81  PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::buildMesh()");
82 
83  using Teuchos::RCP;
84  using Teuchos::rcp;
85 
86  RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
87 
88  // in here you would add your fields...but it is better to use
89  // the two step construction
90 
91  // this calls commit on meta data
92  mesh->initialize(parallelMach,false,doPerceptRefinement());
93 
94  completeMeshConstruction(*mesh,parallelMach);
95 
96  return mesh;
97 }
98 
103 Teuchos::RCP<STK_Interface> STK_ExodusReaderFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
104 {
105  PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::buildUncomittedMesh()");
106 
107  using Teuchos::RCP;
108  using Teuchos::rcp;
109 
110  // read in meta data
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));
114 
115  // add in "FAMILY_TREE" entity for doing refinement
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);
119 
120 #ifdef PANZER_HAVE_UMR
121  // this line registers Ioumr with Ioss
122  Ioumr::IOFactory::factory();
123 
124  meshData->property_add(Ioss::Property("GEOMETRY_FILE", geometryName_));
125  meshData->property_add(Ioss::Property("NUMBER_REFINEMENTS", levelsOfRefinement_));
126 #endif
127 
128  meshData->add_mesh_database(fileName_, fileTypeToIOSSType(fileType_), stk::io::READ_MESH);
129 
130  meshData->create_input_mesh();
131  RCP<stk::mesh::MetaData> metaData = Teuchos::rcp(meshData->meta_data_ptr());
132 
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()));
137 
138  // read in other transient fields, these will be useful later when
139  // trying to read other fields for use in solve
140  meshData->add_all_mesh_fields_as_input_fields();
141 
142  // store mesh data pointer for later use in initializing
143  // bulk data
144  mesh->getMetaData()->declare_attribute_with_delete(meshData);
145 
146  // build element blocks
147  registerElementBlocks(*mesh,*meshData);
148  registerSidesets(*mesh);
149  registerNodesets(*mesh);
150 
151  if (createEdgeBlocks_) {
152  registerEdgeBlocks(*mesh,*meshData);
153  }
154  if (createFaceBlocks_ && mesh->getMetaData()->spatial_dimension() > 2) {
155  registerFaceBlocks(*mesh,*meshData);
156  }
157 
158  buildMetaData(parallelMach, *mesh);
159 
160  mesh->addPeriodicBCs(periodicBCVec_);
161  mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
162 
163  return mesh;
164 }
165 
166 void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
167 {
168  PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::completeMeshConstruction()");
169 
170  using Teuchos::RCP;
171  using Teuchos::rcp;
172 
173 
174  if(not mesh.isInitialized()) {
175  mesh.initialize(parallelMach,true,doPerceptRefinement());
176  }
177 
178  // grab mesh data pointer to build the bulk data
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>());
183  // if const_cast is wrong ... why does it feel so right?
184  // I believe this is safe since we are basically hiding this object under the covers
185  // until the mesh construction can be completed...below I cleanup the object myself.
186  TEUCHOS_ASSERT(metaData.remove_attribute(meshData));
187  // remove the MeshData attribute
188 
189  // build mesh bulk data
190  meshData->populate_bulk_data();
191 
192  if (doPerceptRefinement()) {
193  const bool deleteParentElements = !keepPerceptParentElements_;
194  mesh.refineMesh(levelsOfRefinement_,deleteParentElements);
195  }
196 
197  // The following section of code is applicable if mesh scaling is
198  // turned on from the input file.
199  if (userMeshScaling_)
200  {
201  stk::mesh::Field<double>* coord_field = metaData.get_field<double>(stk::topology::NODE_RANK, "coordinates");
202 
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);
205 
206  int mesh_dim = mesh.getDimension();
207 
208  // Scale the mesh
209  const double inv_msf = 1.0/meshScaleFactor_;
210  for (size_t i=0; i < my_node_buckets.size(); ++i)
211  {
212  stk::mesh::Bucket& b = *(my_node_buckets[i]);
213  double* coordinate_data = field_data( *coord_field, b );
214 
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;
218  }
219  }
220  }
221  }
222 
223  // put in a negative index and (like python) the restart will be from the back
224  // (-1 is the last time step)
225  int restartIndex = restartIndex_;
226  if(restartIndex<0) {
227  std::pair<int,double> lastTimeStep = meshData->get_input_ioss_region()->get_max_time();
228  restartIndex = 1+restartIndex+lastTimeStep.first;
229  }
230 
231  // populate mesh fields with specific index
232  meshData->read_defined_input_fields(restartIndex);
233 
234  mesh.buildSubcells();
235  mesh.buildLocalElementIDs();
236 
237  mesh.beginModification();
238  if (createEdgeBlocks_) {
239  mesh.buildLocalEdgeIDs();
240  addEdgeBlocks(mesh);
241  }
242  if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
243  mesh.buildLocalFaceIDs();
244  addFaceBlocks(mesh);
245  }
246 
247  {
248  const bool find_and_set_shared_nodes_in_stk = false;
249  mesh.endModification(find_and_set_shared_nodes_in_stk);
250  }
251 
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);
256 
257  stk::mesh::communicate_field_data(bulkData, fields);
258  }
259 
260  if(restartIndex>0) // process_input_request is a no-op if restartIndex<=0 ... thus there would be no inital time
261  mesh.setInitialStateTime(meshData->get_input_ioss_region()->get_state_time(restartIndex));
262  else
263  mesh.setInitialStateTime(0.0); // no initial time to speak, might as well use 0.0
264 
265  // clean up mesh data object
266  delete meshData;
267 
268  if(rebalancing_ == "default")
269  // calls Stk_MeshFactory::rebalance
270  this->rebalance(mesh);
271  else if(rebalancing_ != "none")
272  {
273  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
274  "ERROR: Rebalancing was not set to a valid choice");
275  }
276 }
277 
279 void STK_ExodusReaderFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
280 {
281  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(!paramList->isParameter("File Name"),
283  "Error, the parameter {name=\"File Name\","
284  "type=\"string\""
285  "\nis required in parameter (sub)list \""<< paramList->name() <<"\"."
286  "\n\nThe parsed parameter parameter list is: \n" << paramList->currentParametersString()
287  );
288 
289  // Set default values here. Not all the params should be set so this
290  // has to be done manually as opposed to using
291  // validateParametersAndSetDefaults().
292  if(!paramList->isParameter("Restart Index"))
293  paramList->set<int>("Restart Index", -1);
294 
295  if(!paramList->isParameter("File Type"))
296  paramList->set("File Type", "Exodus");
297 
298  if(!paramList->isSublist("Periodic BCs"))
299  paramList->sublist("Periodic BCs");
300 
301  Teuchos::ParameterList& p_bcs = paramList->sublist("Periodic BCs");
302  if (!p_bcs.isParameter("Count"))
303  p_bcs.set<int>("Count", 0);
304 
305  if(!paramList->isParameter("Levels of Uniform Refinement"))
306  paramList->set<int>("Levels of Uniform Refinement", 0);
307 
308  if(!paramList->isParameter("Keep Percept Data"))
309  paramList->set<bool>("Keep Percept Data", false);
310 
311  if(!paramList->isParameter("Keep Percept Parent Elements"))
312  paramList->set<bool>("Keep Percept Parent Elements", false);
313 
314  if(!paramList->isParameter("Rebalancing"))
315  paramList->set<std::string>("Rebalancing", "default");
316 
317  if(!paramList->isParameter("Create Edge Blocks"))
318  // default to false to prevent massive exodiff test failures
319  paramList->set<bool>("Create Edge Blocks", false);
320 
321  if(!paramList->isParameter("Create Face Blocks"))
322  // default to false to prevent massive exodiff test failures
323  paramList->set<bool>("Create Face Blocks", false);
324 
325  if(!paramList->isParameter("Geometry File Name"))
326  paramList->set("Geometry File Name", "");
327 
328  paramList->validateParameters(*getValidParameters(),0);
329 
330  setMyParamList(paramList);
331 
332  fileName_ = paramList->get<std::string>("File Name");
333 
334  geometryName_ = paramList->get<std::string>("Geometry File Name");
335 
336  restartIndex_ = paramList->get<int>("Restart Index");
337 
338  fileType_ = paramList->get<std::string>("File Type");
339 
340  // get any mesh scale factor
341  if (paramList->isParameter("Scale Factor"))
342  {
343  meshScaleFactor_ = paramList->get<double>("Scale Factor");
344  userMeshScaling_ = true;
345  }
346 
347  // read in periodic boundary conditions
348  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
349 
350  keepPerceptData_ = paramList->get<bool>("Keep Percept Data");
351 
352  keepPerceptParentElements_ = paramList->get<bool>("Keep Percept Parent Elements");
353 
354  rebalancing_ = paramList->get<std::string>("Rebalancing");
355 
356  levelsOfRefinement_ = paramList->get<int>("Levels of Uniform Refinement");
357 
358  createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
359  createFaceBlocks_ = paramList->get<bool>("Create Face Blocks");
360 }
361 
363 Teuchos::RCP<const Teuchos::ParameterList> STK_ExodusReaderFactory::getValidParameters() const
364 {
365  static Teuchos::RCP<Teuchos::ParameterList> validParams;
366 
367  if(validParams==Teuchos::null) {
368  validParams = Teuchos::rcp(new Teuchos::ParameterList);
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",
373 
374  validParams->set<int>("Restart Index",-1,"Index of solution to read in",
375  Teuchos::rcp(new Teuchos::AnyNumberParameterEntryValidator(Teuchos::AnyNumberParameterEntryValidator::PREFER_INT,Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes(true))));
376 
377  validParams->set<std::string>("File Type", "Exodus",
378  "Choose input file type - either \"Exodus\", \"Exodus Refinement\" or \"Pamgen\"",
379  rcp(new Teuchos::StringValidator(Teuchos::tuple<std::string>("Exodus", "Pamgen"
380 #ifdef PANZER_HAVE_UMR
381  ,"Exodus Refinement"
382 #endif
383  ))));
384 
385  validParams->set<double>("Scale Factor", 1.0, "Scale factor to apply to mesh after read",
386  Teuchos::rcp(new Teuchos::AnyNumberParameterEntryValidator(Teuchos::AnyNumberParameterEntryValidator::PREFER_DOUBLE,Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes(true))));
387 
388  Teuchos::ParameterList & bcs = validParams->sublist("Periodic BCs");
389  bcs.set<int>("Count",0); // no default periodic boundary conditions
390 
391  validParams->set("Levels of Uniform Refinement",0,"Number of levels of inline uniform mesh refinement");
392 
393  validParams->set("Keep Percept Data",false,"Keep the Percept mesh after uniform refinement is applied");
394 
395  validParams->set("Keep Percept Parent Elements",false,"Keep the parent element information in the Percept data");
396 
397  validParams->set("Rebalancing","default","The type of rebalancing to be performed on the mesh after creation (default, none)");
398 
399  // default to false for backward compatibility
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");
402  }
403 
404  return validParams.getConst();
405 }
406 
407 void STK_ExodusReaderFactory::registerElementBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
408 {
409  using Teuchos::RCP;
410 
411  RCP<stk::mesh::MetaData> femMetaData = mesh.getMetaData();
412 
413  // here we use the Ioss interface because they don't add
414  // "bonus" element blocks and its easier to determine
415  // "real" element blocks versus STK-only blocks
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();
420 
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();
424 
425  TEUCHOS_ASSERT(ct!=0);
426  mesh.addElementBlock(part->name(),ct);
427 
428  if (createEdgeBlocks_) {
429  createUniqueEdgeTopologyMap(mesh, part);
430  }
431  if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
432  createUniqueFaceTopologyMap(mesh, part);
433  }
434  }
435 }
436 
437 template <typename SetType>
438 void buildSetNames(const SetType & setData,std::vector<std::string> & names)
439 {
440  // pull out all names for this set
441  for(typename SetType::const_iterator itr=setData.begin();itr!=setData.end();++itr) {
442  Ioss::GroupingEntity * entity = *itr;
443  names.push_back(entity->name());
444  }
445 }
446 
447 void STK_ExodusReaderFactory::registerSidesets(STK_Interface & mesh) const
448 {
449  using Teuchos::RCP;
450 
451  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
452  const stk::mesh::PartVector & parts = metaData->get_parts();
453 
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();
460 
461  // if a side part ==> this is a sideset: now storage is recursive
462  // on part contains all sub parts with consistent topology
463  if(part->primary_entity_rank()==mesh.getSideRank() && ct==0 && subsets.size()>0) {
464  TEUCHOS_TEST_FOR_EXCEPTION(subsets.size()!=1,std::runtime_error,
465  "STK_ExodusReaderFactory::registerSidesets error - part \"" << part->name() <<
466  "\" has more than one subset");
467 
468  // grab cell topology and name of subset part
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();
472 
473  // only add subset parts that have no topology
474  if(ss_ct!=0)
475  mesh.addSideset(part->name(),ss_ct);
476  }
477  }
478 }
479 
480 void STK_ExodusReaderFactory::registerNodesets(STK_Interface & mesh) const
481 {
482  using Teuchos::RCP;
483 
484  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
485  const stk::mesh::PartVector & parts = metaData->get_parts();
486 
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();
492 
493  // if a side part ==> this is a sideset: now storage is recursive
494  // on part contains all sub parts with consistent topology
495  if(part->primary_entity_rank()==mesh.getNodeRank() && ct==0) {
496 
497  // only add subset parts that have no topology
498  if(part->name()!=STK_Interface::nodesString)
499  mesh.addNodeset(part->name());
500  }
501  }
502 }
503 
504 void STK_ExodusReaderFactory::registerEdgeBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
505 {
506  using Teuchos::RCP;
507 
508  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
509 
510  /* For each edge block in the exodus file, check it's topology
511  * against the list of edge topologies for each element block.
512  * If it matches, add the edge block for that element block.
513  * This will add the edge block as a subset of the element
514  * block in the STK mesh.
515  */
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);
521 
522  for (auto ebuet_iter : elemBlockUniqueEdgeTopologies_) {
523  std::string elemBlockName = ebuet_iter.first;
524  std::vector<stk::topology> uniqueEdgeTopologies = ebuet_iter.second;
525 
526  auto find_result = std::find(uniqueEdgeTopologies.begin(), uniqueEdgeTopologies.end(), edgeBlockTopo);
527  if (find_result != uniqueEdgeTopologies.end()) {
528  mesh.addEdgeBlock(elemBlockName, edgeBlockPart->name(), edgeBlockTopo);
529  }
530  }
531  }
532 }
533 
534 void STK_ExodusReaderFactory::registerFaceBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
535 {
536  using Teuchos::RCP;
537 
538  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
539 
540  /* For each face block in the exodus file, check it's topology
541  * against the list of face topologies for each element block.
542  * If it matches, add the face block for that element block.
543  * This will add the face block as a subset of the element
544  * block in the STK mesh.
545  */
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);
551 
552  for (auto ebuft_iter : elemBlockUniqueFaceTopologies_) {
553  std::string elemBlockName = ebuft_iter.first;
554  std::vector<stk::topology> uniqueFaceTopologies = ebuft_iter.second;
555 
556  auto find_result = std::find(uniqueFaceTopologies.begin(), uniqueFaceTopologies.end(), faceBlockTopo);
557  if (find_result != uniqueFaceTopologies.end()) {
558  mesh.addFaceBlock(elemBlockName, faceBlockPart->name(), faceBlockTopo);
559  }
560  }
561  }
562 }
563 
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()); }
566 
567 void STK_ExodusReaderFactory::createUniqueEdgeTopologyMap(STK_Interface & mesh, const stk::mesh::Part *elemBlockPart) const
568 {
569  using Teuchos::RCP;
570 
571  /* For a given element block, add it's edge topologies to a vector,
572  * sort it, dedupe it and save it to the "unique topo" map.
573  */
574  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
575  const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
576 
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));
580  }
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) );
585 
586  elemBlockUniqueEdgeTopologies_[elemBlockPart->name()] = edge_topologies;
587 }
588 
589 void STK_ExodusReaderFactory::createUniqueFaceTopologyMap(STK_Interface & mesh, const stk::mesh::Part *elemBlockPart) const
590 {
591  using Teuchos::RCP;
592 
593  /* For a given element block, add it's face topologies to a vector,
594  * sort it, dedupe it and save it to the "unique topo" map.
595  */
596  RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
597  const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
598 
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));
602  }
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) );
607 
608  elemBlockUniqueFaceTopologies_[elemBlockPart->name()] = face_topologies;
609 }
610 
611 // Pre-Condition: call beginModification() before entry
612 // Post-Condition: call endModification() after exit
613 void STK_ExodusReaderFactory::addEdgeBlocks(STK_Interface & mesh) const
614 {
615  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
616  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
617 
618  /* For each element block, iterate over it's edge topologies.
619  * For each edge topology, get the matching edge block and
620  * add all edges of that topology to the edge block. Make sure to include
621  * aura values - when they are marked shared, the part membership
622  * gets updated, if only ghosted the part memberships are not
623  * automatically updated.
624  */
625  for (auto iter : elemBlockUniqueEdgeTopologies_) {
626  std::string elemBlockName = iter.first;
627  std::vector<stk::topology> uniqueEdgeTopologies = iter.second;
628 
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);
632 
633  stk::mesh::Selector owned_block;
634  owned_block = *elemBlockPart;
635  owned_block &= edgeTopoPart;
636 
637  std::string edge_block_name = mkBlockName(panzer_stk::STK_Interface::edgeBlockString, topo.name());
638  stk::mesh::Part * edge_block = mesh.getEdgeBlock(edge_block_name);
639 
640  bulkData->change_entity_parts(owned_block, mesh.getEdgeRank(), stk::mesh::PartVector{edge_block});
641  }
642  }
643 }
644 
645 // Pre-Condition: call beginModification() before entry
646 // Post-Condition: call endModification() after exit
647 void STK_ExodusReaderFactory::addFaceBlocks(STK_Interface & mesh) const
648 {
649  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
650  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
651 
652  /* For each element block, iterate over it's face topologies. For
653  * each face topology, get the matching face block and add all
654  * faces of that topology to the face block. Make sure to include
655  * aura values - when they are marked shared, the part membership
656  * gets updated, if only ghosted the part memberships are not
657  * automatically updated.
658  */
659  for (auto iter : elemBlockUniqueFaceTopologies_) {
660  std::string elemBlockName = iter.first;
661  std::vector<stk::topology> uniqueFaceTopologies = iter.second;
662 
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);
666 
667  stk::mesh::Selector owned_block;
668  owned_block = *elemBlockPart;
669  owned_block &= faceTopoPart;
670 
671  std::string face_block_name = mkBlockName(panzer_stk::STK_Interface::faceBlockString, topo.name());
672  stk::mesh::Part * face_block = mesh.getFaceBlock(face_block_name);
673 
674  bulkData->change_entity_parts(owned_block, mesh.getFaceRank(), stk::mesh::PartVector{face_block});
675  }
676  }
677 }
678 
679 void STK_ExodusReaderFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
680 {
681  if (createEdgeBlocks_) {
682  /* For each element block, iterate over it's edge topologies.
683  * For each edge topology, create an edge block for that topology.
684  */
685  for (auto iter : elemBlockUniqueEdgeTopologies_) {
686  std::string elemBlockName = iter.first;
687  std::vector<stk::topology> uniqueEdgeTopologies = iter.second;
688 
689  for (auto topo : uniqueEdgeTopologies ) {
690  std::string edge_block_name = mkBlockName(panzer_stk::STK_Interface::edgeBlockString, topo.name());
691  mesh.addEdgeBlock(elemBlockName, edge_block_name, topo);
692  }
693  }
694  }
695  if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
696  /* For each element block, iterate over it's face topologies.
697  * For each face topology, create a face block for that topology.
698  */
699  for (auto iter : elemBlockUniqueFaceTopologies_) {
700  std::string elemBlockName = iter.first;
701  std::vector<stk::topology> uniqueFaceTopologies = iter.second;
702 
703  for (auto topo : uniqueFaceTopologies ) {
704  std::string face_block_name = mkBlockName(panzer_stk::STK_Interface::faceBlockString, topo.name());
705  mesh.addFaceBlock(elemBlockName, face_block_name, topo);
706  }
707  }
708  }
709 }
710 
711 bool STK_ExodusReaderFactory::doPerceptRefinement() const
712 {
713  return (fileType_!="Exodus Refinement") && (levelsOfRefinement_ > 0);
714 }
715 
716 std::string STK_ExodusReaderFactory::mkBlockName(std::string base, std::string topo_name) const
717 {
718  std::string name;
719  name = topo_name+"_"+base;
720  std::transform(name.begin(), name.end(), name.begin(),
721  [](const char c)
722  { return char(std::tolower(c)); });
723  return name;
724 }
725 
726 }
727 
728 #endif
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
TransListIter iter
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)