Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_ExodusReaderFactory.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43
44#include <Teuchos_TimeMonitor.hpp>
45#include <PanzerAdaptersSTK_config.hpp>
46
49
50#ifdef PANZER_HAVE_IOSS
51
52#include <Ionit_Initializer.h>
53#include <Ioss_ElementBlock.h>
54#include <Ioss_EdgeBlock.h>
55#include <Ioss_FaceBlock.h>
56#include <Ioss_Region.h>
57#include <stk_mesh/base/GetBuckets.hpp>
58#include <stk_io/StkMeshIoBroker.hpp>
59#include <stk_io/IossBridge.hpp>
60#include <stk_mesh/base/FieldParallel.hpp>
61
62#ifdef PANZER_HAVE_UMR
63#include <Ioumr_DatabaseIO.hpp>
64#endif
65
66#include "Teuchos_StandardParameterEntryValidators.hpp"
67
68namespace panzer_stk {
69
70int getMeshDimension(const std::string & meshStr,
71 stk::ParallelMachine parallelMach,
72 const std::string & typeStr)
73{
74 stk::io::StkMeshIoBroker meshData(parallelMach);
75 meshData.property_add(Ioss::Property("LOWER_CASE_VARIABLE_NAMES", false));
76 meshData.add_mesh_database(meshStr, fileTypeToIOSSType(typeStr), stk::io::READ_MESH);
77 meshData.create_input_mesh();
78 return Teuchos::as<int>(meshData.meta_data_ptr()->spatial_dimension());
79}
80
81std::string fileTypeToIOSSType(const std::string & fileType)
82{
83 std::string IOSSType;
84 if (fileType=="Exodus")
85 IOSSType = "exodusii";
86#ifdef PANZER_HAVE_UMR
87 else if (fileType=="Exodus Refinement")
88 IOSSType = "Refinement";
89#endif
90 else if (fileType=="Pamgen")
91 IOSSType = "pamgen";
92
93 return IOSSType;
94}
95
96STK_ExodusReaderFactory::STK_ExodusReaderFactory()
97 : fileName_(""), fileType_(""), restartIndex_(0), userMeshScaling_(false), keepPerceptData_(false),
98 keepPerceptParentElements_(false), rebalancing_("default"),
99meshScaleFactor_(0.0), levelsOfRefinement_(0),
100 createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_("")
101{ }
102
103STK_ExodusReaderFactory::STK_ExodusReaderFactory(const std::string & fileName,
104 const int restartIndex)
105 : fileName_(fileName), fileType_("Exodus"), restartIndex_(restartIndex), userMeshScaling_(false),
106 keepPerceptData_(false), keepPerceptParentElements_(false), rebalancing_("default"),
107 meshScaleFactor_(0.0), levelsOfRefinement_(0), createEdgeBlocks_(false), createFaceBlocks_(false), geometryName_("")
108{ }
109
110Teuchos::RCP<STK_Interface> STK_ExodusReaderFactory::buildMesh(stk::ParallelMachine parallelMach) const
111{
112 PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::buildMesh()");
113
114 using Teuchos::RCP;
115 using Teuchos::rcp;
116
117 RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
118
119 // in here you would add your fields...but it is better to use
120 // the two step construction
121
122 // this calls commit on meta data
123 mesh->initialize(parallelMach,false,doPerceptRefinement());
124
125 completeMeshConstruction(*mesh,parallelMach);
126
127 return mesh;
128}
129
134Teuchos::RCP<STK_Interface> STK_ExodusReaderFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
135{
136 PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::buildUncomittedMesh()");
137
138 using Teuchos::RCP;
139 using Teuchos::rcp;
140
141 // read in meta data
142 stk::io::StkMeshIoBroker* meshData = new stk::io::StkMeshIoBroker(parallelMach);
143 meshData->property_add(Ioss::Property("LOWER_CASE_VARIABLE_NAMES", false));
144
145 // add in "FAMILY_TREE" entity for doing refinement
146 std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
147 entity_rank_names.push_back("FAMILY_TREE");
148 meshData->set_rank_name_vector(entity_rank_names);
149
150#ifdef PANZER_HAVE_UMR
151 // this line registers Ioumr with Ioss
152 Ioumr::IOFactory::factory();
153
154 meshData->property_add(Ioss::Property("GEOMETRY_FILE", geometryName_));
155 meshData->property_add(Ioss::Property("NUMBER_REFINEMENTS", levelsOfRefinement_));
156#endif
157
158 meshData->add_mesh_database(fileName_, fileTypeToIOSSType(fileType_), stk::io::READ_MESH);
159
160 meshData->create_input_mesh();
161 RCP<stk::mesh::MetaData> metaData = Teuchos::rcp(meshData->meta_data_ptr());
162
163 RCP<STK_Interface> mesh = rcp(new STK_Interface(metaData));
164 mesh->initializeFromMetaData();
165 mesh->instantiateBulkData(parallelMach);
166 meshData->set_bulk_data(Teuchos::get_shared_ptr(mesh->getBulkData()));
167
168 // read in other transient fields, these will be useful later when
169 // trying to read other fields for use in solve
170 meshData->add_all_mesh_fields_as_input_fields();
171
172 // store mesh data pointer for later use in initializing
173 // bulk data
174 mesh->getMetaData()->declare_attribute_with_delete(meshData);
175
176 // build element blocks
177 registerElementBlocks(*mesh,*meshData);
178 registerSidesets(*mesh);
179 registerNodesets(*mesh);
180
181 if (createEdgeBlocks_) {
182 registerEdgeBlocks(*mesh,*meshData);
183 }
184 if (createFaceBlocks_ && mesh->getMetaData()->spatial_dimension() > 2) {
185 registerFaceBlocks(*mesh,*meshData);
186 }
187
188 buildMetaData(parallelMach, *mesh);
189
190 mesh->addPeriodicBCs(periodicBCVec_);
191 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
192
193 return mesh;
194}
195
196void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
197{
198 PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::completeMeshConstruction()");
199
200 using Teuchos::RCP;
201 using Teuchos::rcp;
202
203
204 if(not mesh.isInitialized()) {
205 mesh.initialize(parallelMach,true,doPerceptRefinement());
206 }
207
208 // grab mesh data pointer to build the bulk data
209 stk::mesh::MetaData & metaData = *mesh.getMetaData();
210 stk::mesh::BulkData & bulkData = *mesh.getBulkData();
211 stk::io::StkMeshIoBroker * meshData =
212 const_cast<stk::io::StkMeshIoBroker *>(metaData.get_attribute<stk::io::StkMeshIoBroker>());
213 // if const_cast is wrong ... why does it feel so right?
214 // I believe this is safe since we are basically hiding this object under the covers
215 // until the mesh construction can be completed...below I cleanup the object myself.
216 TEUCHOS_ASSERT(metaData.remove_attribute(meshData));
217 // remove the MeshData attribute
218
219 // build mesh bulk data
220 meshData->populate_bulk_data();
221
222 if (doPerceptRefinement()) {
223 const bool deleteParentElements = !keepPerceptParentElements_;
224 mesh.refineMesh(levelsOfRefinement_,deleteParentElements);
225 }
226
227 // The following section of code is applicable if mesh scaling is
228 // turned on from the input file.
229 if (userMeshScaling_)
230 {
231 stk::mesh::Field<double,stk::mesh::Cartesian>* coord_field =
232 metaData.get_field<stk::mesh::Field<double, stk::mesh::Cartesian> >(stk::topology::NODE_RANK, "coordinates");
233
234 stk::mesh::Selector select_all_local = metaData.locally_owned_part() | metaData.globally_shared_part();
235 stk::mesh::BucketVector const& my_node_buckets = bulkData.get_buckets(stk::topology::NODE_RANK, select_all_local);
236
237 int mesh_dim = mesh.getDimension();
238
239 // Scale the mesh
240 const double inv_msf = 1.0/meshScaleFactor_;
241 for (size_t i=0; i < my_node_buckets.size(); ++i)
242 {
243 stk::mesh::Bucket& b = *(my_node_buckets[i]);
244 double* coordinate_data = field_data( *coord_field, b );
245
246 for (size_t j=0; j < b.size(); ++j) {
247 for (int k=0; k < mesh_dim; ++k) {
248 coordinate_data[mesh_dim*j + k] *= inv_msf;
249 }
250 }
251 }
252 }
253
254 // put in a negative index and (like python) the restart will be from the back
255 // (-1 is the last time step)
256 int restartIndex = restartIndex_;
257 if(restartIndex<0) {
258 std::pair<int,double> lastTimeStep = meshData->get_input_io_region()->get_max_time();
259 restartIndex = 1+restartIndex+lastTimeStep.first;
260 }
261
262 // populate mesh fields with specific index
263 meshData->read_defined_input_fields(restartIndex);
264
265 mesh.buildSubcells();
266 mesh.buildLocalElementIDs();
267
268 mesh.beginModification();
269 if (createEdgeBlocks_) {
270 mesh.buildLocalEdgeIDs();
271 addEdgeBlocks(mesh);
272 }
273 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
274 mesh.buildLocalFaceIDs();
275 addFaceBlocks(mesh);
276 }
277 mesh.endModification();
278
279 if (userMeshScaling_) {
280 stk::mesh::Field<double,stk::mesh::Cartesian>* coord_field =
281 metaData.get_field<stk::mesh::Field<double, stk::mesh::Cartesian> >(stk::topology::NODE_RANK, "coordinates");
282 std::vector< const stk::mesh::FieldBase *> fields;
283 fields.push_back(coord_field);
284
285 stk::mesh::communicate_field_data(bulkData, fields);
286 }
287
288 if(restartIndex>0) // process_input_request is a no-op if restartIndex<=0 ... thus there would be no inital time
289 mesh.setInitialStateTime(meshData->get_input_io_region()->get_state_time(restartIndex));
290 else
291 mesh.setInitialStateTime(0.0); // no initial time to speak, might as well use 0.0
292
293 // clean up mesh data object
294 delete meshData;
295
296 if(rebalancing_ == "default")
297 // calls Stk_MeshFactory::rebalance
298 this->rebalance(mesh);
299 else if(rebalancing_ != "none")
300 {
301 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
302 "ERROR: Rebalancing was not set to a valid choice");
303 }
304}
305
307void STK_ExodusReaderFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
308{
309 TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(!paramList->isParameter("File Name"),
310 Teuchos::Exceptions::InvalidParameterName,
311 "Error, the parameter {name=\"File Name\","
312 "type=\"string\""
313 "\nis required in parameter (sub)list \""<< paramList->name() <<"\"."
314 "\n\nThe parsed parameter parameter list is: \n" << paramList->currentParametersString()
315 );
316
317 // Set default values here. Not all the params should be set so this
318 // has to be done manually as opposed to using
319 // validateParametersAndSetDefaults().
320 if(!paramList->isParameter("Restart Index"))
321 paramList->set<int>("Restart Index", -1);
322
323 if(!paramList->isParameter("File Type"))
324 paramList->set("File Type", "Exodus");
325
326 if(!paramList->isSublist("Periodic BCs"))
327 paramList->sublist("Periodic BCs");
328
329 Teuchos::ParameterList& p_bcs = paramList->sublist("Periodic BCs");
330 if (!p_bcs.isParameter("Count"))
331 p_bcs.set<int>("Count", 0);
332
333 if(!paramList->isParameter("Levels of Uniform Refinement"))
334 paramList->set<int>("Levels of Uniform Refinement", 0);
335
336 if(!paramList->isParameter("Keep Percept Data"))
337 paramList->set<bool>("Keep Percept Data", false);
338
339 if(!paramList->isParameter("Keep Percept Parent Elements"))
340 paramList->set<bool>("Keep Percept Parent Elements", false);
341
342 if(!paramList->isParameter("Rebalancing"))
343 paramList->set<std::string>("Rebalancing", "default");
344
345 if(!paramList->isParameter("Create Edge Blocks"))
346 // default to false to prevent massive exodiff test failures
347 paramList->set<bool>("Create Edge Blocks", false);
348
349 if(!paramList->isParameter("Create Face Blocks"))
350 // default to false to prevent massive exodiff test failures
351 paramList->set<bool>("Create Face Blocks", false);
352
353 if(!paramList->isParameter("Geometry File Name"))
354 paramList->set("Geometry File Name", "");
355
356 paramList->validateParameters(*getValidParameters(),0);
357
358 setMyParamList(paramList);
359
360 fileName_ = paramList->get<std::string>("File Name");
361
362 geometryName_ = paramList->get<std::string>("Geometry File Name");
363
364 restartIndex_ = paramList->get<int>("Restart Index");
365
366 fileType_ = paramList->get<std::string>("File Type");
367
368 // get any mesh scale factor
369 if (paramList->isParameter("Scale Factor"))
370 {
371 meshScaleFactor_ = paramList->get<double>("Scale Factor");
372 userMeshScaling_ = true;
373 }
374
375 // read in periodic boundary conditions
376 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
377
378 keepPerceptData_ = paramList->get<bool>("Keep Percept Data");
379
380 keepPerceptParentElements_ = paramList->get<bool>("Keep Percept Parent Elements");
381
382 rebalancing_ = paramList->get<std::string>("Rebalancing");
383
384 levelsOfRefinement_ = paramList->get<int>("Levels of Uniform Refinement");
385
386 createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
387 createFaceBlocks_ = paramList->get<bool>("Create Face Blocks");
388}
389
391Teuchos::RCP<const Teuchos::ParameterList> STK_ExodusReaderFactory::getValidParameters() const
392{
393 static Teuchos::RCP<Teuchos::ParameterList> validParams;
394
395 if(validParams==Teuchos::null) {
396 validParams = Teuchos::rcp(new Teuchos::ParameterList);
397 validParams->set<std::string>("File Name","<file name not set>","Name of exodus file to be read",
398 Teuchos::rcp(new Teuchos::FileNameValidator));
399 validParams->set<std::string>("Geometry File Name","<file name not set>","Name of geometry file for refinement",
400 Teuchos::rcp(new Teuchos::FileNameValidator));
401
402 validParams->set<int>("Restart Index",-1,"Index of solution to read in",
403 Teuchos::rcp(new Teuchos::AnyNumberParameterEntryValidator(Teuchos::AnyNumberParameterEntryValidator::PREFER_INT,Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes(true))));
404
405 Teuchos::setStringToIntegralParameter<int>("File Type",
406 "Exodus",
407 "Choose input file type - either \"Exodus\", \"Exodus Refinement\" or \"Pamgen\"",
408 Teuchos::tuple<std::string>("Exodus","Pamgen"
409#ifdef PANZER_HAVE_UMR
410 ,"Exodus Refinement"
411#endif
412 ),
413 validParams.get()
414 );
415
416 validParams->set<double>("Scale Factor", 1.0, "Scale factor to apply to mesh after read",
417 Teuchos::rcp(new Teuchos::AnyNumberParameterEntryValidator(Teuchos::AnyNumberParameterEntryValidator::PREFER_DOUBLE,Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes(true))));
418
419 Teuchos::ParameterList & bcs = validParams->sublist("Periodic BCs");
420 bcs.set<int>("Count",0); // no default periodic boundary conditions
421
422 validParams->set("Levels of Uniform Refinement",0,"Number of levels of inline uniform mesh refinement");
423
424 validParams->set("Keep Percept Data",false,"Keep the Percept mesh after uniform refinement is applied");
425
426 validParams->set("Keep Percept Parent Elements",false,"Keep the parent element information in the Percept data");
427
428 validParams->set("Rebalancing","default","The type of rebalancing to be performed on the mesh after creation (default, none)");
429
430 // default to false for backward compatibility
431 validParams->set("Create Edge Blocks",false,"Create or copy edge blocks in the mesh");
432 validParams->set("Create Face Blocks",false,"Create or copy face blocks in the mesh");
433 }
434
435 return validParams.getConst();
436}
437
438void STK_ExodusReaderFactory::registerElementBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
439{
440 using Teuchos::RCP;
441
442 RCP<stk::mesh::MetaData> femMetaData = mesh.getMetaData();
443
444 // here we use the Ioss interface because they don't add
445 // "bonus" element blocks and its easier to determine
446 // "real" element blocks versus STK-only blocks
447 const Ioss::ElementBlockContainer & elem_blocks = meshData.get_input_io_region()->get_element_blocks();
448 for(Ioss::ElementBlockContainer::const_iterator itr=elem_blocks.begin();itr!=elem_blocks.end();++itr) {
449 Ioss::GroupingEntity * entity = *itr;
450 const std::string & name = entity->name();
451
452 const stk::mesh::Part * part = femMetaData->get_part(name);
453 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(femMetaData->get_topology(*part));
454 const CellTopologyData * ct = cellTopo.getCellTopologyData();
455
456 TEUCHOS_ASSERT(ct!=0);
457 mesh.addElementBlock(part->name(),ct);
458
459 if (createEdgeBlocks_) {
460 createUniqueEdgeTopologyMap(mesh, part);
461 }
462 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
463 createUniqueFaceTopologyMap(mesh, part);
464 }
465 }
466}
467
468template <typename SetType>
469void buildSetNames(const SetType & setData,std::vector<std::string> & names)
470{
471 // pull out all names for this set
472 for(typename SetType::const_iterator itr=setData.begin();itr!=setData.end();++itr) {
473 Ioss::GroupingEntity * entity = *itr;
474 names.push_back(entity->name());
475 }
476}
477
478void STK_ExodusReaderFactory::registerSidesets(STK_Interface & mesh) const
479{
480 using Teuchos::RCP;
481
482 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
483 const stk::mesh::PartVector & parts = metaData->get_parts();
484
485 stk::mesh::PartVector::const_iterator partItr;
486 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
487 const stk::mesh::Part * part = *partItr;
488 const stk::mesh::PartVector & subsets = part->subsets();
489 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
490 const CellTopologyData * ct = cellTopo.getCellTopologyData();
491
492 // if a side part ==> this is a sideset: now storage is recursive
493 // on part contains all sub parts with consistent topology
494 if(part->primary_entity_rank()==mesh.getSideRank() && ct==0 && subsets.size()>0) {
495 TEUCHOS_TEST_FOR_EXCEPTION(subsets.size()!=1,std::runtime_error,
496 "STK_ExodusReaderFactory::registerSidesets error - part \"" << part->name() <<
497 "\" has more than one subset");
498
499 // grab cell topology and name of subset part
500 const stk::mesh::Part * ss_part = subsets[0];
501 shards::CellTopology ss_cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*ss_part));
502 const CellTopologyData * ss_ct = ss_cellTopo.getCellTopologyData();
503
504 // only add subset parts that have no topology
505 if(ss_ct!=0)
506 mesh.addSideset(part->name(),ss_ct);
507 }
508 }
509}
510
511void STK_ExodusReaderFactory::registerNodesets(STK_Interface & mesh) const
512{
513 using Teuchos::RCP;
514
515 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
516 const stk::mesh::PartVector & parts = metaData->get_parts();
517
518 stk::mesh::PartVector::const_iterator partItr;
519 for(partItr=parts.begin();partItr!=parts.end();++partItr) {
520 const stk::mesh::Part * part = *partItr;
521 shards::CellTopology cellTopo = stk::mesh::get_cell_topology(metaData->get_topology(*part));
522 const CellTopologyData * ct = cellTopo.getCellTopologyData();
523
524 // if a side part ==> this is a sideset: now storage is recursive
525 // on part contains all sub parts with consistent topology
526 if(part->primary_entity_rank()==mesh.getNodeRank() && ct==0) {
527
528 // only add subset parts that have no topology
529 if(part->name()!=STK_Interface::nodesString)
530 mesh.addNodeset(part->name());
531 }
532 }
533}
534
535void STK_ExodusReaderFactory::registerEdgeBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
536{
537 using Teuchos::RCP;
538
539 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
540
541 /* For each edge block in the exodus file, check it's topology
542 * against the list of edge topologies for each element block.
543 * If it matches, add the edge block for that element block.
544 * This will add the edge block as a subset of the element
545 * block in the STK mesh.
546 */
547 const Ioss::EdgeBlockContainer & edge_blocks = meshData.get_input_io_region()->get_edge_blocks();
548 for(Ioss::EdgeBlockContainer::const_iterator ebc_iter=edge_blocks.begin();ebc_iter!=edge_blocks.end();++ebc_iter) {
549 Ioss::GroupingEntity * entity = *ebc_iter;
550 const stk::mesh::Part * edgeBlockPart = metaData->get_part(entity->name());
551 const stk::topology edgeBlockTopo = metaData->get_topology(*edgeBlockPart);
552
553 for (auto ebuet_iter : elemBlockUniqueEdgeTopologies_) {
554 std::string elemBlockName = ebuet_iter.first;
555 std::vector<stk::topology> uniqueEdgeTopologies = ebuet_iter.second;
556
557 auto find_result = std::find(uniqueEdgeTopologies.begin(), uniqueEdgeTopologies.end(), edgeBlockTopo);
558 if (find_result != uniqueEdgeTopologies.end()) {
559 mesh.addEdgeBlock(elemBlockName, edgeBlockPart->name(), edgeBlockTopo);
560 }
561 }
562 }
563}
564
565void STK_ExodusReaderFactory::registerFaceBlocks(STK_Interface & mesh,stk::io::StkMeshIoBroker & meshData) const
566{
567 using Teuchos::RCP;
568
569 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
570
571 /* For each face block in the exodus file, check it's topology
572 * against the list of face topologies for each element block.
573 * If it matches, add the face block for that element block.
574 * This will add the face block as a subset of the element
575 * block in the STK mesh.
576 */
577 const Ioss::FaceBlockContainer & face_blocks = meshData.get_input_io_region()->get_face_blocks();
578 for(Ioss::FaceBlockContainer::const_iterator fbc_itr=face_blocks.begin();fbc_itr!=face_blocks.end();++fbc_itr) {
579 Ioss::GroupingEntity * entity = *fbc_itr;
580 const stk::mesh::Part * faceBlockPart = metaData->get_part(entity->name());
581 const stk::topology faceBlockTopo = metaData->get_topology(*faceBlockPart);
582
583 for (auto ebuft_iter : elemBlockUniqueFaceTopologies_) {
584 std::string elemBlockName = ebuft_iter.first;
585 std::vector<stk::topology> uniqueFaceTopologies = ebuft_iter.second;
586
587 auto find_result = std::find(uniqueFaceTopologies.begin(), uniqueFaceTopologies.end(), faceBlockTopo);
588 if (find_result != uniqueFaceTopologies.end()) {
589 mesh.addFaceBlock(elemBlockName, faceBlockPart->name(), faceBlockTopo);
590 }
591 }
592 }
593}
594
595bool topo_less (stk::topology &i,stk::topology &j) { return (i.value() < j.value()); }
596bool topo_equal (stk::topology &i,stk::topology &j) { return (i.value() == j.value()); }
597
598void STK_ExodusReaderFactory::createUniqueEdgeTopologyMap(STK_Interface & mesh, const stk::mesh::Part *elemBlockPart) const
599{
600 using Teuchos::RCP;
601
602 /* For a given element block, add it's edge topologies to a vector,
603 * sort it, dedupe it and save it to the "unique topo" map.
604 */
605 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
606 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
607
608 std::vector<stk::topology> edge_topologies;
609 for (unsigned i=0;i<elemBlockTopo.num_edges();i++) {
610 edge_topologies.push_back(elemBlockTopo.edge_topology(i));
611 }
612 std::sort(edge_topologies.begin(), edge_topologies.end(), topo_less);
613 std::vector<stk::topology>::iterator new_end;
614 new_end = std::unique(edge_topologies.begin(), edge_topologies.end(), topo_equal);
615 edge_topologies.resize( std::distance(edge_topologies.begin(),new_end) );
616
617 elemBlockUniqueEdgeTopologies_[elemBlockPart->name()] = edge_topologies;
618}
619
620void STK_ExodusReaderFactory::createUniqueFaceTopologyMap(STK_Interface & mesh, const stk::mesh::Part *elemBlockPart) const
621{
622 using Teuchos::RCP;
623
624 /* For a given element block, add it's face topologies to a vector,
625 * sort it, dedupe it and save it to the "unique topo" map.
626 */
627 RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
628 const stk::topology elemBlockTopo = metaData->get_topology(*elemBlockPart);
629
630 std::vector<stk::topology> face_topologies;
631 for (unsigned i=0;i<elemBlockTopo.num_faces();i++) {
632 face_topologies.push_back(elemBlockTopo.face_topology(i));
633 }
634 std::sort(face_topologies.begin(), face_topologies.end(), topo_less);
635 std::vector<stk::topology>::iterator new_end;
636 new_end = std::unique(face_topologies.begin(), face_topologies.end(), topo_equal);
637 face_topologies.resize( std::distance(face_topologies.begin(),new_end) );
638
639 elemBlockUniqueFaceTopologies_[elemBlockPart->name()] = face_topologies;
640}
641
642// Pre-Condition: call beginModification() before entry
643// Post-Condition: call endModification() after exit
644void STK_ExodusReaderFactory::addEdgeBlocks(STK_Interface & mesh) const
645{
646 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
647 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
648
649 /* For each element block, iterate over it's edge topologies.
650 * For each edge topology, get the matching edge block and
651 * add all edges of that topology to the edge block.
652 */
653 for (auto iter : elemBlockUniqueEdgeTopologies_) {
654 std::string elemBlockName = iter.first;
655 std::vector<stk::topology> uniqueEdgeTopologies = iter.second;
656
657 for (auto topo : uniqueEdgeTopologies ) {
658 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
659 const stk::mesh::Part & edgeTopoPart = metaData->get_topology_root_part(topo);
660
661 stk::mesh::Selector owned_block;
662 owned_block = *elemBlockPart;
663 owned_block &= edgeTopoPart;
664 owned_block &= metaData->locally_owned_part();
665
666 std::string edge_block_name = mkBlockName(panzer_stk::STK_Interface::edgeBlockString, topo.name());
667 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edge_block_name);
668
669 std::vector<stk::mesh::Entity> all_edges_for_topo;
670 bulkData->get_entities(mesh.getEdgeRank(),owned_block,all_edges_for_topo);
671 mesh.addEntitiesToEdgeBlock(all_edges_for_topo, edge_block);
672 }
673 }
674}
675
676// Pre-Condition: call beginModification() before entry
677// Post-Condition: call endModification() after exit
678void STK_ExodusReaderFactory::addFaceBlocks(STK_Interface & mesh) const
679{
680 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
681 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
682
683 /* For each element block, iterate over it's face topologies.
684 * For each face topology, get the matching face block and
685 * add all faces of that topology to the face block.
686 */
687 for (auto iter : elemBlockUniqueFaceTopologies_) {
688 std::string elemBlockName = iter.first;
689 std::vector<stk::topology> uniqueFaceTopologies = iter.second;
690
691 for (auto topo : uniqueFaceTopologies ) {
692 const stk::mesh::Part * elemBlockPart = metaData->get_part(elemBlockName);
693 const stk::mesh::Part & faceTopoPart = metaData->get_topology_root_part(topo);
694
695 stk::mesh::Selector owned_block;
696 owned_block = *elemBlockPart;
697 owned_block &= faceTopoPart;
698 owned_block &= metaData->locally_owned_part();
699
700 std::string face_block_name = mkBlockName(panzer_stk::STK_Interface::faceBlockString, topo.name());
701 stk::mesh::Part * face_block = mesh.getFaceBlock(face_block_name);
702
703 std::vector<stk::mesh::Entity> all_faces_for_topo;
704 bulkData->get_entities(mesh.getFaceRank(),owned_block,all_faces_for_topo);
705 mesh.addEntitiesToFaceBlock(all_faces_for_topo, face_block);
706 }
707 }
708}
709
710void STK_ExodusReaderFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
711{
712 if (createEdgeBlocks_) {
713 /* For each element block, iterate over it's edge topologies.
714 * For each edge topology, create an edge block for that topology.
715 */
716 for (auto iter : elemBlockUniqueEdgeTopologies_) {
717 std::string elemBlockName = iter.first;
718 std::vector<stk::topology> uniqueEdgeTopologies = iter.second;
719
720 for (auto topo : uniqueEdgeTopologies ) {
721 std::string edge_block_name = mkBlockName(panzer_stk::STK_Interface::edgeBlockString, topo.name());
722 mesh.addEdgeBlock(elemBlockName, edge_block_name, topo);
723 }
724 }
725 }
726 if (createFaceBlocks_ && mesh.getMetaData()->spatial_dimension() > 2) {
727 /* For each element block, iterate over it's face topologies.
728 * For each face topology, create a face block for that topology.
729 */
730 for (auto iter : elemBlockUniqueFaceTopologies_) {
731 std::string elemBlockName = iter.first;
732 std::vector<stk::topology> uniqueFaceTopologies = iter.second;
733
734 for (auto topo : uniqueFaceTopologies ) {
735 std::string face_block_name = mkBlockName(panzer_stk::STK_Interface::faceBlockString, topo.name());
736 mesh.addFaceBlock(elemBlockName, face_block_name, topo);
737 }
738 }
739 }
740}
741
742bool STK_ExodusReaderFactory::doPerceptRefinement() const
743{
744 return (fileType_!="Exodus Refinement") && (levelsOfRefinement_ > 0);
745}
746
747std::string STK_ExodusReaderFactory::mkBlockName(std::string base, std::string topo_name) const
748{
749 std::string name;
750 name = topo_name+"_"+base;
751 std::transform(name.begin(), name.end(), name.begin(),
752 [](const char c)
753 { return char(std::tolower(c)); });
754 return name;
755}
756
757}
758
759#endif
static const std::string edgeBlockString
static const std::string faceBlockString