Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_Interface.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 <PanzerAdaptersSTK_config.hpp>
12 #include <Panzer_STK_Interface.hpp>
13 
14 #include <Teuchos_as.hpp>
15 #include <Teuchos_RCPStdSharedPtrConversions.hpp>
16 
17 #include <limits>
18 
19 #include <stk_mesh/base/FieldBase.hpp>
20 #include <stk_mesh/base/Comm.hpp>
21 #include <stk_mesh/base/Selector.hpp>
22 #include <stk_mesh/base/GetEntities.hpp>
23 #include <stk_mesh/base/MeshBuilder.hpp>
24 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
25 
26 // #include <stk_rebalance/Rebalance.hpp>
27 // #include <stk_rebalance/Partition.hpp>
28 // #include <stk_rebalance/ZoltanPartition.hpp>
29 // #include <stk_rebalance_utils/RebalanceUtils.hpp>
30 
31 #include <stk_util/parallel/ParallelReduce.hpp>
32 #include <stk_util/parallel/CommSparse.hpp>
33 
34 #include <stk_tools/mesh_tools/FixNodeSharingViaSearch.hpp>
35 
36 #ifdef PANZER_HAVE_IOSS
37 #include <Ionit_Initializer.h>
38 #include <stk_io/IossBridge.hpp>
39 #include <stk_io/WriteMesh.hpp>
40 #endif
41 
42 #ifdef PANZER_HAVE_PERCEPT
43 #include <percept/PerceptMesh.hpp>
44 #include <adapt/UniformRefinerPattern.hpp>
45 #include <adapt/UniformRefiner.hpp>
46 #endif
47 
49 
50 #include <set>
51 
52 using Teuchos::RCP;
53 using Teuchos::rcp;
54 
55 namespace panzer_stk {
56 
58 ElementDescriptor::ElementDescriptor(stk::mesh::EntityId gid,const std::vector<stk::mesh::EntityId> & nodes)
59  : gid_(gid), nodes_(nodes) {}
61 
65 buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes)
66 {
67  return Teuchos::rcp(new ElementDescriptor(elmtId,nodes));
68 }
69 
70 const std::string STK_Interface::coordsString = "coordinates";
71 const std::string STK_Interface::nodesString = "nodes";
72 const std::string STK_Interface::edgesString = "edges";
73 const std::string STK_Interface::facesString = "faces";
74 const std::string STK_Interface::edgeBlockString = "edge_block";
75 const std::string STK_Interface::faceBlockString = "face_block";
76 
78  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
79 {
80  metaData_ = rcp(new stk::mesh::MetaData());
81  metaData_->use_simple_fields();
82 }
83 
85  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
86 {
87  metaData_ = metaData;
88  metaData_->use_simple_fields();
89 }
90 
92  : dimension_(dim), initialized_(false), currentLocalId_(0), useFieldCoordinates_(false)
93 {
94  std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
95  entity_rank_names.push_back("FAMILY_TREE");
96 
97  metaData_ = rcp(new stk::mesh::MetaData(dimension_,entity_rank_names));
98  metaData_->use_simple_fields();
99 
101 }
102 
103 void STK_Interface::addSideset(const std::string & name,const CellTopologyData * ctData)
104 {
107 
108  stk::mesh::Part * sideset = metaData_->get_part(name);
109  if(sideset==nullptr)
110  sideset = &metaData_->declare_part_with_topology(name,
111  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
112  sidesets_.insert(std::make_pair(name,sideset));
113 }
114 
115 void STK_Interface::addNodeset(const std::string & name)
116 {
119 
120  stk::mesh::Part * nodeset = metaData_->get_part(name);
121  if(nodeset==nullptr) {
122  const CellTopologyData * ctData = shards::getCellTopologyData<shards::Node>();
123  nodeset = &metaData_->declare_part_with_topology(name,
124  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
125  }
126  nodesets_.insert(std::make_pair(name,nodeset));
127 }
128 
129 void STK_Interface::addSolutionField(const std::string & fieldName,const std::string & blockId)
130 {
132  "Unknown element block \"" << blockId << "\"");
133  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
134 
135  // add & declare field if not already added...currently assuming linears
136  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
137  SolutionFieldType * field = metaData_->get_field<double>(stk::topology::NODE_RANK, fieldName);
138  if(field==0)
139  field = &metaData_->declare_field<double>(stk::topology::NODE_RANK, fieldName);
140  if ( initialized_ ) {
141  metaData_->enable_late_fields();
142  stk::mesh::put_field_on_mesh(*field, metaData_->universal_part(), nullptr);
143  }
145  }
146 }
147 
148 void STK_Interface::addCellField(const std::string & fieldName,const std::string & blockId)
149 {
151  "Unknown element block \"" << blockId << "\"");
152  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
153 
154  // add & declare field if not already added...currently assuming linears
155  if(fieldNameToCellField_.find(key)==fieldNameToCellField_.end()) {
156  SolutionFieldType * field = metaData_->get_field<double>(stk::topology::ELEMENT_RANK, fieldName);
157  if(field==0)
158  field = &metaData_->declare_field<double>(stk::topology::ELEMENT_RANK, fieldName);
159 
160  if ( initialized_ ) {
161  metaData_->enable_late_fields();
162  stk::mesh::put_field_on_mesh(*field, metaData_->universal_part(), nullptr);
163  }
165  }
166 }
167 
168 void STK_Interface::addEdgeField(const std::string & fieldName,const std::string & blockId)
169 {
171  "Unknown element block \"" << blockId << "\"");
172  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
173 
174  // add & declare field if not already added...currently assuming linears
175  if(fieldNameToEdgeField_.find(key)==fieldNameToEdgeField_.end()) {
176  SolutionFieldType * field = metaData_->get_field<double>(stk::topology::EDGE_RANK, fieldName);
177  if(field==0) {
178  field = &metaData_->declare_field<double>(stk::topology::EDGE_RANK, fieldName);
179  }
180 
181  if ( initialized_ ) {
182  metaData_->enable_late_fields();
183  stk::mesh::put_field_on_mesh(*field, metaData_->universal_part(), nullptr);
184  }
186  }
187 }
188 
189 void STK_Interface::addFaceField(const std::string & fieldName,const std::string & blockId)
190 {
192  "Unknown element block \"" << blockId << "\"");
193  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
194 
195  // add & declare field if not already added...currently assuming linears
196  if(fieldNameToFaceField_.find(key)==fieldNameToFaceField_.end()) {
197  SolutionFieldType * field = metaData_->get_field<double>(stk::topology::FACE_RANK, fieldName);
198  if(field==0) {
199  field = &metaData_->declare_field<double>(stk::topology::FACE_RANK, fieldName);
200  }
201 
202  if ( initialized_ ) {
203  metaData_->enable_late_fields();
204  stk::mesh::put_field_on_mesh(*field, metaData_->universal_part(), nullptr);
205  }
207  }
208 }
209 
210 void STK_Interface::addMeshCoordFields(const std::string & blockId,
211  const std::vector<std::string> & coordNames,
212  const std::string & dispPrefix)
213 {
215  TEUCHOS_ASSERT(dimension_==coordNames.size());
218  "Unknown element block \"" << blockId << "\"");
219 
220  // we only allow one alternative coordinate field
221  TEUCHOS_TEST_FOR_EXCEPTION(meshCoordFields_.find(blockId)!=meshCoordFields_.end(),std::invalid_argument,
222  "STK_Interface::addMeshCoordFields: Can't set more than one set of coordinate "
223  "fields for element block \""+blockId+"\".");
224 
225  // Note that there is a distinction between the key which is used for lookups
226  // and the field that lives on the mesh, which is used for printing the displacement.
227 
228  // just copy the coordinate names
229  meshCoordFields_[blockId] = coordNames;
230 
231  // must fill in the displacement fields
232  std::vector<std::string> & dispFields = meshDispFields_[blockId];
233  dispFields.resize(dimension_);
234 
235  for(unsigned i=0;i<dimension_;i++) {
236  std::pair<std::string,std::string> key = std::make_pair(coordNames[i],blockId);
237  std::string dispName = dispPrefix+coordNames[i];
238 
239  dispFields[i] = dispName; // record this field as a
240  // displacement field
241 
242  // add & declare field if not already added...currently assuming linears
243  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
244 
245  SolutionFieldType * field = metaData_->get_field<double>(stk::topology::NODE_RANK, dispName);
246  if(field==0) {
247  field = &metaData_->declare_field<double>(stk::topology::NODE_RANK, dispName);
248  }
250  }
251  }
252 }
253 
254 void STK_Interface::addInformationRecords(const std::vector<std::string> & info_records)
255 {
256  informationRecords_.insert(info_records.begin(), info_records.end());
257 }
258 
259 void STK_Interface::initialize(stk::ParallelMachine parallelMach,bool setupIO,
260  const bool buildRefinementSupport)
261 {
263  TEUCHOS_ASSERT(dimension_!=0); // no zero dimensional meshes!
264 
265  if(mpiComm_==Teuchos::null)
266  mpiComm_ = getSafeCommunicator(parallelMach);
267 
268  procRank_ = stk::parallel_machine_rank(*mpiComm_->getRawMpiComm());
269 
270  // associating the field with a part: universal part!
271  stk::mesh::put_field_on_mesh( *coordinatesField_ , metaData_->universal_part(), getDimension(), nullptr);
272  stk::mesh::put_field_on_mesh( *edgesField_ , metaData_->universal_part(), getDimension(), nullptr);
273  if (dimension_ > 2)
274  stk::mesh::put_field_on_mesh( *facesField_ , metaData_->universal_part(), getDimension(), nullptr);
275  stk::mesh::put_field_on_mesh( *processorIdField_ , metaData_->universal_part(), nullptr);
276  stk::mesh::put_field_on_mesh( *loadBalField_ , metaData_->universal_part(), nullptr);
277 
282 
283 #ifdef PANZER_HAVE_IOSS
284  if(setupIO) {
285  // setup Exodus file IO
287 
288  // add element blocks
289  {
290  std::map<std::string, stk::mesh::Part*>::iterator itr;
291  for(itr=elementBlocks_.begin();itr!=elementBlocks_.end();++itr)
292  if(!stk::io::is_part_io_part(*itr->second))
293  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
294  }
295 
296  // add edge blocks
297  {
298  std::map<std::string, stk::mesh::Part*>::iterator itr;
299  for(itr=edgeBlocks_.begin();itr!=edgeBlocks_.end();++itr)
300  if(!stk::io::is_part_edge_block_io_part(*itr->second))
301  stk::io::put_edge_block_io_part_attribute(*itr->second); // this can only be called once per part
302  }
303 
304  // add face blocks
305  {
306  std::map<std::string, stk::mesh::Part*>::iterator itr;
307  for(itr=faceBlocks_.begin();itr!=faceBlocks_.end();++itr)
308  if(!stk::io::is_part_face_block_io_part(*itr->second))
309  stk::io::put_face_block_io_part_attribute(*itr->second); // this can only be called once per part
310  }
311 
312  // add side sets
313  {
314  std::map<std::string, stk::mesh::Part*>::iterator itr;
315  for(itr=sidesets_.begin();itr!=sidesets_.end();++itr)
316  if(!stk::io::is_part_io_part(*itr->second))
317  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
318  }
319 
320  // add node sets
321  {
322  std::map<std::string, stk::mesh::Part*>::iterator itr;
323  for(itr=nodesets_.begin();itr!=nodesets_.end();++itr)
324  if(!stk::io::is_part_io_part(*itr->second))
325  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
326  }
327 
328  // add nodes
329  if(!stk::io::is_part_io_part(*nodesPart_))
330  stk::io::put_io_part_attribute(*nodesPart_);
331 
332  stk::io::set_field_role(*coordinatesField_, Ioss::Field::MESH);
333  stk::io::set_field_output_type(*coordinatesField_, stk::io::FieldOutputType::VECTOR_3D);
334  stk::io::set_field_role(*edgesField_, Ioss::Field::MESH);
335  stk::io::set_field_output_type(*edgesField_, stk::io::FieldOutputType::VECTOR_3D);
336  if (dimension_ > 2) {
337  stk::io::set_field_role(*facesField_, Ioss::Field::MESH);
338  stk::io::set_field_output_type(*facesField_, stk::io::FieldOutputType::VECTOR_3D);
339  }
340  stk::io::set_field_role(*processorIdField_, Ioss::Field::TRANSIENT);
341  // stk::io::set_field_role(*loadBalField_, Ioss::Field::TRANSIENT);
342 
343  }
344 #endif
345 
346  if (buildRefinementSupport) {
347 #ifdef PANZER_HAVE_PERCEPT
348  refinedMesh_ = Teuchos::rcp(new percept::PerceptMesh(this->getMetaData().get(),
349  this->getBulkData().get(),
350  true));
351 
352  breakPattern_ = Teuchos::rcp(new percept::URP_Heterogeneous_3D(*refinedMesh_));
353 #else
354  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
355  "ERROR: Uniform refinement requested. This requires the Percept package to be enabled in Trilinos!");
356 #endif
357  }
358 
359  if(bulkData_==Teuchos::null)
360  instantiateBulkData(*mpiComm_->getRawMpiComm());
361 
362  metaData_->commit();
363 
364  initialized_ = true;
365 }
366 
367 void STK_Interface::initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
368  bool setupIO)
369 {
370  std::set<SolutionFieldType*> uniqueFields;
371  std::map<std::pair<std::string,std::string>,SolutionFieldType*>::const_iterator fieldIter;
372  for(fieldIter=nameToField.begin();fieldIter!=nameToField.end();++fieldIter)
373  uniqueFields.insert(fieldIter->second); // this makes setting up IO easier!
374 
375  {
376  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
377  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
378  stk::mesh::put_field_on_mesh(*(*uniqueFieldIter),metaData_->universal_part(),nullptr);
379  }
380 
381 #ifdef PANZER_HAVE_IOSS
382  if(setupIO) {
383  // add solution fields
384  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
385  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
386  stk::io::set_field_role(*(*uniqueFieldIter), Ioss::Field::TRANSIENT);
387  }
388 #endif
389 }
390 
391 void STK_Interface::instantiateBulkData(stk::ParallelMachine parallelMach)
392 {
393  TEUCHOS_ASSERT(bulkData_==Teuchos::null);
394  if(mpiComm_==Teuchos::null)
395  mpiComm_ = getSafeCommunicator(parallelMach);
396 
397  std::unique_ptr<stk::mesh::BulkData> bulkUPtr = stk::mesh::MeshBuilder(*mpiComm_->getRawMpiComm()).create(Teuchos::get_shared_ptr(metaData_));
398  bulkData_ = rcp(bulkUPtr.release());
399 }
400 
402 {
403  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
404  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"beginModification\"");
405 
406  bulkData_->modification_begin();
407 }
408 
409 void STK_Interface::endModification(const bool find_and_set_shared_nodes_in_stk)
410 {
411  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
412  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"endModification\"");
413 
414  // Resolving sharing this way comes at a cost in performance. The
415  // STK mesh database requires users to declare their own node
416  // sharing. We should find where shared entities are being created
417  // in Panzer (the inline mesh factories) and declare it. Once this
418  // is done, the extra code below can be deleted. Note that the
419  // exodus reader already has shared nodes set up properly, so only
420  // the inline mesh factories need this.
421  if (find_and_set_shared_nodes_in_stk) {
422  const double radius = 10.0 * std::numeric_limits<double>::epsilon();
423  stk::tools::fix_node_sharing_via_search(*bulkData_,radius);
424  }
425 
426  bulkData_->modification_end();
427 
430 }
431 
432 void STK_Interface::addNode(stk::mesh::EntityId gid, const std::vector<double> & coord)
433 {
434  TEUCHOS_TEST_FOR_EXCEPTION(not isModifiable(),std::logic_error,
435  "STK_Interface::addNode: STK_Interface must be modifiable to add a node");
436  TEUCHOS_TEST_FOR_EXCEPTION(coord.size()!=getDimension(),std::logic_error,
437  "STK_Interface::addNode: number of coordinates in vector must mation dimension");
438  TEUCHOS_TEST_FOR_EXCEPTION(gid==0,std::logic_error,
439  "STK_Interface::addNode: STK has STUPID restriction of no zero GIDs, pick something else");
440  stk::mesh::EntityRank nodeRank = getNodeRank();
441 
442  stk::mesh::Entity node = bulkData_->declare_entity(nodeRank,gid,nodesPartVec_);
443 
444  // set coordinate vector
445  double * fieldCoords = stk::mesh::field_data(*coordinatesField_,node);
446  for(std::size_t i=0;i<coord.size();++i)
447  fieldCoords[i] = coord[i];
448 }
449 
450 void STK_Interface::addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset)
451 {
452  std::vector<stk::mesh::Part*> sidesetV;
453  sidesetV.push_back(sideset);
454 
455  bulkData_->change_entity_parts(entity,sidesetV);
456 }
457 
458 void STK_Interface::addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset)
459 {
460  std::vector<stk::mesh::Part*> nodesetV;
461  nodesetV.push_back(nodeset);
462 
463  bulkData_->change_entity_parts(entity,nodesetV);
464 }
465 
466 void STK_Interface::addEntityToEdgeBlock(stk::mesh::Entity entity,stk::mesh::Part * edgeblock)
467 {
468  std::vector<stk::mesh::Part*> edgeblockV;
469  edgeblockV.push_back(edgeblock);
470 
471  bulkData_->change_entity_parts(entity,edgeblockV);
472 }
473 void STK_Interface::addEntitiesToEdgeBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * edgeblock)
474 {
475  if (entities.size() > 0) {
476  std::vector<stk::mesh::Part*> edgeblockV;
477  edgeblockV.push_back(edgeblock);
478 
479  bulkData_->change_entity_parts(entities,edgeblockV);
480  }
481 }
482 
483 void STK_Interface::addEntityToFaceBlock(stk::mesh::Entity entity,stk::mesh::Part * faceblock)
484 {
485  std::vector<stk::mesh::Part*> faceblockV;
486  faceblockV.push_back(faceblock);
487 
488  bulkData_->change_entity_parts(entity,faceblockV);
489 }
490 void STK_Interface::addEntitiesToFaceBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * faceblock)
491 {
492  if (entities.size() > 0) {
493  std::vector<stk::mesh::Part*> faceblockV;
494  faceblockV.push_back(faceblock);
495 
496  bulkData_->change_entity_parts(entities,faceblockV);
497  }
498 }
499 
500 void STK_Interface::addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block)
501 {
502  std::vector<stk::mesh::Part*> blockVec;
503  blockVec.push_back(block);
504 
505  stk::mesh::EntityRank elementRank = getElementRank();
506  stk::mesh::EntityRank nodeRank = getNodeRank();
507  stk::mesh::Entity element = bulkData_->declare_entity(elementRank,ed->getGID(),blockVec);
508 
509  // build relations that give the mesh structure
510  const std::vector<stk::mesh::EntityId> & nodes = ed->getNodes();
511  for(std::size_t i=0;i<nodes.size();++i) {
512  // add element->node relation
513  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodes[i]);
514  TEUCHOS_ASSERT(bulkData_->is_valid(node));
515  bulkData_->declare_relation(element,node,i);
516  }
517 
518  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
519  procId[0] = Teuchos::as<ProcIdData>(procRank_);
520 }
521 
522 
524 {
525  // loop over elements
526  stk::mesh::EntityRank edgeRank = getEdgeRank();
527  std::vector<stk::mesh::Entity> localElmts;
528  getMyElements(localElmts);
529  std::vector<stk::mesh::Entity>::const_iterator itr;
530  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
531  stk::mesh::Entity element = (*itr);
532  stk::mesh::EntityId gid = bulkData_->identifier(element);
533  std::vector<stk::mesh::EntityId> subcellIds;
534  getSubcellIndices(edgeRank,gid,subcellIds);
535 
536  for(std::size_t i=0;i<subcellIds.size();++i) {
537  stk::mesh::Entity edge = bulkData_->get_entity(edgeRank,subcellIds[i]);
538  stk::mesh::Entity const* relations = bulkData_->begin_nodes(edge);
539 
540  double * node_coord_1 = stk::mesh::field_data(*coordinatesField_,relations[0]);
541  double * node_coord_2 = stk::mesh::field_data(*coordinatesField_,relations[1]);
542 
543  // set coordinate vector
544  double * edgeCoords = stk::mesh::field_data(*edgesField_,edge);
545  for(std::size_t j=0;j<getDimension();++j)
546  edgeCoords[j] = (node_coord_1[j]+node_coord_2[j])/2.0;
547  }
548  }
549 }
550 
552 {
553  // loop over elements
554  stk::mesh::EntityRank faceRank = getFaceRank();
555  std::vector<stk::mesh::Entity> localElmts;
556  getMyElements(localElmts);
557  std::vector<stk::mesh::Entity>::const_iterator itr;
558  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
559  stk::mesh::Entity element = (*itr);
560  stk::mesh::EntityId gid = elementGlobalId(element);
561  std::vector<stk::mesh::EntityId> subcellIds;
562  getSubcellIndices(faceRank,gid,subcellIds);
563 
564  for(std::size_t i=0;i<subcellIds.size();++i) {
565  stk::mesh::Entity face = bulkData_->get_entity(faceRank,subcellIds[i]);
566  stk::mesh::Entity const* relations = bulkData_->begin_nodes(face);
567  const size_t num_relations = bulkData_->num_nodes(face);
568 
569  // set coordinate vector
570  double * faceCoords = stk::mesh::field_data(*facesField_,face);
571  for(std::size_t j=0;j<getDimension();++j){
572  faceCoords[j] = 0.0;
573  for(std::size_t k=0;k<num_relations;++k)
574  faceCoords[j] += stk::mesh::field_data(*coordinatesField_,relations[k])[j];
575  faceCoords[j] /= double(num_relations);
576  }
577  }
578  }
579 }
580 
581 stk::mesh::Entity STK_Interface::findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
582 {
583  const size_t num_rels = bulkData_->num_connectivity(src, tgt_rank);
584  stk::mesh::Entity const* relations = bulkData_->begin(src, tgt_rank);
585  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData_->begin_ordinals(src, tgt_rank);
586  for (size_t i = 0; i < num_rels; ++i) {
587  if (ordinals[i] == static_cast<stk::mesh::ConnectivityOrdinal>(rel_id)) {
588  return relations[i];
589  }
590  }
591 
592  return stk::mesh::Entity();
593 }
594 
596 //
597 // writeToExodus()
598 //
600 void
602 writeToExodus(const std::string& filename,
603  const bool append)
604 {
605  setupExodusFile(filename,append);
606  writeToExodus(0.0);
607 } // end of writeToExodus()
608 
610 //
611 // setupExodusFile()
612 //
614 void
616 setupExodusFile(const std::string& filename,
617  const bool append,
618  const bool append_after_restart_time,
619  const double restart_time)
620 {
621  std::vector<Ioss::Property> ioss_properties;
622  setupExodusFile(filename,
623  ioss_properties,
624  append,
625  append_after_restart_time,
626  restart_time);
627 }
628 
629 void
631 setupExodusFile(const std::string& filename,
632  const std::vector<Ioss::Property>& ioss_properties,
633  const bool append,
634  const bool append_after_restart_time,
635  const double restart_time)
636 {
637  using std::runtime_error;
638  using stk::io::StkMeshIoBroker;
639  using stk::mesh::FieldVector;
640  using stk::ParallelMachine;
641  using Teuchos::rcp;
642  PANZER_FUNC_TIME_MONITOR("STK_Interface::setupExodusFile(filename)");
643 #ifdef PANZER_HAVE_IOSS
645 
646  ParallelMachine comm = *mpiComm_->getRawMpiComm();
647  meshData_ = rcp(new StkMeshIoBroker(comm));
648  meshData_->set_bulk_data(Teuchos::get_shared_ptr(bulkData_));
649  Ioss::PropertyManager props;
650  props.add(Ioss::Property("LOWER_CASE_VARIABLE_NAMES", "FALSE"));
651  for ( const auto& p : ioss_properties ) {
652  props.add(p);
653  }
654  if (append) {
655  if (append_after_restart_time) {
656  meshIndex_ = meshData_->create_output_mesh(filename, stk::io::APPEND_RESULTS,
657  props, restart_time);
658  }
659  else // Append results to the end of the file
660  meshIndex_ = meshData_->create_output_mesh(filename, stk::io::APPEND_RESULTS, props);
661  }
662  else
663  meshIndex_ = meshData_->create_output_mesh(filename, stk::io::WRITE_RESULTS, props);
664  const FieldVector& fields = metaData_->get_fields();
665  for (size_t i(0); i < fields.size(); ++i) {
666  // Do NOT add MESH type stk fields to exodus io, but do add everything
667  // else. This allows us to avoid having to catch a throw for
668  // re-registering coordinates, sidesets, etc... Note that some
669  // fields like LOAD_BAL don't always have a role assigned, so for
670  // roles that point to null, we need to register them as well.
671  auto role = stk::io::get_field_role(*fields[i]);
672  if (role != nullptr) {
673  if (*role != Ioss::Field::MESH)
674  meshData_->add_field(meshIndex_, *fields[i]);
675  } else {
676  meshData_->add_field(meshIndex_, *fields[i]);
677  }
678  }
679 
680  // convert the set to a vector
681  std::vector<std::string> deduped_info_records(informationRecords_.begin(),informationRecords_.end());
682  meshData_->add_info_records(meshIndex_, deduped_info_records);
683 #else
684  TEUCHOS_ASSERT(false)
685 #endif
686 } // end of setupExodusFile()
687 
689 //
690 // writeToExodus()
691 //
693 void
696  double timestep)
697 {
698  using std::cout;
699  using std::endl;
700  using Teuchos::FancyOStream;
701  using Teuchos::rcpFromRef;
702  PANZER_FUNC_TIME_MONITOR("STK_Interface::writeToExodus(timestep)");
703 #ifdef PANZER_HAVE_IOSS
704  if (not meshData_.is_null())
705  {
706  currentStateTime_ = timestep;
707  globalToExodus(GlobalVariable::ADD);
708  meshData_->begin_output_step(meshIndex_, currentStateTime_);
709  meshData_->write_defined_output_fields(meshIndex_);
710  globalToExodus(GlobalVariable::WRITE);
711  meshData_->end_output_step(meshIndex_);
712  }
713  else // if (meshData_.is_null())
714  {
715  FancyOStream out(rcpFromRef(cout));
716  out.setOutputToRootOnly(0);
717  out << "WARNING: Exodus I/O has been disabled or not setup properly; "
718  << "not writing to Exodus." << endl;
719  } // end if (meshData_.is_null()) or not
720 #else
721  TEUCHOS_ASSERT(false);
722 #endif
723 } // end of writeToExodus()
724 
726 //
727 // globalToExodus()
728 //
730 void
731 STK_Interface::
732 globalToExodus(
733  const GlobalVariable& flag)
734 {
735  using std::invalid_argument;
736  using std::string;
737  using Teuchos::Array;
738 
739  // Loop over all the global variables to be added to the Exodus output file.
740  // For each global variable, we determine the data type, and then add or
741  // write it accordingly, depending on the value of flag.
742  for (auto i = globalData_.begin(); i != globalData_.end(); ++i)
743  {
744  const string& name = globalData_.name(i);
745 
746  // Integers.
747  if (globalData_.isType<int>(name))
748  {
749  const auto& value = globalData_.get<int>(name);
750  if (flag == GlobalVariable::ADD)
751  {
752  try
753  {
754  meshData_->add_global(meshIndex_, name, value,
755  stk::util::ParameterType::INTEGER);
756  }
757  catch (...)
758  {
759  return;
760  }
761  }
762  else // if (flag == GlobalVariable::WRITE)
763  meshData_->write_global(meshIndex_, name, value,
764  stk::util::ParameterType::INTEGER);
765  }
766 
767  // Doubles.
768  else if (globalData_.isType<double>(name))
769  {
770  const auto& value = globalData_.get<double>(name);
771  if (flag == GlobalVariable::ADD)
772  {
773  try
774  {
775  meshData_->add_global(meshIndex_, name, value,
776  stk::util::ParameterType::DOUBLE);
777  }
778  catch (...)
779  {
780  return;
781  }
782  }
783  else // if (flag == GlobalVariable::WRITE)
784  meshData_->write_global(meshIndex_, name, value,
785  stk::util::ParameterType::DOUBLE);
786  }
787 
788  // Vectors of integers.
789  else if (globalData_.isType<Array<int>>(name))
790  {
791  const auto& value = globalData_.get<Array<int>>(name).toVector();
792  if (flag == GlobalVariable::ADD)
793  {
794  try
795  {
796  meshData_->add_global(meshIndex_, name, value,
797  stk::util::ParameterType::INTEGERVECTOR);
798  }
799  catch (...)
800  {
801  return;
802  }
803  }
804  else // if (flag == GlobalVariable::WRITE)
805  meshData_->write_global(meshIndex_, name, value,
806  stk::util::ParameterType::INTEGERVECTOR);
807  }
808 
809  // Vectors of doubles.
810  else if (globalData_.isType<Array<double>>(name))
811  {
812  const auto& value = globalData_.get<Array<double>>(name).toVector();
813  if (flag == GlobalVariable::ADD)
814  {
815  try
816  {
817  meshData_->add_global(meshIndex_, name, value,
818  stk::util::ParameterType::DOUBLEVECTOR);
819  }
820  catch (...)
821  {
822  return;
823  }
824  }
825  else // if (flag == GlobalVariable::WRITE)
826  meshData_->write_global(meshIndex_, name, value,
827  stk::util::ParameterType::DOUBLEVECTOR);
828  }
829 
830  // If the data type is something else, throw an exception.
831  else
832  TEUCHOS_TEST_FOR_EXCEPTION(true, invalid_argument,
833  "STK_Interface::globalToExodus(): The global variable to be added " \
834  "to the Exodus output file is of an invalid type. Valid types are " \
835  "int and double, along with std::vectors of those types.")
836  } // end loop over globalData_
837 } // end of globalToExodus()
838 
840 //
841 // addGlobalToExodus()
842 //
844 void
847  const std::string& key,
848  const int& value)
849 {
850  globalData_.set(key, value);
851 } // end of addGlobalToExodus()
852 
854 //
855 // addGlobalToExodus()
856 //
858 void
861  const std::string& key,
862  const double& value)
863 {
864  globalData_.set(key, value);
865 } // end of addGlobalToExodus()
866 
868 //
869 // addGlobalToExodus()
870 //
872 void
875  const std::string& key,
876  const std::vector<int>& value)
877 {
878  using Teuchos::Array;
879  globalData_.set(key, Array<int>(value));
880 } // end of addGlobalToExodus()
881 
883 //
884 // addGlobalToExodus()
885 //
887 void
890  const std::string& key,
891  const std::vector<double>& value)
892 {
893  using Teuchos::Array;
894  globalData_.set(key, Array<double>(value));
895 } // end of addGlobalToExodus()
896 
898 {
899  #ifdef PANZER_HAVE_IOSS
900  return true;
901  #else
902  return false;
903  #endif
904 }
905 
906 void STK_Interface::getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const
907 {
908  stk::mesh::EntityRank nodeRank = getNodeRank();
909 
910  // get all relations for node
911  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodeId);
912  const size_t numElements = bulkData_->num_elements(node);
913  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
914 
915  // extract elements sharing nodes
916  elements.insert(elements.end(), relations, relations + numElements);
917 }
918 
919 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
920  std::vector<int> & relIds) const
921 {
922  // get all relations for node
923  const size_t numElements = bulkData_->num_elements(node);
924  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
925  stk::mesh::ConnectivityOrdinal const* rel_ids = bulkData_->begin_element_ordinals(node);
926 
927  // extract elements sharing nodes
928  for (size_t i = 0; i < numElements; ++i) {
929  stk::mesh::Entity element = relations[i];
930 
931  // if owned by this processor
932  if(bulkData_->parallel_owner_rank(element) == static_cast<int>(procRank_)) {
933  elements.push_back(element);
934  relIds.push_back(rel_ids[i]);
935  }
936  }
937 }
938 
939 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
940  std::vector<int> & relIds, unsigned int matchType) const
941 {
942  stk::mesh::EntityRank rank;
943  if(matchType == 0)
944  rank = getNodeRank();
945  else if(matchType == 1)
946  rank = getEdgeRank();
947  else if(matchType == 2)
948  rank = getFaceRank();
949  else
950  TEUCHOS_ASSERT(false);
951 
952  stk::mesh::Entity node = bulkData_->get_entity(rank,nodeId);
953 
954  getOwnedElementsSharingNode(node,elements,relIds);
955 }
956 
957 void STK_Interface::getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeIds,std::vector<stk::mesh::Entity> & elements) const
958 {
959  std::vector<stk::mesh::Entity> current;
960 
961  getElementsSharingNode(nodeIds[0],current); // fill it with elements touching first node
962  std::sort(current.begin(),current.end()); // sort for intersection on the pointer
963 
964  // find intersection with remaining nodes
965  for(std::size_t n=1;n<nodeIds.size();++n) {
966  // get elements associated with next node
967  std::vector<stk::mesh::Entity> nextNode;
968  getElementsSharingNode(nodeIds[n],nextNode); // fill it with elements touching first node
969  std::sort(nextNode.begin(),nextNode.end()); // sort for intersection on the pointer ID
970 
971  // intersect next node elements with current element list
972  std::vector<stk::mesh::Entity> intersection(std::min(nextNode.size(),current.size()));
973  std::vector<stk::mesh::Entity>::const_iterator endItr
974  = std::set_intersection(current.begin(),current.end(),
975  nextNode.begin(),nextNode.end(),
976  intersection.begin());
977  std::size_t newLength = endItr-intersection.begin();
978  intersection.resize(newLength);
979 
980  // store intersection
981  current.clear();
982  current = intersection;
983  }
984 
985  // return the elements computed
986  elements = current;
987 }
988 
989 void STK_Interface::getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const
990 {
991  stk::mesh::Entity const* nodeRel = getBulkData()->begin_nodes(element);
992  const size_t numNodes = getBulkData()->num_nodes(element);
993 
994  nodeIds.reserve(numNodes);
995  for(size_t i = 0; i < numNodes; ++i) {
996  nodeIds.push_back(elementGlobalId(nodeRel[i]));
997  }
998 }
999 
1001 {
1002  entityCounts_.clear();
1003  stk::mesh::comm_mesh_counts(*bulkData_,entityCounts_);
1004 }
1005 
1007 {
1008  // developed to mirror "comm_mesh_counts" in stk_mesh/base/Comm.cpp
1009 
1010  const auto entityRankCount = metaData_->entity_rank_count();
1011  const size_t commCount = 10; // entityRankCount
1012 
1013  TEUCHOS_ASSERT(entityRankCount<10);
1014 
1015  // stk::ParallelMachine mach = bulkData_->parallel();
1016  stk::ParallelMachine mach = *mpiComm_->getRawMpiComm();
1017 
1018  std::vector<stk::mesh::EntityId> local(commCount,0);
1019 
1020  // determine maximum ID for this processor for each entity type
1021  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
1022  for(stk::mesh::EntityRank i=stk::topology::NODE_RANK;
1023  i < static_cast<stk::mesh::EntityRank>(entityRankCount); ++i) {
1024  std::vector<stk::mesh::Entity> entities;
1025 
1026  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(i),entities);
1027 
1028  // determine maximum ID for this processor
1029  std::vector<stk::mesh::Entity>::const_iterator itr;
1030  for(itr=entities.begin();itr!=entities.end();++itr) {
1031  stk::mesh::EntityId id = bulkData_->identifier(*itr);
1032  if(id>local[i])
1033  local[i] = id;
1034  }
1035  }
1036 
1037  // get largest IDs across processors
1038  stk::all_reduce(mach,stk::ReduceMax<10>(&local[0]));
1039  maxEntityId_.assign(local.begin(),local.begin()+entityRankCount+1);
1040 }
1041 
1042 std::size_t STK_Interface::getEntityCounts(unsigned entityRank) const
1043 {
1044  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=entityCounts_.size(),std::logic_error,
1045  "STK_Interface::getEntityCounts: Entity counts do not include rank: " << entityRank);
1046 
1047  return entityCounts_[entityRank];
1048 }
1049 
1050 stk::mesh::EntityId STK_Interface::getMaxEntityId(unsigned entityRank) const
1051 {
1052  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=maxEntityId_.size(),std::logic_error,
1053  "STK_Interface::getMaxEntityId: Max entity ids do not include rank: " << entityRank);
1054 
1055  return maxEntityId_[entityRank];
1056 }
1057 
1059 {
1060  stk::mesh::PartVector emptyPartVector;
1061  stk::mesh::create_adjacent_entities(*bulkData_,emptyPartVector);
1062 
1065 
1066  addEdges();
1067  if (dimension_ > 2)
1068  addFaces();
1069 }
1070 
1071 const double * STK_Interface::getNodeCoordinates(stk::mesh::EntityId nodeId) const
1072 {
1073  stk::mesh::Entity node = bulkData_->get_entity(getNodeRank(),nodeId);
1074  return stk::mesh::field_data(*coordinatesField_,node);
1075 }
1076 
1077 const double * STK_Interface::getNodeCoordinates(stk::mesh::Entity node) const
1078 {
1079  return stk::mesh::field_data(*coordinatesField_,node);
1080 }
1081 
1082 void STK_Interface::getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
1083  std::vector<stk::mesh::EntityId> & subcellIds) const
1084 {
1085  stk::mesh::EntityRank elementRank = getElementRank();
1086  stk::mesh::Entity cell = bulkData_->get_entity(elementRank,elementId);
1087 
1088  TEUCHOS_TEST_FOR_EXCEPTION(!bulkData_->is_valid(cell),std::logic_error,
1089  "STK_Interface::getSubcellIndices: could not find element requested (GID = " << elementId << ")");
1090 
1091  const size_t numSubcells = bulkData_->num_connectivity(cell, static_cast<stk::mesh::EntityRank>(entityRank));
1092  stk::mesh::Entity const* subcells = bulkData_->begin(cell, static_cast<stk::mesh::EntityRank>(entityRank));
1093  subcellIds.clear();
1094  subcellIds.resize(numSubcells,0);
1095 
1096  // loop over relations and fill subcell vector
1097  for(size_t i = 0; i < numSubcells; ++i) {
1098  stk::mesh::Entity subcell = subcells[i];
1099  subcellIds[i] = bulkData_->identifier(subcell);
1100  }
1101 }
1102 
1103 void STK_Interface::getMyElements(std::vector<stk::mesh::Entity> & elements) const
1104 {
1105  // setup local ownership
1106  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
1107 
1108  // grab elements
1109  stk::mesh::EntityRank elementRank = getElementRank();
1110  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(elementRank),elements);
1111 }
1112 
1113 void STK_Interface::getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
1114 {
1115  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
1116 
1117  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
1118 
1119  // setup local ownership
1120  // stk::mesh::Selector block = *elementBlock;
1121  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & (*elementBlock);
1122 
1123  // grab elements
1124  stk::mesh::EntityRank elementRank = getElementRank();
1125  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(elementRank),elements);
1126 }
1127 
1128 void STK_Interface::getNeighborElements(std::vector<stk::mesh::Entity> & elements) const
1129 {
1130  // setup local ownership
1131  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part());
1132 
1133  // grab elements
1134  stk::mesh::EntityRank elementRank = getElementRank();
1135  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
1136 }
1137 
1138 void STK_Interface::getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
1139 {
1140  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
1141 
1142  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
1143 
1144  // setup local ownership
1145  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part()) & (*elementBlock);
1146 
1147  // grab elements
1148  stk::mesh::EntityRank elementRank = getElementRank();
1149  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
1150 }
1151 
1152 void STK_Interface::getMyEdges(std::vector<stk::mesh::Entity> & edges) const
1153 {
1154  // setup local ownership
1155  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
1156 
1157  // grab elements
1158  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(getEdgeRank()),edges);
1159 }
1160 
1161 void STK_Interface::getMyEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const
1162 {
1163  stk::mesh::Part * edgeBlockPart = getEdgeBlock(edgeBlockName);
1164  TEUCHOS_TEST_FOR_EXCEPTION(edgeBlockPart==0,std::logic_error,
1165  "Unknown edge block \"" << edgeBlockName << "\"");
1166 
1167  stk::mesh::Selector edge_block = *edgeBlockPart;
1168  stk::mesh::Selector owned_block = metaData_->locally_owned_part() & edge_block;
1169 
1170  // grab elements
1171  stk::mesh::get_selected_entities(owned_block,bulkData_->buckets(getEdgeRank()),edges);
1172 }
1173 
1174 void STK_Interface::getMyEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const
1175 {
1176  stk::mesh::Part * edgeBlockPart = getEdgeBlock(edgeBlockName);
1177  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1179  "Unknown edge block \"" << edgeBlockName << "\"");
1181  "Unknown element block \"" << blockName << "\"");
1182 
1183  stk::mesh::Selector edge_block = *edgeBlockPart;
1184  stk::mesh::Selector element_block = *elmtPart;
1185  stk::mesh::Selector owned_block = metaData_->locally_owned_part() & element_block & edge_block;
1186 
1187  // grab elements
1188  stk::mesh::get_selected_entities(owned_block,bulkData_->buckets(getEdgeRank()),edges);
1189 }
1190 
1191 void STK_Interface::getAllEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const
1192 {
1193  stk::mesh::Part * edgeBlockPart = getEdgeBlock(edgeBlockName);
1194  TEUCHOS_TEST_FOR_EXCEPTION(edgeBlockPart==0,std::logic_error,
1195  "Unknown edge block \"" << edgeBlockName << "\"");
1196 
1197  stk::mesh::Selector edge_block = *edgeBlockPart;
1198 
1199  // grab elements
1200  stk::mesh::get_selected_entities(edge_block,bulkData_->buckets(getEdgeRank()),edges);
1201 }
1202 
1203 void STK_Interface::getAllEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const
1204 {
1205  stk::mesh::Part * edgeBlockPart = getEdgeBlock(edgeBlockName);
1206  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1208  "Unknown edge block \"" << edgeBlockName << "\"");
1210  "Unknown element block \"" << blockName << "\"");
1211 
1212  stk::mesh::Selector edge_block = *edgeBlockPart;
1213  stk::mesh::Selector element_block = *elmtPart;
1214  stk::mesh::Selector element_edge_block = element_block & edge_block;
1215 
1216  // grab elements
1217  stk::mesh::get_selected_entities(element_edge_block,bulkData_->buckets(getEdgeRank()),edges);
1218 }
1219 
1220 void STK_Interface::getMyFaces(std::vector<stk::mesh::Entity> & faces) const
1221 {
1222  // setup local ownership
1223  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
1224 
1225  // grab elements
1226  stk::mesh::EntityRank faceRank = getFaceRank();
1227  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(faceRank),faces);
1228 }
1229 
1230 void STK_Interface::getMyFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const
1231 {
1232  stk::mesh::Part * faceBlockPart = getFaceBlock(faceBlockName);
1233  TEUCHOS_TEST_FOR_EXCEPTION(faceBlockPart==0,std::logic_error,
1234  "Unknown face block \"" << faceBlockName << "\"");
1235 
1236  stk::mesh::Selector face_block = *faceBlockPart;
1237  stk::mesh::Selector owned_block = metaData_->locally_owned_part() & face_block;
1238 
1239  // grab elements
1240  stk::mesh::get_selected_entities(owned_block,bulkData_->buckets(getFaceRank()),faces);
1241 }
1242 
1243 void STK_Interface::getMyFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const
1244 {
1245  stk::mesh::Part * faceBlockPart = getFaceBlock(faceBlockName);
1246  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1248  "Unknown face block \"" << faceBlockName << "\"");
1250  "Unknown element block \"" << blockName << "\"");
1251 
1252  stk::mesh::Selector face_block = *faceBlockPart;
1253  stk::mesh::Selector element_block = *elmtPart;
1254  stk::mesh::Selector owned_block = metaData_->locally_owned_part() & element_block & face_block;
1255 
1256  // grab elements
1257  stk::mesh::get_selected_entities(owned_block,bulkData_->buckets(getFaceRank()),faces);
1258 }
1259 
1260 void STK_Interface::getAllFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const
1261 {
1262  stk::mesh::Part * faceBlockPart = getFaceBlock(faceBlockName);
1263  TEUCHOS_TEST_FOR_EXCEPTION(faceBlockPart==0,std::logic_error,
1264  "Unknown face block \"" << faceBlockName << "\"");
1265 
1266  stk::mesh::Selector face_block = *faceBlockPart;
1267 
1268  // grab elements
1269  stk::mesh::get_selected_entities(face_block,bulkData_->buckets(getFaceRank()),faces);
1270 }
1271 
1272 void STK_Interface::getAllFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const
1273 {
1274  stk::mesh::Part * faceBlockPart = getFaceBlock(faceBlockName);
1275  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1277  "Unknown face block \"" << faceBlockName << "\"");
1279  "Unknown element block \"" << blockName << "\"");
1280 
1281  stk::mesh::Selector face_block = *faceBlockPart;
1282  stk::mesh::Selector element_block = *elmtPart;
1283  stk::mesh::Selector element_face_block = element_block & face_block;
1284 
1285  // grab elements
1286  stk::mesh::get_selected_entities(element_face_block,bulkData_->buckets(getFaceRank()),faces);
1287 }
1288 
1289 void STK_Interface::getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
1290 {
1291  stk::mesh::Part * sidePart = getSideset(sideName);
1292  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
1293  "Unknown side set \"" << sideName << "\"");
1294 
1295  stk::mesh::Selector side = *sidePart;
1296  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & side;
1297 
1298  // grab elements
1299  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
1300 }
1301 
1302 void STK_Interface::getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
1303 {
1304  stk::mesh::Part * sidePart = getSideset(sideName);
1305  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1307  "Unknown side set \"" << sideName << "\"");
1309  "Unknown element block \"" << blockName << "\"");
1310 
1311  stk::mesh::Selector side = *sidePart;
1312  stk::mesh::Selector block = *elmtPart;
1313  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & side;
1314 
1315  // grab elements
1316  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
1317 }
1318 
1319 void STK_Interface::getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
1320 {
1321  stk::mesh::Part * sidePart = getSideset(sideName);
1322  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
1323  "Unknown side set \"" << sideName << "\"");
1324 
1325  stk::mesh::Selector side = *sidePart;
1326 
1327  // grab elements
1328  stk::mesh::get_selected_entities(side,bulkData_->buckets(getSideRank()),sides);
1329 }
1330 
1331 void STK_Interface::getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
1332 {
1333  stk::mesh::Part * sidePart = getSideset(sideName);
1334  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1336  "Unknown side set \"" << sideName << "\"");
1338  "Unknown element block \"" << blockName << "\"");
1339 
1340  stk::mesh::Selector side = *sidePart;
1341  stk::mesh::Selector block = *elmtPart;
1342  stk::mesh::Selector sideBlock = block & side;
1343 
1344  // grab elements
1345  stk::mesh::get_selected_entities(sideBlock,bulkData_->buckets(getSideRank()),sides);
1346 }
1347 
1348 void STK_Interface::getMyNodes(const std::string & nodesetName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const
1349 {
1350  stk::mesh::Part * nodePart = getNodeset(nodesetName);
1351  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1353  "Unknown node set \"" << nodesetName << "\"");
1355  "Unknown element block \"" << blockName << "\"");
1356 
1357  stk::mesh::Selector nodeset = *nodePart;
1358  stk::mesh::Selector block = *elmtPart;
1359  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & nodeset;
1360 
1361  // grab elements
1362  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getNodeRank()),nodes);
1363 }
1364 
1365 void STK_Interface::getElementBlockNames(std::vector<std::string> & names) const
1366 {
1367  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
1368 
1369  names.clear();
1370 
1371  // fill vector with automagically ordered string values
1372  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr; // Element blocks
1373  for(blkItr=elementBlocks_.begin();blkItr!=elementBlocks_.end();++blkItr)
1374  names.push_back(blkItr->first);
1375 }
1376 
1377 void STK_Interface::getSidesetNames(std::vector<std::string> & names) const
1378 {
1379  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
1380 
1381  names.clear();
1382 
1383  // fill vector with automagically ordered string values
1384  std::map<std::string, stk::mesh::Part*>::const_iterator sideItr; // Side sets
1385  for(sideItr=sidesets_.begin();sideItr!=sidesets_.end();++sideItr)
1386  names.push_back(sideItr->first);
1387 }
1388 
1389 void STK_Interface::getNodesetNames(std::vector<std::string> & names) const
1390 {
1391  names.clear();
1392 
1393  // fill vector with automagically ordered string values
1394  std::map<std::string, stk::mesh::Part*>::const_iterator nodeItr; // Node sets
1395  for(nodeItr=nodesets_.begin();nodeItr!=nodesets_.end();++nodeItr)
1396  names.push_back(nodeItr->first);
1397 }
1398 
1399 void STK_Interface::getEdgeBlockNames(std::vector<std::string> & names) const
1400 {
1401  names.clear();
1402 
1403  // fill vector with automagically ordered string values
1404  std::map<std::string, stk::mesh::Part*>::const_iterator edgeBlockItr; // Edge blocks
1405  for(edgeBlockItr=edgeBlocks_.begin();edgeBlockItr!=edgeBlocks_.end();++edgeBlockItr)
1406  names.push_back(edgeBlockItr->first);
1407 }
1408 
1409 void STK_Interface::getFaceBlockNames(std::vector<std::string> & names) const
1410 {
1411  names.clear();
1412 
1413  // fill vector with automagically ordered string values
1414  std::map<std::string, stk::mesh::Part*>::const_iterator faceBlockItr; // Face blocks
1415  for(faceBlockItr=faceBlocks_.begin();faceBlockItr!=faceBlocks_.end();++faceBlockItr)
1416  names.push_back(faceBlockItr->first);
1417 }
1418 
1419 std::size_t STK_Interface::elementLocalId(stk::mesh::Entity elmt) const
1420 {
1421  return elementLocalId(bulkData_->identifier(elmt));
1422  // const std::size_t * fieldCoords = stk::mesh::field_data(*localIdField_,*elmt);
1423  // return fieldCoords[0];
1424 }
1425 
1426 std::size_t STK_Interface::elementLocalId(stk::mesh::EntityId gid) const
1427 {
1428  // stk::mesh::EntityRank elementRank = getElementRank();
1429  // stk::mesh::Entity elmt = bulkData_->get_entity(elementRank,gid);
1430  // TEUCHOS_ASSERT(elmt->owner_rank()==procRank_);
1431  // return elementLocalId(elmt);
1432  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localIDHash_.find(gid);
1433  TEUCHOS_ASSERT(itr!=localIDHash_.end());
1434  return itr->second;
1435 }
1436 
1437 bool STK_Interface::isEdgeLocal(stk::mesh::Entity edge) const
1438 {
1439  return isEdgeLocal(bulkData_->identifier(edge));
1440 }
1441 
1442 bool STK_Interface::isEdgeLocal(stk::mesh::EntityId gid) const
1443 {
1444  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localEdgeIDHash_.find(gid);
1445  if (itr==localEdgeIDHash_.end()) {
1446  return false;
1447  }
1448  return true;
1449 }
1450 
1451 std::size_t STK_Interface::edgeLocalId(stk::mesh::Entity edge) const
1452 {
1453  return edgeLocalId(bulkData_->identifier(edge));
1454 }
1455 
1456 std::size_t STK_Interface::edgeLocalId(stk::mesh::EntityId gid) const
1457 {
1458  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localEdgeIDHash_.find(gid);
1459  TEUCHOS_ASSERT(itr!=localEdgeIDHash_.end());
1460  return itr->second;
1461 }
1462 
1463 bool STK_Interface::isFaceLocal(stk::mesh::Entity face) const
1464 {
1465  return isFaceLocal(bulkData_->identifier(face));
1466 }
1467 
1468 bool STK_Interface::isFaceLocal(stk::mesh::EntityId gid) const
1469 {
1470  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localFaceIDHash_.find(gid);
1471  if (itr==localFaceIDHash_.end()) {
1472  return false;
1473  }
1474  return true;
1475 }
1476 
1477 std::size_t STK_Interface::faceLocalId(stk::mesh::Entity face) const
1478 {
1479  return faceLocalId(bulkData_->identifier(face));
1480 }
1481 
1482 std::size_t STK_Interface::faceLocalId(stk::mesh::EntityId gid) const
1483 {
1484  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localFaceIDHash_.find(gid);
1485  TEUCHOS_ASSERT(itr!=localFaceIDHash_.end());
1486  return itr->second;
1487 }
1488 
1489 
1490 std::string STK_Interface::containingBlockId(stk::mesh::Entity elmt) const
1491 {
1492  for(const auto & eb_pair : elementBlocks_)
1493  if(bulkData_->bucket(elmt).member(*(eb_pair.second)))
1494  return eb_pair.first;
1495  return "";
1496 }
1497 
1498 stk::mesh::Field<double> * STK_Interface::getSolutionField(const std::string & fieldName,
1499  const std::string & blockId) const
1500 {
1501  // look up field in map
1502  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1503  iter = fieldNameToSolution_.find(std::make_pair(fieldName,blockId));
1504 
1505  // check to make sure field was actually found
1506  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToSolution_.end(),std::runtime_error,
1507  "Solution field name \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1508 
1509  return iter->second;
1510 }
1511 
1512 stk::mesh::Field<double> * STK_Interface::getCellField(const std::string & fieldName,
1513  const std::string & blockId) const
1514 {
1515  // look up field in map
1516  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1517  iter = fieldNameToCellField_.find(std::make_pair(fieldName,blockId));
1518 
1519  // check to make sure field was actually found
1520  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToCellField_.end(),std::runtime_error,
1521  "Cell field named \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1522 
1523  return iter->second;
1524 }
1525 
1526 stk::mesh::Field<double> * STK_Interface::getEdgeField(const std::string & fieldName,
1527  const std::string & blockId) const
1528 {
1529  // look up field in map
1530  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1531  iter = fieldNameToEdgeField_.find(std::make_pair(fieldName,blockId));
1532 
1533  // check to make sure field was actually found
1534  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToEdgeField_.end(),std::runtime_error,
1535  "Edge field named \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1536 
1537  return iter->second;
1538 }
1539 
1540 stk::mesh::Field<double> * STK_Interface::getFaceField(const std::string & fieldName,
1541  const std::string & blockId) const
1542 {
1543  // look up field in map
1544  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1545  iter = fieldNameToFaceField_.find(std::make_pair(fieldName,blockId));
1546 
1547  // check to make sure field was actually found
1548  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToFaceField_.end(),std::runtime_error,
1549  "Face field named \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1550 
1551  return iter->second;
1552 }
1553 
1555 {
1556  using Teuchos::RCP;
1557  using Teuchos::rcp;
1558 
1559  if(orderedElementVector_==Teuchos::null) {
1560  // safe because essentially this is a call to modify a mutable object
1561  const_cast<STK_Interface*>(this)->buildLocalElementIDs();
1562  }
1563 
1565 }
1566 
1567 void STK_Interface::addElementBlock(const std::string & name,const CellTopologyData * ctData)
1568 {
1570 
1571  stk::mesh::Part * block = metaData_->get_part(name);
1572  if(block==0) {
1573  block = &metaData_->declare_part_with_topology(name, stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
1574  }
1575 
1576  // construct cell topology object for this block
1578  = Teuchos::rcp(new shards::CellTopology(ctData));
1579 
1580  // add element block part and cell topology
1581  elementBlocks_.insert(std::make_pair(name,block));
1582  elementBlockCT_.insert(std::make_pair(name,ct));
1583 }
1584 
1586 {
1587  using Teuchos::RCP;
1588  using Teuchos::rcp;
1589 
1590  if(orderedEdgeVector_==Teuchos::null) {
1591  // safe because essentially this is a call to modify a mutable object
1592  const_cast<STK_Interface*>(this)->buildLocalEdgeIDs();
1593  }
1594 
1595  return orderedEdgeVector_.getConst();
1596 }
1597 
1598 void STK_Interface::addEdgeBlock(const std::string & elemBlockName,
1599  const std::string & edgeBlockName,
1600  const stk::topology & topology)
1601 {
1603 
1604  stk::mesh::Part * edge_block = metaData_->get_part(edgeBlockName);
1605  if(edge_block==0) {
1606  edge_block = &metaData_->declare_part_with_topology(edgeBlockName, topology);
1607  }
1608 
1609  /* There is only one edge block for each edge topology, so declare it
1610  * as a subset of the element block even if it wasn't just created.
1611  */
1612  stk::mesh::Part * elem_block = metaData_->get_part(elemBlockName);
1613  metaData_->declare_part_subset(*elem_block, *edge_block);
1614 
1615  // add edge block part
1616  edgeBlocks_.insert(std::make_pair(edgeBlockName,edge_block));
1617 }
1618 
1619 void STK_Interface::addEdgeBlock(const std::string & elemBlockName,
1620  const std::string & edgeBlockName,
1621  const CellTopologyData * ctData)
1622 {
1624 
1625  addEdgeBlock(elemBlockName,
1626  edgeBlockName,
1627  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
1628 }
1629 
1631 {
1632  using Teuchos::RCP;
1633  using Teuchos::rcp;
1634 
1635  if(orderedFaceVector_==Teuchos::null) {
1636  // safe because essentially this is a call to modify a mutable object
1637  const_cast<STK_Interface*>(this)->buildLocalFaceIDs();
1638  }
1639 
1640  return orderedFaceVector_.getConst();
1641 }
1642 
1643 void STK_Interface::addFaceBlock(const std::string & elemBlockName,
1644  const std::string & faceBlockName,
1645  const stk::topology & topology)
1646 {
1648 
1649  stk::mesh::Part * face_block = metaData_->get_part(faceBlockName);
1650  if(face_block==0) {
1651  face_block = &metaData_->declare_part_with_topology(faceBlockName, topology);
1652  }
1653 
1654  /* There is only one face block for each edge topology, so declare it
1655  * as a subset of the element block even if it wasn't just created.
1656  */
1657  stk::mesh::Part * elem_block = metaData_->get_part(elemBlockName);
1658  metaData_->declare_part_subset(*elem_block, *face_block);
1659 
1660  // add face block part
1661  faceBlocks_.insert(std::make_pair(faceBlockName,face_block));
1662 }
1663 
1664 void STK_Interface::addFaceBlock(const std::string & elemBlockName,
1665  const std::string & faceBlockName,
1666  const CellTopologyData * ctData)
1667 {
1669 
1670  addFaceBlock(elemBlockName,
1671  faceBlockName,
1672  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
1673 }
1674 
1676 {
1677  dimension_ = metaData_->spatial_dimension();
1678 
1679  // declare coordinates and node parts
1680  coordinatesField_ = &metaData_->declare_field<double>(stk::topology::NODE_RANK, coordsString);
1681  edgesField_ = &metaData_->declare_field<double>(stk::topology::EDGE_RANK, edgesString);
1682  if (dimension_ > 2)
1683  facesField_ = &metaData_->declare_field<double>(stk::topology::FACE_RANK, facesString);
1684  processorIdField_ = &metaData_->declare_field<ProcIdData>(stk::topology::ELEMENT_RANK, "PROC_ID");
1685  loadBalField_ = &metaData_->declare_field<double>(stk::topology::ELEMENT_RANK, "LOAD_BAL");
1686 
1687  // stk::mesh::put_field( *coordinatesField_ , getNodeRank() , metaData_->universal_part() );
1688 
1689  nodesPart_ = &metaData_->declare_part(nodesString,getNodeRank());
1690  nodesPartVec_.push_back(nodesPart_);
1691  edgesPart_ = &metaData_->declare_part(edgesString,getEdgeRank());
1692  edgesPartVec_.push_back(edgesPart_);
1693  if (dimension_ > 2) {
1694  facesPart_ = &metaData_->declare_part(facesString,getFaceRank());
1695  facesPartVec_.push_back(facesPart_);
1696  }
1697 }
1698 
1700 {
1701  currentLocalId_ = 0;
1702 
1703  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
1704 
1705  // might be better (faster) to do this by buckets
1706  std::vector<stk::mesh::Entity> elements;
1707  getMyElements(elements);
1708 
1709  for(std::size_t index=0;index<elements.size();++index) {
1710  stk::mesh::Entity element = elements[index];
1711 
1712  // set processor rank
1713  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1714  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1715 
1716  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1717 
1718  currentLocalId_++;
1719  }
1720 
1721  // copy elements into the ordered element vector
1722  orderedElementVector_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>(elements));
1723 
1724  elements.clear();
1725  getNeighborElements(elements);
1726 
1727  for(std::size_t index=0;index<elements.size();++index) {
1728  stk::mesh::Entity element = elements[index];
1729 
1730  // set processor rank
1731  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1732  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1733 
1734  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1735 
1736  currentLocalId_++;
1737  }
1738 
1739  orderedElementVector_->insert(orderedElementVector_->end(),elements.begin(),elements.end());
1740 }
1741 
1743 {
1744  std::vector<std::string> names;
1745  getElementBlockNames(names);
1746 
1747  for(std::size_t b=0;b<names.size();b++) {
1748  // find user specified block weight, otherwise use 1.0
1749  std::map<std::string,double>::const_iterator bw_itr = blockWeights_.find(names[b]);
1750  double blockWeight = (bw_itr!=blockWeights_.end()) ? bw_itr->second : 1.0;
1751 
1752  std::vector<stk::mesh::Entity> elements;
1753  getMyElements(names[b],elements);
1754 
1755  for(std::size_t index=0;index<elements.size();++index) {
1756  // set local element ID
1757  double * loadBal = stk::mesh::field_data(*loadBalField_,elements[index]);
1758  loadBal[0] = blockWeight;
1759  }
1760  }
1761 }
1762 
1764 {
1765  currentLocalId_ = 0;
1766 
1767  orderedEdgeVector_ = Teuchos::null; // forces rebuild of ordered lists
1768 
1769  // might be better (faster) to do this by buckets
1770  std::vector<stk::mesh::Entity> edges;
1771  getMyEdges(edges);
1772 
1773  for(std::size_t index=0;index<edges.size();++index) {
1774  stk::mesh::Entity edge = edges[index];
1775  localEdgeIDHash_[bulkData_->identifier(edge)] = currentLocalId_;
1776  currentLocalId_++;
1777  }
1778 
1779  // copy edges into the ordered edge vector
1780  orderedEdgeVector_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>(edges));
1781 }
1782 
1784 {
1785  currentLocalId_ = 0;
1786 
1787  orderedFaceVector_ = Teuchos::null; // forces rebuild of ordered lists
1788 
1789  // might be better (faster) to do this by buckets
1790  std::vector<stk::mesh::Entity> faces;
1791  getMyFaces(faces);
1792 
1793  for(std::size_t index=0;index<faces.size();++index) {
1794  stk::mesh::Entity face = faces[index];
1795  localFaceIDHash_[bulkData_->identifier(face)] = currentLocalId_;
1796  currentLocalId_++;
1797  }
1798 
1799  // copy faces into the ordered face vector
1800  orderedFaceVector_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>(faces));
1801 }
1802 
1803 bool
1804 STK_Interface::isMeshCoordField(const std::string & eBlock,
1805  const std::string & fieldName,
1806  int & axis) const
1807 {
1808  std::map<std::string,std::vector<std::string> >::const_iterator blkItr = meshCoordFields_.find(eBlock);
1809  if(blkItr==meshCoordFields_.end()) {
1810  return false;
1811  }
1812 
1813  axis = 0;
1814  for(axis=0;axis<Teuchos::as<int>(blkItr->second.size());axis++) {
1815  if(blkItr->second[axis]==fieldName)
1816  break; // found axis, break
1817  }
1818 
1819  if(axis>=Teuchos::as<int>(blkItr->second.size()))
1820  return false;
1821 
1822  return true;
1823 }
1824 
1825 std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
1827 {
1829  Teuchos::RCP<std::vector<unsigned int > > type_vec = rcp(new std::vector<unsigned int>);
1830  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & matchers = getPeriodicBCVector();
1831  const bool & useBBoxSearch = useBoundingBoxSearch();
1832  std::vector<std::vector<std::string> > matchedSides(3); // (coord,edge,face)
1833 
1834  // build up the vectors by looping over the matched pair
1835  for(std::size_t m=0;m<matchers.size();m++){
1836  unsigned int type;
1837  if(matchers[m]->getType() == "coord")
1838  type = 0;
1839  else if(matchers[m]->getType() == "edge")
1840  type = 1;
1841  else if(matchers[m]->getType() == "face")
1842  type = 2;
1843  else
1844  TEUCHOS_ASSERT(false);
1845 #ifdef PANZER_HAVE_STKSEARCH
1846 
1847  if (useBBoxSearch) {
1848  vec = matchers[m]->getMatchedPair(*this,matchedSides[type],vec);
1849  } else {
1850  vec = matchers[m]->getMatchedPair(*this,vec);
1851  }
1852 #else
1853  TEUCHOS_TEST_FOR_EXCEPTION(useBBoxSearch,std::logic_error,
1854  "panzer::STK_Interface::getPeriodicNodePairing(): Requested bounding box search, but "
1855  "did not compile with STK_SEARCH enabled.");
1856  vec = matchers[m]->getMatchedPair(*this,vec);
1857 #endif
1858  type_vec->insert(type_vec->begin(),vec->size()-type_vec->size(),type);
1859  matchedSides[type].push_back(matchers[m]->getLeftSidesetName());
1860  }
1861 
1862  return std::make_pair(vec,type_vec);
1863 }
1864 
1865 bool STK_Interface::validBlockId(const std::string & blockId) const
1866 {
1867  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr = elementBlocks_.find(blockId);
1868 
1869  return blkItr!=elementBlocks_.end();
1870 }
1871 
1872 void STK_Interface::print(std::ostream & os) const
1873 {
1874  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1875 
1876  getElementBlockNames(blockNames);
1877  getSidesetNames(sidesetNames);
1878  getNodesetNames(nodesetNames);
1879 
1880  os << "STK Mesh data:\n";
1881  os << " Spatial dim = " << getDimension() << "\n";
1882  if(getDimension()==2)
1883  os << " Entity counts (Nodes, Edges, Cells) = ( "
1884  << getEntityCounts(getNodeRank()) << ", "
1885  << getEntityCounts(getEdgeRank()) << ", "
1886  << getEntityCounts(getElementRank()) << " )\n";
1887  else if(getDimension()==3)
1888  os << " Entity counts (Nodes, Edges, Faces, Cells) = ( "
1889  << getEntityCounts(getNodeRank()) << ", "
1890  << getEntityCounts(getEdgeRank()) << ", "
1891  << getEntityCounts(getSideRank()) << ", "
1892  << getEntityCounts(getElementRank()) << " )\n";
1893  else
1894  os << " Entity counts (Nodes, Cells) = ( "
1895  << getEntityCounts(getNodeRank()) << ", "
1896  << getEntityCounts(getElementRank()) << " )\n";
1897 
1898  os << " Element blocks = ";
1899  for(std::size_t i=0;i<blockNames.size();i++)
1900  os << "\"" << blockNames[i] << "\" ";
1901  os << "\n";
1902  os << " Sidesets = ";
1903  for(std::size_t i=0;i<sidesetNames.size();i++)
1904  os << "\"" << sidesetNames[i] << "\" ";
1905  os << "\n";
1906  os << " Nodesets = ";
1907  for(std::size_t i=0;i<nodesetNames.size();i++)
1908  os << "\"" << nodesetNames[i] << "\" ";
1909  os << std::endl;
1910 
1911  // print out periodic boundary conditions
1912  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1913  = getPeriodicBCVector();
1914  if(bcVector.size()!=0) {
1915  os << " Periodic BCs:\n";
1916  for(std::size_t i=0;i<bcVector.size();i++)
1917  os << " " << bcVector[i]->getString() << "\n";
1918  os << std::endl;
1919  }
1920 }
1921 
1922 void STK_Interface::printMetaData(std::ostream & os) const
1923 {
1924  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1925 
1926  getElementBlockNames(blockNames);
1927  getSidesetNames(sidesetNames);
1928  getNodesetNames(nodesetNames);
1929 
1930  os << "STK Meta data:\n";
1931  os << " Element blocks = ";
1932  for(std::size_t i=0;i<blockNames.size();i++)
1933  os << "\"" << blockNames[i] << "\" ";
1934  os << "\n";
1935  os << " Sidesets = ";
1936  for(std::size_t i=0;i<sidesetNames.size();i++)
1937  os << "\"" << sidesetNames[i] << "\" ";
1938  os << "\n";
1939  os << " Nodesets = ";
1940  for(std::size_t i=0;i<nodesetNames.size();i++)
1941  os << "\"" << nodesetNames[i] << "\" ";
1942  os << std::endl;
1943 
1944  // print out periodic boundary conditions
1945  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1946  = getPeriodicBCVector();
1947  if(bcVector.size()!=0) {
1948  os << " Periodic BCs:\n";
1949  for(std::size_t i=0;i<bcVector.size();i++)
1950  os << " " << bcVector[i]->getString() << "\n";
1951  os << std::endl;
1952  }
1953 
1954  // print all available fields in meta data
1955  os << " Fields = ";
1956  const stk::mesh::FieldVector & fv = metaData_->get_fields();
1957  for(std::size_t i=0;i<fv.size();i++)
1958  os << "\"" << fv[i]->name() << "\" ";
1959  os << std::endl;
1960 }
1961 
1963 {
1964  std::map<std::string, Teuchos::RCP<shards::CellTopology> >::const_iterator itr;
1965  itr = elementBlockCT_.find(eBlock);
1966 
1967  if(itr==elementBlockCT_.end()) {
1968  std::stringstream ss;
1969  printMetaData(ss);
1970  TEUCHOS_TEST_FOR_EXCEPTION(itr==elementBlockCT_.end(),std::logic_error,
1971  "STK_Interface::getCellTopology: No such element block \"" +eBlock +"\" available.\n\n"
1972  << "STK Meta Data follows: \n" << ss.str());
1973  }
1974 
1975  return itr->second;
1976 }
1977 
1979 {
1980  MPI_Comm newComm;
1981  const int err = MPI_Comm_dup (parallelMach, &newComm);
1982  TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error,
1983  "panzer::STK_Interface: MPI_Comm_dup failed with error \""
1984  << Teuchos::mpiErrorCodeToString (err) << "\".");
1985 
1986  return Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper (newComm,MPI_Comm_free)));
1987 }
1988 
1990 {
1991  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Rebalance not currently supported");
1992 #if 0
1993  // make sure weights have been set (a local operation)
1995 
1996  stk::mesh::Selector selector(getMetaData()->universal_part());
1997  stk::mesh::Selector owned_selector(getMetaData()->locally_owned_part());
1998 
1999  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
2000  out.setOutputToRootOnly(0);
2001  out << "Load balance before: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
2002 
2003  // perform reblance
2004  Teuchos::ParameterList graph;
2005  if(params.begin()!=params.end())
2006  graph.sublist(stk::rebalance::Zoltan::default_parameters_name()) = params;
2007  stk::rebalance::Zoltan zoltan_partition(*mpiComm_->getRawMpiComm(), getDimension(), graph);
2008  stk::rebalance::rebalance(*getBulkData(), owned_selector, &getCoordinatesField(), loadBalField_, zoltan_partition);
2009 
2010  out << "Load balance after: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
2011 
2012  currentLocalId_ = 0;
2013  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
2014 #endif
2015 }
2016 
2019 {
2020  TEUCHOS_ASSERT(this->isInitialized());
2021  return mpiComm_;
2022 }
2023 
2024 void STK_Interface::refineMesh(const int numberOfLevels, const bool deleteParentElements) {
2025 #ifdef PANZER_HAVE_PERCEPT
2026  TEUCHOS_TEST_FOR_EXCEPTION(numberOfLevels < 1,std::runtime_error,
2027  "ERROR: Number of levels for uniform refinement must be greater than 0");
2028  TEUCHOS_ASSERT(nonnull(refinedMesh_));
2029  TEUCHOS_ASSERT(nonnull(breakPattern_));
2030 
2031  refinedMesh_->setCoordinatesField();
2032 
2033  percept::UniformRefiner breaker(*refinedMesh_,*breakPattern_);
2034  breaker.setRemoveOldElements(deleteParentElements);
2035 
2036  for (int i=0; i < numberOfLevels; ++i)
2037  breaker.doBreak();
2038 
2039 #else
2040  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
2041  "ERROR: Uniform refinement requested. This requires the Percept package to be enabled in Trilinos!");
2042 #endif
2043 }
2044 
2045 
2046 } // end namespace panzer_stk
Teuchos::RCP< ElementDescriptor > buildElementDescriptor(stk::mesh::EntityId elmtId, std::vector< stk::mesh::EntityId > &nodes)
Teuchos::RCP< Teuchos::MpiComm< int > > getSafeCommunicator(stk::ParallelMachine parallelMach) const
void getSidesetNames(std::vector< std::string > &name) const
RCP< const T > getConst() const
std::pair< Teuchos::RCP< std::vector< std::pair< std::size_t, std::size_t > > >, Teuchos::RCP< std::vector< unsigned int > > > getPeriodicNodePairing() const
const bool & useBoundingBoxSearch() const
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedEdgeVector_
void addNodeset(const std::string &name)
Teuchos::RCP< stk::mesh::BulkData > bulkData_
void getFaceBlockNames(std::vector< std::string > &names) const
void print(std::ostream &os) const
void getMyEdges(std::vector< stk::mesh::Entity > &edges) const
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getElementsOrderedByLID() const
std::vector< stk::mesh::Part * > facesPartVec_
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToSolution_
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void addGlobalToExodus(const std::string &key, const int &value)
Add an int global variable to the information to be written to the Exodus output file.
void getElementBlockNames(std::vector< std::string > &names) const
void getSubcellIndices(unsigned entityRank, stk::mesh::EntityId elementId, std::vector< stk::mesh::EntityId > &subcellIds) const
void getElementsSharingNodes(const std::vector< stk::mesh::EntityId > nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing multiple nodes
stk::mesh::Field< double > SolutionFieldType
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
std::vector< stk::mesh::Part * > edgesPartVec_
bool isFaceLocal(stk::mesh::Entity face) const
std::map< std::string, double > blockWeights_
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getEdgesOrderedByLID() const
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::size_t getEntityCounts(unsigned entityRank) const
get the global counts for the entity of specified rank
void printMetaData(std::ostream &os) const
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
std::map< std::string, stk::mesh::Part * > faceBlocks_
void addSolutionField(const std::string &fieldName, const std::string &blockId)
std::map< std::string, Teuchos::RCP< shards::CellTopology > > elementBlockCT_
void addEntityToFaceBlock(stk::mesh::Entity entity, stk::mesh::Part *faceblock)
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToFaceField_
std::map< std::string, stk::mesh::Part * > sidesets_
std::vector< std::size_t > entityCounts_
bool validBlockId(const std::string &blockId) const
stk::mesh::Part * getNodeset(const std::string &name) const
void getEdgeBlockNames(std::vector< std::string > &names) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
void addInformationRecords(const std::vector< std::string > &info_records)
void endModification(const bool find_and_set_shared_nodes_in_stk=true)
void getAllFaces(const std::string &faceBlockName, std::vector< stk::mesh::Entity > &faces) const
stk::mesh::Field< double > * getSolutionField(const std::string &fieldName, const std::string &blockId) const
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const
stk::mesh::Part * getSideset(const std::string &name) const
static const std::string nodesString
void instantiateBulkData(stk::ParallelMachine parallelMach)
std::set< std::string > informationRecords_
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
stk::mesh::Field< double > * getEdgeField(const std::string &fieldName, const std::string &blockId) const
std::size_t elementLocalId(stk::mesh::Entity elmt) const
unsigned getDimension() const
get the dimension
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
std::map< std::string, std::vector< std::string > > meshDispFields_
void setupExodusFile(const std::string &filename, const bool append=false, const bool append_after_restart_time=false, const double restart_time=0.0)
Set up an output Exodus file for writing results.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
stk::mesh::EntityRank getSideRank() const
void getAllEdges(const std::string &edgeBlockName, std::vector< stk::mesh::Entity > &edges) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
static const std::string coordsString
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
std::map< std::string, stk::mesh::Part * > edgeBlocks_
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
std::map< std::string, stk::mesh::Part * > elementBlocks_
std::vector< stk::mesh::Part * > nodesPartVec_
static const std::string edgeBlockString
static const std::string edgesString
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
std::string containingBlockId(stk::mesh::Entity elmt) const
stk::mesh::EntityRank getFaceRank() const
std::size_t edgeLocalId(stk::mesh::Entity elmt) const
void getOwnedElementsSharingNode(stk::mesh::Entity node, std::vector< stk::mesh::Entity > &elements, std::vector< int > &relIds) const
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void buildSubcells()
force the mesh to build subcells: edges and faces
Teuchos::RCP< Teuchos::MpiComm< int > > mpiComm_
Teuchos::RCP< stk::mesh::MetaData > metaData_
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedElementVector_
bool isInitialized() const
Has initialize been called on this mesh object?
TransListIter iter
void writeToExodus(const std::string &filename, const bool append=false)
Write this mesh and associated fields to the given output file.
void getMyNodes(const std::string &sideName, const std::string &blockName, std::vector< stk::mesh::Entity > &nodes) const
bool nonnull(const boost::shared_ptr< T > &p)
void addEdgeField(const std::string &fieldName, const std::string &blockId)
void addMeshCoordFields(const std::string &blockId, const std::vector< std::string > &coordField, const std::string &dispPrefix)
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
bool isMeshCoordField(const std::string &eBlock, const std::string &fieldName, int &axis) const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
static const std::string faceBlockString
std::unordered_map< stk::mesh::EntityId, std::size_t > localFaceIDHash_
void addEntityToEdgeBlock(stk::mesh::Entity entity, stk::mesh::Part *edgeblock)
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToEdgeField_
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::unordered_map< stk::mesh::EntityId, std::size_t > localEdgeIDHash_
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
void rebalance(const Teuchos::ParameterList &params)
std::map< std::string, std::vector< std::string > > meshCoordFields_
stk::mesh::EntityRank getEdgeRank() const
static const std::string facesString
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
void initializeFieldsInSTK(const std::map< std::pair< std::string, std::string >, SolutionFieldType * > &nameToField, bool setupIO)
const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector() const
void refineMesh(const int numberOfLevels, const bool deleteParentElements)
stk::mesh::Field< double > * getFaceField(const std::string &fieldName, const std::string &blockId) const
std::size_t faceLocalId(stk::mesh::Entity elmt) const
void addFaceBlock(const std::string &elemBlockName, const std::string &faceBlockName, const stk::topology &topology)
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedFaceVector_
std::unordered_map< stk::mesh::EntityId, std::size_t > localIDHash_
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
std::vector< stk::mesh::EntityId > maxEntityId_
void getElementsSharingNode(stk::mesh::EntityId nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing a single node
void addFaceField(const std::string &fieldName, const std::string &blockId)
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
#define TEUCHOS_ASSERT(assertion_test)
stk::mesh::EntityRank getNodeRank() const
void getMyFaces(std::vector< stk::mesh::Entity > &faces) const
void addCellField(const std::string &fieldName, const std::string &blockId)
void getNodeIdsForElement(stk::mesh::Entity element, std::vector< stk::mesh::EntityId > &nodeIds) const
get a list of node ids for nodes connected to an element
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getFacesOrderedByLID() const
stk::mesh::EntityRank getElementRank() const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
const VectorFieldType & getCoordinatesField() const
bool isEdgeLocal(stk::mesh::Entity edge) const
void getNeighborElements(std::vector< stk::mesh::Entity > &elements) const
void getNodesetNames(std::vector< std::string > &name) const
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToCellField_
std::map< std::string, stk::mesh::Part * > nodesets_
bool is_null() const