Panzer  Version of the Day
Panzer_STK_Interface.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 #include <PanzerAdaptersSTK_config.hpp>
44 #include <Panzer_STK_Interface.hpp>
45 
46 #include <Teuchos_as.hpp>
47 
48 #include <limits>
49 
50 #include <stk_mesh/base/FieldBase.hpp>
51 #include <stk_mesh/base/Comm.hpp>
52 #include <stk_mesh/base/Selector.hpp>
53 #include <stk_mesh/base/GetEntities.hpp>
54 #include <stk_mesh/base/GetBuckets.hpp>
55 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
56 
57 // #include <stk_rebalance/Rebalance.hpp>
58 // #include <stk_rebalance/Partition.hpp>
59 // #include <stk_rebalance/ZoltanPartition.hpp>
60 // #include <stk_rebalance_utils/RebalanceUtils.hpp>
61 
62 #include <stk_util/parallel/ParallelReduce.hpp>
63 #include <stk_util/parallel/CommSparse.hpp>
64 
65 #ifdef PANZER_HAVE_IOSS
66 #include <Ionit_Initializer.h>
67 #include <stk_io/IossBridge.hpp>
68 #endif
69 
71 
72 #include <set>
73 
74 using Teuchos::RCP;
75 using Teuchos::rcp;
76 
77 namespace panzer_stk {
78 
80 ElementDescriptor::ElementDescriptor(stk::mesh::EntityId gid,const std::vector<stk::mesh::EntityId> & nodes)
81  : gid_(gid), nodes_(nodes) {}
83 
87 buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes)
88 {
89  return Teuchos::rcp(new ElementDescriptor(elmtId,nodes));
90 }
91 
92 const std::string STK_Interface::coordsString = "coordinates";
93 const std::string STK_Interface::nodesString = "nodes";
94 const std::string STK_Interface::edgesString = "edges";
95 const std::string STK_Interface::facesString = "faces";
96 
98  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
99 {
100  metaData_ = rcp(new stk::mesh::MetaData());
101 }
102 
104  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
105 {
106  metaData_ = metaData;
107 }
108 
110  : dimension_(dim), initialized_(false), currentLocalId_(0), useFieldCoordinates_(false)
111 {
112  std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
113  entity_rank_names.push_back("FAMILY_TREE");
114 
115  metaData_ = rcp(new stk::mesh::MetaData(dimension_,entity_rank_names));
116 
118 }
119 
120 void STK_Interface::addSideset(const std::string & name,const CellTopologyData * ctData)
121 {
124 
125  stk::mesh::Part * sideset = metaData_->get_part(name);
126  if(sideset==NULL)
127  sideset = &metaData_->declare_part_with_topology(name,
128  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
129  sidesets_.insert(std::make_pair(name,sideset));
130 }
131 
132 void STK_Interface::addNodeset(const std::string & name)
133 {
136 
137  stk::mesh::Part * nodeset = metaData_->get_part(name);
138  if(nodeset==NULL) {
139  const CellTopologyData * ctData = shards::getCellTopologyData<shards::Node>();
140  nodeset = &metaData_->declare_part_with_topology(name,
141  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
142  }
143  nodesets_.insert(std::make_pair(name,nodeset));
144 }
145 
146 void STK_Interface::addSolutionField(const std::string & fieldName,const std::string & blockId)
147 {
150  "Unknown element block \"" << blockId << "\"");
151  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
152 
153  // add & declare field if not already added...currently assuming linears
154  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
155  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::NODE_RANK, fieldName);
156  if(field==0)
157  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::NODE_RANK, fieldName);
159  }
160 }
161 
162 void STK_Interface::addCellField(const std::string & fieldName,const std::string & blockId)
163 {
166  "Unknown element block \"" << blockId << "\"");
167  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
168 
169  // add & declare field if not already added...currently assuming linears
170  if(fieldNameToCellField_.find(key)==fieldNameToCellField_.end()) {
171  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, fieldName);
172  if(field==0)
173  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, fieldName);
175  }
176 }
177 
178 void STK_Interface::addMeshCoordFields(const std::string & blockId,
179  const std::vector<std::string> & coordNames,
180  const std::string & dispPrefix)
181 {
183  TEUCHOS_ASSERT(dimension_==coordNames.size());
186  "Unknown element block \"" << blockId << "\"");
187 
188  // we only allow one alternative coordinate field
189  TEUCHOS_TEST_FOR_EXCEPTION(meshCoordFields_.find(blockId)!=meshCoordFields_.end(),std::invalid_argument,
190  "STK_Interface::addMeshCoordFields: Can't set more than one set of coordinate "
191  "fields for element block \""+blockId+"\".");
192 
193  // Note that there is a distinction between the key which is used for lookups
194  // and the field that lives on the mesh, which is used for printing the displacement.
195 
196  // just copy the coordinate names
197  meshCoordFields_[blockId] = coordNames;
198 
199  // must fill in the displacement fields
200  std::vector<std::string> & dispFields = meshDispFields_[blockId];
201  dispFields.resize(dimension_);
202 
203  for(unsigned i=0;i<dimension_;i++) {
204  std::pair<std::string,std::string> key = std::make_pair(coordNames[i],blockId);
205  std::string dispName = dispPrefix+coordNames[i];
206 
207  dispFields[i] = dispName; // record this field as a
208  // displacement field
209 
210  // add & declare field if not already added...currently assuming linears
211  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
212 
213  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::NODE_RANK, dispName);
214  if(field==0) {
215  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::NODE_RANK, dispName);
216  }
218  }
219  }
220 }
221 
222 void STK_Interface::initialize(stk::ParallelMachine parallelMach,bool setupIO)
223 {
225  TEUCHOS_ASSERT(dimension_!=0); // no zero dimensional meshes!
226 
227  if(mpiComm_==Teuchos::null)
228  mpiComm_ = getSafeCommunicator(parallelMach);
229 
230  procRank_ = stk::parallel_machine_rank(*mpiComm_->getRawMpiComm());
231 
232  // associating the field with a part: universal part!
233  stk::mesh::put_field( *coordinatesField_ , metaData_->universal_part(), getDimension());
234  stk::mesh::put_field( *edgesField_ , metaData_->universal_part(), getDimension());
235  if (dimension_ > 2)
236  stk::mesh::put_field( *facesField_ , metaData_->universal_part(), getDimension());
237  stk::mesh::put_field( *processorIdField_ , metaData_->universal_part());
238  stk::mesh::put_field( *loadBalField_ , metaData_->universal_part());
239 
242 
243 #ifdef PANZER_HAVE_IOSS
244  if(setupIO) {
245  // setup Exodus file IO
247 
248  // add element blocks
249  {
250  std::map<std::string, stk::mesh::Part*>::iterator itr;
251  for(itr=elementBlocks_.begin();itr!=elementBlocks_.end();++itr)
252  if(!stk::io::is_part_io_part(*itr->second))
253  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
254  }
255 
256  // add side sets
257  {
258  std::map<std::string, stk::mesh::Part*>::iterator itr;
259  for(itr=sidesets_.begin();itr!=sidesets_.end();++itr)
260  if(!stk::io::is_part_io_part(*itr->second))
261  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
262  }
263 
264  // add node sets
265  {
266  std::map<std::string, stk::mesh::Part*>::iterator itr;
267  for(itr=nodesets_.begin();itr!=nodesets_.end();++itr)
268  if(!stk::io::is_part_io_part(*itr->second))
269  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
270  }
271 
272  // add nodes
273  if(!stk::io::is_part_io_part(*nodesPart_))
274  stk::io::put_io_part_attribute(*nodesPart_);
275 
276  stk::io::set_field_role(*coordinatesField_, Ioss::Field::MESH);
277  stk::io::set_field_role(*edgesField_, Ioss::Field::MESH);
278  if (dimension_ > 2)
279  stk::io::set_field_role(*facesField_, Ioss::Field::MESH);
280  stk::io::set_field_role(*processorIdField_, Ioss::Field::TRANSIENT);
281  // stk::io::set_field_role(*loadBalField_, Ioss::Field::TRANSIENT);
282  }
283 #endif
284 
285  metaData_->commit();
286  if(bulkData_==Teuchos::null)
287  instantiateBulkData(*mpiComm_->getRawMpiComm());
288 
289  initialized_ = true;
290 }
291 
292 void STK_Interface::initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
293  bool setupIO)
294 {
295  std::set<SolutionFieldType*> uniqueFields;
296  std::map<std::pair<std::string,std::string>,SolutionFieldType*>::const_iterator fieldIter;
297  for(fieldIter=nameToField.begin();fieldIter!=nameToField.end();++fieldIter)
298  uniqueFields.insert(fieldIter->second); // this makes setting up IO easier!
299 
300  {
301  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
302  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
303  stk::mesh::put_field(*(*uniqueFieldIter),metaData_->universal_part());
304  }
305 
306 #ifdef PANZER_HAVE_IOSS
307  if(setupIO) {
308  // add solution fields
309  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
310  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
311  stk::io::set_field_role(*(*uniqueFieldIter), Ioss::Field::TRANSIENT);
312  }
313 #endif
314 }
315 
316 void STK_Interface::instantiateBulkData(stk::ParallelMachine parallelMach)
317 {
318  TEUCHOS_ASSERT(bulkData_==Teuchos::null);
319  if(mpiComm_==Teuchos::null)
320  mpiComm_ = getSafeCommunicator(parallelMach);
321 
322  bulkData_ = rcp(new stk::mesh::BulkData(*metaData_, *mpiComm_->getRawMpiComm()));
323 }
324 
326 {
327  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
328  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"beginModification\"");
329 
330  bulkData_->modification_begin();
331 }
332 
334 {
335  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
336  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"endModification\"");
337 
338  // TODO: Resolving sharing this way comes at a cost in performance. The STK
339  // team has decided that users need to declare their own sharing. We should
340  // find where shared entities are being created in Panzer and declare it.
341  // Once this is done, the extra code below can be deleted.
342 
343  stk::CommSparse comm(bulkData_->parallel());
344 
345  for (int phase=0;phase<2;++phase) {
346  for (int i=0;i<bulkData_->parallel_size();++i) {
347  if ( i != bulkData_->parallel_rank() ) {
348  const stk::mesh::BucketVector& buckets = bulkData_->buckets(stk::topology::NODE_RANK);
349  for (size_t j=0;j<buckets.size();++j) {
350  const stk::mesh::Bucket& bucket = *buckets[j];
351  if ( bucket.owned() ) {
352  for (size_t k=0;k<bucket.size();++k) {
353  stk::mesh::EntityKey key = bulkData_->entity_key(bucket[k]);
354  comm.send_buffer(i).pack<stk::mesh::EntityKey>(key);
355  }
356  }
357  }
358  }
359  }
360 
361  if (phase == 0 ) {
362  comm.allocate_buffers();
363  }
364  else {
365  comm.communicate();
366  }
367  }
368 
369  for (int i=0;i<bulkData_->parallel_size();++i) {
370  if ( i != bulkData_->parallel_rank() ) {
371  while(comm.recv_buffer(i).remaining()) {
372  stk::mesh::EntityKey key;
373  comm.recv_buffer(i).unpack<stk::mesh::EntityKey>(key);
374  stk::mesh::Entity node = bulkData_->get_entity(key);
375  if ( bulkData_->is_valid(node) ) {
376  bulkData_->add_node_sharing(node, i);
377  }
378  }
379  }
380  }
381 
382 
383  bulkData_->modification_end();
384 
387 }
388 
389 void STK_Interface::addNode(stk::mesh::EntityId gid, const std::vector<double> & coord)
390 {
391  TEUCHOS_TEST_FOR_EXCEPTION(not isModifiable(),std::logic_error,
392  "STK_Interface::addNode: STK_Interface must be modifiable to add a node");
393  TEUCHOS_TEST_FOR_EXCEPTION(coord.size()!=getDimension(),std::logic_error,
394  "STK_Interface::addNode: number of coordinates in vector must mation dimension");
395  TEUCHOS_TEST_FOR_EXCEPTION(gid==0,std::logic_error,
396  "STK_Interface::addNode: STK has STUPID restriction of no zero GIDs, pick something else");
397  stk::mesh::EntityRank nodeRank = getNodeRank();
398 
399  stk::mesh::Entity node = bulkData_->declare_entity(nodeRank,gid,nodesPartVec_);
400 
401  // set coordinate vector
402  double * fieldCoords = stk::mesh::field_data(*coordinatesField_,node);
403  for(std::size_t i=0;i<coord.size();++i)
404  fieldCoords[i] = coord[i];
405 }
406 
407 void STK_Interface::addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset)
408 {
409  std::vector<stk::mesh::Part*> sidesetV;
410  sidesetV.push_back(sideset);
411 
412  bulkData_->change_entity_parts(entity,sidesetV);
413 }
414 
415 void STK_Interface::addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset)
416 {
417  std::vector<stk::mesh::Part*> nodesetV;
418  nodesetV.push_back(nodeset);
419 
420  bulkData_->change_entity_parts(entity,nodesetV);
421 }
422 
423 void STK_Interface::addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block)
424 {
425  std::vector<stk::mesh::Part*> blockVec;
426  blockVec.push_back(block);
427 
428  stk::mesh::EntityRank elementRank = getElementRank();
429  stk::mesh::EntityRank nodeRank = getNodeRank();
430  stk::mesh::Entity element = bulkData_->declare_entity(elementRank,ed->getGID(),blockVec);
431 
432  // build relations that give the mesh structure
433  const std::vector<stk::mesh::EntityId> & nodes = ed->getNodes();
434  for(std::size_t i=0;i<nodes.size();++i) {
435  // add element->node relation
436  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodes[i]);
437  TEUCHOS_ASSERT(bulkData_->is_valid(node));
438  bulkData_->declare_relation(element,node,i);
439  }
440 
441  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
442  procId[0] = Teuchos::as<ProcIdData>(procRank_);
443 }
444 
445 
447 {
448  // loop over elements
449  stk::mesh::EntityRank edgeRank = getEdgeRank();
450  std::vector<stk::mesh::Entity> localElmts;
451  getMyElements(localElmts);
452  std::vector<stk::mesh::Entity>::const_iterator itr;
453  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
454  stk::mesh::Entity element = (*itr);
455  stk::mesh::EntityId gid = bulkData_->identifier(element);
456  std::vector<stk::mesh::EntityId> subcellIds;
457  getSubcellIndices(edgeRank,gid,subcellIds);
458 
459  for(std::size_t i=0;i<subcellIds.size();++i) {
460  stk::mesh::Entity edge = bulkData_->get_entity(edgeRank,subcellIds[i]);
461  stk::mesh::Entity const* relations = bulkData_->begin_nodes(edge);
462 
463  double * node_coord_1 = stk::mesh::field_data(*coordinatesField_,relations[0]);
464  double * node_coord_2 = stk::mesh::field_data(*coordinatesField_,relations[1]);
465 
466  // set coordinate vector
467  double * edgeCoords = stk::mesh::field_data(*edgesField_,edge);
468  for(std::size_t i=0;i<getDimension();++i)
469  edgeCoords[i] = (node_coord_1[i]+node_coord_2[i])/2.0;
470  }
471  }
472 }
473 
475 {
476  // loop over elements
477  stk::mesh::EntityRank faceRank = getFaceRank();
478  std::vector<stk::mesh::Entity> localElmts;
479  getMyElements(localElmts);
480  std::vector<stk::mesh::Entity>::const_iterator itr;
481  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
482  stk::mesh::Entity element = (*itr);
483  stk::mesh::EntityId gid = elementGlobalId(element);
484  std::vector<stk::mesh::EntityId> subcellIds;
485  getSubcellIndices(faceRank,gid,subcellIds);
486 
487  for(std::size_t i=0;i<subcellIds.size();++i) {
488  stk::mesh::Entity face = bulkData_->get_entity(faceRank,subcellIds[i]);
489  stk::mesh::Entity const* relations = bulkData_->begin_nodes(face);
490  const size_t num_relations = bulkData_->num_nodes(face);
491 
492  // set coordinate vector
493  double * faceCoords = stk::mesh::field_data(*facesField_,face);
494  for(std::size_t i=0;i<getDimension();++i){
495  faceCoords[i] = 0.0;
496  for(std::size_t j=0;j<num_relations;++j)
497  faceCoords[i] += stk::mesh::field_data(*coordinatesField_,relations[j])[i];
498  faceCoords[i] /= double(num_relations);
499  }
500  }
501  }
502 }
503 
504 stk::mesh::Entity STK_Interface::findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
505 {
506  const size_t num_rels = bulkData_->num_connectivity(src, tgt_rank);
507  stk::mesh::Entity const* relations = bulkData_->begin(src, tgt_rank);
508  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData_->begin_ordinals(src, tgt_rank);
509  for (size_t i = 0; i < num_rels; ++i) {
510  if (ordinals[i] == static_cast<stk::mesh::ConnectivityOrdinal>(rel_id)) {
511  return relations[i];
512  }
513  }
514 
515  return stk::mesh::Entity();
516 }
517 
519 //
520 // writeToExodus()
521 //
523 void
526  const std::string& filename)
527 {
528  setupExodusFile(filename);
529  writeToExodus(0.0);
530 } // end of writeToExodus()
531 
533 //
534 // setupExodusFile()
535 //
537 void
540  const std::string& filename)
541 {
542  using std::runtime_error;
543  using stk::io::StkMeshIoBroker;
544  using stk::mesh::FieldVector;
545  using stk::ParallelMachine;
546  using Teuchos::rcp;
547  PANZER_FUNC_TIME_MONITOR("STK_Interface::setupExodusFile(filename)");
548 #ifdef PANZER_HAVE_IOSS
550  ParallelMachine comm = *mpiComm_->getRawMpiComm();
551  meshData_ = rcp(new StkMeshIoBroker(comm));
552  meshData_->set_bulk_data(bulkData_);
553  meshIndex_ = meshData_->create_output_mesh(filename, stk::io::WRITE_RESULTS);
554  const FieldVector& fields = metaData_->get_fields();
555  for (size_t i(0); i < fields.size(); ++i) {
556  // Do NOT add MESH type stk fields to exodus io, but do add everything
557  // else. This allows us to avoid having to catch a throw for
558  // re-registering coordinates, sidesets, etc... Note that some
559  // fields like LOAD_BAL don't always have a role assigned, so for
560  // roles that point to null, we need to register them as well.
561  auto role = stk::io::get_field_role(*fields[i]);
562  if (role != nullptr) {
563  if (*role != Ioss::Field::MESH)
564  meshData_->add_field(meshIndex_, *fields[i]);
565  } else {
566  meshData_->add_field(meshIndex_, *fields[i]);
567  }
568  }
569 #else
570  TEUCHOS_ASSERT(false)
571 #endif
572 } // end of setupExodusFile()
573 
575 //
576 // writeToExodus()
577 //
579 void
582  double timestep)
583 {
584  using std::cout;
585  using std::endl;
586  using Teuchos::FancyOStream;
587  using Teuchos::rcpFromRef;
588  PANZER_FUNC_TIME_MONITOR("STK_Interface::writeToExodus(timestep)");
589 #ifdef PANZER_HAVE_IOSS
590  if (not meshData_.is_null())
591  {
592  currentStateTime_ = timestep;
593  globalToExodus(GlobalVariable::ADD);
594  meshData_->begin_output_step(meshIndex_, currentStateTime_);
595  meshData_->write_defined_output_fields(meshIndex_);
596  globalToExodus(GlobalVariable::WRITE);
597  meshData_->end_output_step(meshIndex_);
598  }
599  else // if (meshData_.is_null())
600  {
601  FancyOStream out(rcpFromRef(cout));
602  out.setOutputToRootOnly(0);
603  out << "WARNING: Exodus I/O has been disabled or not setup properly; "
604  << "not writing to Exodus." << endl;
605  } // end if (meshData_.is_null()) or not
606 #else
607  TEUCHOS_ASSERT(false);
608 #endif
609 } // end of writeToExodus()
610 
612 //
613 // globalToExodus()
614 //
616 void
617 STK_Interface::
618 globalToExodus(
619  const GlobalVariable& flag)
620 {
621  using std::invalid_argument;
622  using std::string;
623  using Teuchos::Array;
624 
625  // Loop over all the global variables to be added to the Exodus output file.
626  // For each global variable, we determine the data type, and then add or
627  // write it accordingly, depending on the value of flag.
628  for (auto i = globalData_.begin(); i != globalData_.end(); ++i)
629  {
630  const string& name = globalData_.name(i);
631 
632  // Integers.
633  if (globalData_.isType<int>(name))
634  {
635  const auto& value = globalData_.get<int>(name);
636  if (flag == GlobalVariable::ADD)
637  {
638  try
639  {
640  meshData_->add_global(meshIndex_, name, value,
641  stk::util::ParameterType::INTEGER);
642  }
643  catch (...)
644  {
645  return;
646  }
647  }
648  else // if (flag == GlobalVariable::WRITE)
649  meshData_->write_global(meshIndex_, name, value,
650  stk::util::ParameterType::INTEGER);
651  }
652 
653  // Doubles.
654  else if (globalData_.isType<double>(name))
655  {
656  const auto& value = globalData_.get<double>(name);
657  if (flag == GlobalVariable::ADD)
658  {
659  try
660  {
661  meshData_->add_global(meshIndex_, name, value,
662  stk::util::ParameterType::DOUBLE);
663  }
664  catch (...)
665  {
666  return;
667  }
668  }
669  else // if (flag == GlobalVariable::WRITE)
670  meshData_->write_global(meshIndex_, name, value,
671  stk::util::ParameterType::DOUBLE);
672  }
673 
674  // Vectors of integers.
675  else if (globalData_.isType<Array<int>>(name))
676  {
677  const auto& value = globalData_.get<Array<int>>(name).toVector();
678  if (flag == GlobalVariable::ADD)
679  {
680  try
681  {
682  meshData_->add_global(meshIndex_, name, value,
683  stk::util::ParameterType::INTEGERVECTOR);
684  }
685  catch (...)
686  {
687  return;
688  }
689  }
690  else // if (flag == GlobalVariable::WRITE)
691  meshData_->write_global(meshIndex_, name, value,
692  stk::util::ParameterType::INTEGERVECTOR);
693  }
694 
695  // Vectors of doubles.
696  else if (globalData_.isType<Array<double>>(name))
697  {
698  const auto& value = globalData_.get<Array<double>>(name).toVector();
699  if (flag == GlobalVariable::ADD)
700  {
701  try
702  {
703  meshData_->add_global(meshIndex_, name, value,
704  stk::util::ParameterType::DOUBLEVECTOR);
705  }
706  catch (...)
707  {
708  return;
709  }
710  }
711  else // if (flag == GlobalVariable::WRITE)
712  meshData_->write_global(meshIndex_, name, value,
713  stk::util::ParameterType::DOUBLEVECTOR);
714  }
715 
716  // If the data type is something else, throw an exception.
717  else
718  TEUCHOS_TEST_FOR_EXCEPTION(true, invalid_argument,
719  "STK_Interface::globalToExodus(): The global variable to be added " \
720  "to the Exodus output file is of an invalid type. Valid types are " \
721  "int and double, along with std::vectors of those types.")
722  } // end loop over globalData_
723 } // end of globalToExodus()
724 
726 //
727 // addGlobalToExodus()
728 //
730 void
733  const std::string& key,
734  const int& value)
735 {
736  globalData_.set(key, value);
737 } // end of addGlobalToExodus()
738 
740 //
741 // addGlobalToExodus()
742 //
744 void
747  const std::string& key,
748  const double& value)
749 {
750  globalData_.set(key, value);
751 } // end of addGlobalToExodus()
752 
754 //
755 // addGlobalToExodus()
756 //
758 void
761  const std::string& key,
762  const std::vector<int>& value)
763 {
764  using Teuchos::Array;
765  globalData_.set(key, Array<int>(value));
766 } // end of addGlobalToExodus()
767 
769 //
770 // addGlobalToExodus()
771 //
773 void
776  const std::string& key,
777  const std::vector<double>& value)
778 {
779  using Teuchos::Array;
780  globalData_.set(key, Array<double>(value));
781 } // end of addGlobalToExodus()
782 
784 {
785  #ifdef PANZER_HAVE_IOSS
786  return true;
787  #else
788  return false;
789  #endif
790 }
791 
792 void STK_Interface::getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const
793 {
794  stk::mesh::EntityRank nodeRank = getNodeRank();
795 
796  // get all relations for node
797  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodeId);
798  const size_t numElements = bulkData_->num_elements(node);
799  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
800 
801  // extract elements sharing nodes
802  elements.insert(elements.end(), relations, relations + numElements);
803 }
804 
805 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
806  std::vector<int> & relIds) const
807 {
808  // get all relations for node
809  const size_t numElements = bulkData_->num_elements(node);
810  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
811  stk::mesh::ConnectivityOrdinal const* rel_ids = bulkData_->begin_element_ordinals(node);
812 
813  // extract elements sharing nodes
814  for (size_t i = 0; i < numElements; ++i) {
815  stk::mesh::Entity element = relations[i];
816 
817  // if owned by this processor
818  if(bulkData_->parallel_owner_rank(element) == static_cast<int>(procRank_)) {
819  elements.push_back(element);
820  relIds.push_back(rel_ids[i]);
821  }
822  }
823 }
824 
825 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
826  std::vector<int> & relIds, unsigned int matchType) const
827 {
828  stk::mesh::EntityRank rank;
829  if(matchType == 0)
830  rank = getNodeRank();
831  else if(matchType == 1)
832  rank = getEdgeRank();
833  else if(matchType == 2)
834  rank = getFaceRank();
835  else
836  TEUCHOS_ASSERT(false);
837 
838  stk::mesh::Entity node = bulkData_->get_entity(rank,nodeId);
839 
840  getOwnedElementsSharingNode(node,elements,relIds);
841 }
842 
843 void STK_Interface::getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeIds,std::vector<stk::mesh::Entity> & elements) const
844 {
845  std::vector<stk::mesh::Entity> current;
846 
847  getElementsSharingNode(nodeIds[0],current); // fill it with elements touching first node
848  std::sort(current.begin(),current.end()); // sort for intersection on the pointer
849 
850  // find intersection with remaining nodes
851  for(std::size_t n=1;n<nodeIds.size();++n) {
852  // get elements associated with next node
853  std::vector<stk::mesh::Entity> nextNode;
854  getElementsSharingNode(nodeIds[n],nextNode); // fill it with elements touching first node
855  std::sort(nextNode.begin(),nextNode.end()); // sort for intersection on the pointer ID
856 
857  // intersect next node elements with current element list
858  std::vector<stk::mesh::Entity> intersection(std::min(nextNode.size(),current.size()));
859  std::vector<stk::mesh::Entity>::const_iterator endItr
860  = std::set_intersection(current.begin(),current.end(),
861  nextNode.begin(),nextNode.end(),
862  intersection.begin());
863  std::size_t newLength = endItr-intersection.begin();
864  intersection.resize(newLength);
865 
866  // store intersection
867  current.clear();
868  current = intersection;
869  }
870 
871  // return the elements computed
872  elements = current;
873 }
874 
875 void STK_Interface::getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const
876 {
877  stk::mesh::Entity const* nodeRel = getBulkData()->begin_nodes(element);
878  const size_t numNodes = getBulkData()->num_nodes(element);
879 
880  nodeIds.reserve(numNodes);
881  for(size_t i = 0; i < numNodes; ++i) {
882  nodeIds.push_back(elementGlobalId(nodeRel[i]));
883  }
884 }
885 
887 {
888  entityCounts_.clear();
889  stk::mesh::comm_mesh_counts(*bulkData_,entityCounts_);
890 }
891 
893 {
894  // developed to mirror "comm_mesh_counts" in stk_mesh/base/Comm.cpp
895 
896  const unsigned entityRankCount = metaData_->entity_rank_count();
897  const size_t commCount = 10; // entityRankCount
898 
899  TEUCHOS_ASSERT(entityRankCount<10);
900 
901  // stk::ParallelMachine mach = bulkData_->parallel();
902  stk::ParallelMachine mach = *mpiComm_->getRawMpiComm();
903 
904  std::vector<stk::mesh::EntityId> local(commCount,0);
905 
906  // determine maximum ID for this processor for each entity type
907  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
908  for(stk::mesh::EntityRank i=stk::topology::NODE_RANK; i<entityRankCount; ++i) {
909  std::vector<stk::mesh::Entity> entities;
910 
911  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(i),entities);
912 
913  // determine maximum ID for this processor
914  std::vector<stk::mesh::Entity>::const_iterator itr;
915  for(itr=entities.begin();itr!=entities.end();++itr) {
916  stk::mesh::EntityId id = bulkData_->identifier(*itr);
917  if(id>local[i])
918  local[i] = id;
919  }
920  }
921 
922  // get largest IDs across processors
923  stk::all_reduce(mach,stk::ReduceMax<10>(&local[0]));
924  maxEntityId_.assign(local.begin(),local.begin()+entityRankCount+1);
925 }
926 
927 std::size_t STK_Interface::getEntityCounts(unsigned entityRank) const
928 {
929  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=entityCounts_.size(),std::logic_error,
930  "STK_Interface::getEntityCounts: Entity counts do not include rank: " << entityRank);
931 
932  return entityCounts_[entityRank];
933 }
934 
935 stk::mesh::EntityId STK_Interface::getMaxEntityId(unsigned entityRank) const
936 {
937  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=maxEntityId_.size(),std::logic_error,
938  "STK_Interface::getMaxEntityId: Max entity ids do not include rank: " << entityRank);
939 
940  return maxEntityId_[entityRank];
941 }
942 
944 {
945  stk::mesh::PartVector emptyPartVector;
946  stk::mesh::create_adjacent_entities(*bulkData_,emptyPartVector);
947 
950 
951  addEdges();
952  if (dimension_ > 2)
953  addFaces();
954 }
955 
956 const double * STK_Interface::getNodeCoordinates(stk::mesh::EntityId nodeId) const
957 {
958  stk::mesh::Entity node = bulkData_->get_entity(getNodeRank(),nodeId);
959  return stk::mesh::field_data(*coordinatesField_,node);
960 }
961 
962 const double * STK_Interface::getNodeCoordinates(stk::mesh::Entity node) const
963 {
964  return stk::mesh::field_data(*coordinatesField_,node);
965 }
966 
967 void STK_Interface::getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
968  std::vector<stk::mesh::EntityId> & subcellIds) const
969 {
970  stk::mesh::EntityRank elementRank = getElementRank();
971  stk::mesh::Entity cell = bulkData_->get_entity(elementRank,elementId);
972 
973  TEUCHOS_TEST_FOR_EXCEPTION(!bulkData_->is_valid(cell),std::logic_error,
974  "STK_Interface::getSubcellIndices: could not find element requested (GID = " << elementId << ")");
975 
976  const size_t numSubcells = bulkData_->num_connectivity(cell, static_cast<stk::mesh::EntityRank>(entityRank));
977  stk::mesh::Entity const* subcells = bulkData_->begin(cell, static_cast<stk::mesh::EntityRank>(entityRank));
978  subcellIds.clear();
979  subcellIds.resize(numSubcells,0);
980 
981  // loop over relations and fill subcell vector
982  for(size_t i = 0; i < numSubcells; ++i) {
983  stk::mesh::Entity subcell = subcells[i];
984  subcellIds[i] = bulkData_->identifier(subcell);
985  }
986 }
987 
988 void STK_Interface::getMyElements(std::vector<stk::mesh::Entity> & elements) const
989 {
990  // setup local ownership
991  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
992 
993  // grab elements
994  stk::mesh::EntityRank elementRank = getElementRank();
995  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(elementRank),elements);
996 }
997 
998 void STK_Interface::getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
999 {
1000  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
1001 
1002  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
1003 
1004  // setup local ownership
1005  // stk::mesh::Selector block = *elementBlock;
1006  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & (*elementBlock);
1007 
1008  // grab elements
1009  stk::mesh::EntityRank elementRank = getElementRank();
1010  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(elementRank),elements);
1011 }
1012 
1013 void STK_Interface::getNeighborElements(std::vector<stk::mesh::Entity> & elements) const
1014 {
1015  // setup local ownership
1016  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part());
1017 
1018  // grab elements
1019  stk::mesh::EntityRank elementRank = getElementRank();
1020  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
1021 }
1022 
1023 void STK_Interface::getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
1024 {
1025  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
1026 
1027  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
1028 
1029  // setup local ownership
1030  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part()) & (*elementBlock);
1031 
1032  // grab elements
1033  stk::mesh::EntityRank elementRank = getElementRank();
1034  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
1035 }
1036 
1037 void STK_Interface::getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
1038 {
1039  stk::mesh::Part * sidePart = getSideset(sideName);
1040  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
1041  "Unknown side set \"" << sideName << "\"");
1042 
1043  stk::mesh::Selector side = *sidePart;
1044  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & side;
1045 
1046  // grab elements
1047  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
1048 }
1049 
1050 void STK_Interface::getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
1051 {
1052  stk::mesh::Part * sidePart = getSideset(sideName);
1053  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1055  "Unknown side set \"" << sideName << "\"");
1057  "Unknown element block \"" << blockName << "\"");
1058 
1059  stk::mesh::Selector side = *sidePart;
1060  stk::mesh::Selector block = *elmtPart;
1061  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & side;
1062 
1063  // grab elements
1064  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
1065 }
1066 
1067 void STK_Interface::getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
1068 {
1069  stk::mesh::Part * sidePart = getSideset(sideName);
1070  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
1071  "Unknown side set \"" << sideName << "\"");
1072 
1073  stk::mesh::Selector side = *sidePart;
1074 
1075  // grab elements
1076  stk::mesh::get_selected_entities(side,bulkData_->buckets(getSideRank()),sides);
1077 }
1078 
1079 void STK_Interface::getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
1080 {
1081  stk::mesh::Part * sidePart = getSideset(sideName);
1082  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1084  "Unknown side set \"" << sideName << "\"");
1086  "Unknown element block \"" << blockName << "\"");
1087 
1088  stk::mesh::Selector side = *sidePart;
1089  stk::mesh::Selector block = *elmtPart;
1090  stk::mesh::Selector sideBlock = block & side;
1091 
1092  // grab elements
1093  stk::mesh::get_selected_entities(sideBlock,bulkData_->buckets(getSideRank()),sides);
1094 }
1095 
1096 void STK_Interface::getMyNodes(const std::string & nodesetName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const
1097 {
1098  stk::mesh::Part * nodePart = getNodeset(nodesetName);
1099  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
1101  "Unknown node set \"" << nodesetName << "\"");
1103  "Unknown element block \"" << blockName << "\"");
1104 
1105  stk::mesh::Selector nodeset = *nodePart;
1106  stk::mesh::Selector block = *elmtPart;
1107  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & nodeset;
1108 
1109  // grab elements
1110  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getNodeRank()),nodes);
1111 }
1112 
1113 void STK_Interface::getElementBlockNames(std::vector<std::string> & names) const
1114 {
1115  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
1116 
1117  names.clear();
1118 
1119  // fill vector with automagically ordered string values
1120  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr; // Element blocks
1121  for(blkItr=elementBlocks_.begin();blkItr!=elementBlocks_.end();++blkItr)
1122  names.push_back(blkItr->first);
1123 }
1124 
1125 void STK_Interface::getSidesetNames(std::vector<std::string> & names) const
1126 {
1127  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
1128 
1129  names.clear();
1130 
1131  // fill vector with automagically ordered string values
1132  std::map<std::string, stk::mesh::Part*>::const_iterator sideItr; // Element blocks
1133  for(sideItr=sidesets_.begin();sideItr!=sidesets_.end();++sideItr)
1134  names.push_back(sideItr->first);
1135 }
1136 
1137 void STK_Interface::getNodesetNames(std::vector<std::string> & names) const
1138 {
1139  names.clear();
1140 
1141  // fill vector with automagically ordered string values
1142  std::map<std::string, stk::mesh::Part*>::const_iterator nodeItr; // Element blocks
1143  for(nodeItr=nodesets_.begin();nodeItr!=nodesets_.end();++nodeItr)
1144  names.push_back(nodeItr->first);
1145 }
1146 
1147 std::size_t STK_Interface::elementLocalId(stk::mesh::Entity elmt) const
1148 {
1149  return elementLocalId(bulkData_->identifier(elmt));
1150  // const std::size_t * fieldCoords = stk::mesh::field_data(*localIdField_,*elmt);
1151  // return fieldCoords[0];
1152 }
1153 
1154 std::size_t STK_Interface::elementLocalId(stk::mesh::EntityId gid) const
1155 {
1156  // stk::mesh::EntityRank elementRank = getElementRank();
1157  // stk::mesh::Entity elmt = bulkData_->get_entity(elementRank,gid);
1158  // TEUCHOS_ASSERT(elmt->owner_rank()==procRank_);
1159  // return elementLocalId(elmt);
1160  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localIDHash_.find(gid);
1161  TEUCHOS_ASSERT(itr!=localIDHash_.end());
1162  return itr->second;
1163 }
1164 
1165 
1166 std::string STK_Interface::containingBlockId(stk::mesh::Entity elmt) const
1167 {
1168  for(const auto & eb_pair : elementBlocks_)
1169  if(bulkData_->bucket(elmt).member(*(eb_pair.second)))
1170  return eb_pair.first;
1171  return "";
1172 }
1173 
1174 stk::mesh::Field<double> * STK_Interface::getSolutionField(const std::string & fieldName,
1175  const std::string & blockId) const
1176 {
1177  // look up field in map
1178  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1179  iter = fieldNameToSolution_.find(std::make_pair(fieldName,blockId));
1180 
1181  // check to make sure field was actually found
1182  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToSolution_.end(),std::runtime_error,
1183  "Solution field name \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1184 
1185  return iter->second;
1186 }
1187 
1188 stk::mesh::Field<double> * STK_Interface::getCellField(const std::string & fieldName,
1189  const std::string & blockId) const
1190 {
1191  // look up field in map
1192  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
1193  iter = fieldNameToCellField_.find(std::make_pair(fieldName,blockId));
1194 
1195  // check to make sure field was actually found
1196  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToCellField_.end(),std::runtime_error,
1197  "Cell field named \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1198 
1199  return iter->second;
1200 }
1201 
1203 {
1204  using Teuchos::RCP;
1205  using Teuchos::rcp;
1206 
1207  if(orderedElementVector_==Teuchos::null) {
1208  // safe because essentially this is a call to modify a mutable object
1209  const_cast<STK_Interface*>(this)->buildLocalElementIDs();
1210  }
1211 
1213 }
1214 
1215 void STK_Interface::addElementBlock(const std::string & name,const CellTopologyData * ctData)
1216 {
1218 
1219  stk::mesh::Part * block = metaData_->get_part(name);
1220  if(block==0) {
1221  block = &metaData_->declare_part_with_topology(name, stk::mesh::get_topology(shards::CellTopology(ctData)));
1222  }
1223 
1224  // construct cell topology object for this block
1226  = Teuchos::rcp(new shards::CellTopology(ctData));
1227 
1228  // add element block part and cell topology
1229  elementBlocks_.insert(std::make_pair(name,block));
1230  elementBlockCT_.insert(std::make_pair(name,ct));
1231 }
1232 
1234 {
1235  dimension_ = metaData_->spatial_dimension();
1236 
1237  // declare coordinates and node parts
1238  coordinatesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::NODE_RANK, coordsString);
1239  edgesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::EDGE_RANK, edgesString);
1240  if (dimension_ > 2)
1241  facesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::FACE_RANK, facesString);
1242  processorIdField_ = &metaData_->declare_field<ProcIdFieldType>(stk::topology::ELEMENT_RANK, "PROC_ID");
1243  loadBalField_ = &metaData_->declare_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, "LOAD_BAL");
1244 
1245  // stk::mesh::put_field( *coordinatesField_ , getNodeRank() , metaData_->universal_part() );
1246 
1247  nodesPart_ = &metaData_->declare_part(nodesString,getNodeRank());
1248  nodesPartVec_.push_back(nodesPart_);
1249  edgesPart_ = &metaData_->declare_part(edgesString,getEdgeRank());
1250  edgesPartVec_.push_back(edgesPart_);
1251  if (dimension_ > 2) {
1252  facesPart_ = &metaData_->declare_part(facesString,getFaceRank());
1253  facesPartVec_.push_back(facesPart_);
1254  }
1255 }
1256 
1258 {
1259  currentLocalId_ = 0;
1260 
1261  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
1262 
1263  // might be better (faster) to do this by buckets
1264  std::vector<stk::mesh::Entity> elements;
1265  getMyElements(elements);
1266 
1267  for(std::size_t index=0;index<elements.size();++index) {
1268  stk::mesh::Entity element = elements[index];
1269 
1270  // set processor rank
1271  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1272  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1273 
1274  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1275 
1276  currentLocalId_++;
1277  }
1278 
1279  // copy elements into the ordered element vector
1280  orderedElementVector_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>(elements));
1281 
1282  elements.clear();
1283  getNeighborElements(elements);
1284 
1285  for(std::size_t index=0;index<elements.size();++index) {
1286  stk::mesh::Entity element = elements[index];
1287 
1288  // set processor rank
1289  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1290  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1291 
1292  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1293 
1294  currentLocalId_++;
1295  }
1296 
1297  orderedElementVector_->insert(orderedElementVector_->end(),elements.begin(),elements.end());
1298 }
1299 
1301 {
1302  std::vector<std::string> names;
1303  getElementBlockNames(names);
1304 
1305  for(std::size_t b=0;b<names.size();b++) {
1306  // find user specified block weight, otherwise use 1.0
1307  std::map<std::string,double>::const_iterator bw_itr = blockWeights_.find(names[b]);
1308  double blockWeight = (bw_itr!=blockWeights_.end()) ? bw_itr->second : 1.0;
1309 
1310  std::vector<stk::mesh::Entity> elements;
1311  getMyElements(names[b],elements);
1312 
1313  for(std::size_t index=0;index<elements.size();++index) {
1314  // set local element ID
1315  double * loadBal = stk::mesh::field_data(*loadBalField_,elements[index]);
1316  loadBal[0] = blockWeight;
1317  }
1318  }
1319 }
1320 
1321 bool
1322 STK_Interface::isMeshCoordField(const std::string & eBlock,
1323  const std::string & fieldName,
1324  int & axis) const
1325 {
1326  std::map<std::string,std::vector<std::string> >::const_iterator blkItr = meshCoordFields_.find(eBlock);
1327  if(blkItr==meshCoordFields_.end()) {
1328  return false;
1329  }
1330 
1331  axis = 0;
1332  for(axis=0;axis<Teuchos::as<int>(blkItr->second.size());axis++) {
1333  if(blkItr->second[axis]==fieldName)
1334  break; // found axis, break
1335  }
1336 
1337  if(axis>=Teuchos::as<int>(blkItr->second.size()))
1338  return false;
1339 
1340  return true;
1341 }
1342 
1343 std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
1345 {
1347  Teuchos::RCP<std::vector<unsigned int > > type_vec = rcp(new std::vector<unsigned int>);
1348  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & matchers = getPeriodicBCVector();
1349 
1350  // build up the vectors by looping over the matched pair
1351  for(std::size_t m=0;m<matchers.size();m++){
1352  vec = matchers[m]->getMatchedPair(*this,vec);
1353  unsigned int type;
1354  if(matchers[m]->getType() == "coord")
1355  type = 0;
1356  else if(matchers[m]->getType() == "edge")
1357  type = 1;
1358  else if(matchers[m]->getType() == "face")
1359  type = 2;
1360  else
1361  TEUCHOS_ASSERT(false);
1362  type_vec->insert(type_vec->begin(),vec->size()-type_vec->size(),type);
1363  }
1364 
1365  return std::make_pair(vec,type_vec);
1366 }
1367 
1368 bool STK_Interface::validBlockId(const std::string & blockId) const
1369 {
1370  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr = elementBlocks_.find(blockId);
1371 
1372  return blkItr!=elementBlocks_.end();
1373 }
1374 
1375 void STK_Interface::print(std::ostream & os) const
1376 {
1377  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1378 
1379  getElementBlockNames(blockNames);
1380  getSidesetNames(sidesetNames);
1381  getNodesetNames(nodesetNames);
1382 
1383  os << "STK Mesh data:\n";
1384  os << " Spatial dim = " << getDimension() << "\n";
1385  if(getDimension()==2)
1386  os << " Entity counts (Nodes, Edges, Cells) = ( "
1387  << getEntityCounts(getNodeRank()) << ", "
1388  << getEntityCounts(getEdgeRank()) << ", "
1389  << getEntityCounts(getElementRank()) << " )\n";
1390  else if(getDimension()==3)
1391  os << " Entity counts (Nodes, Edges, Faces, Cells) = ( "
1392  << getEntityCounts(getNodeRank()) << ", "
1393  << getEntityCounts(getEdgeRank()) << ", "
1394  << getEntityCounts(getSideRank()) << ", "
1395  << getEntityCounts(getElementRank()) << " )\n";
1396  else
1397  os << " Entity counts (Nodes, Cells) = ( "
1398  << getEntityCounts(getNodeRank()) << ", "
1399  << getEntityCounts(getElementRank()) << " )\n";
1400 
1401  os << " Element blocks = ";
1402  for(std::size_t i=0;i<blockNames.size();i++)
1403  os << "\"" << blockNames[i] << "\" ";
1404  os << "\n";
1405  os << " Sidesets = ";
1406  for(std::size_t i=0;i<sidesetNames.size();i++)
1407  os << "\"" << sidesetNames[i] << "\" ";
1408  os << "\n";
1409  os << " Nodesets = ";
1410  for(std::size_t i=0;i<nodesetNames.size();i++)
1411  os << "\"" << nodesetNames[i] << "\" ";
1412  os << std::endl;
1413 
1414  // print out periodic boundary conditions
1415  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1416  = getPeriodicBCVector();
1417  if(bcVector.size()!=0) {
1418  os << " Periodic BCs:\n";
1419  for(std::size_t i=0;i<bcVector.size();i++)
1420  os << " " << bcVector[i]->getString() << "\n";
1421  os << std::endl;
1422  }
1423 }
1424 
1425 void STK_Interface::printMetaData(std::ostream & os) const
1426 {
1427  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1428 
1429  getElementBlockNames(blockNames);
1430  getSidesetNames(sidesetNames);
1431  getNodesetNames(nodesetNames);
1432 
1433  os << "STK Meta data:\n";
1434  os << " Element blocks = ";
1435  for(std::size_t i=0;i<blockNames.size();i++)
1436  os << "\"" << blockNames[i] << "\" ";
1437  os << "\n";
1438  os << " Sidesets = ";
1439  for(std::size_t i=0;i<sidesetNames.size();i++)
1440  os << "\"" << sidesetNames[i] << "\" ";
1441  os << "\n";
1442  os << " Nodesets = ";
1443  for(std::size_t i=0;i<nodesetNames.size();i++)
1444  os << "\"" << nodesetNames[i] << "\" ";
1445  os << std::endl;
1446 
1447  // print out periodic boundary conditions
1448  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1449  = getPeriodicBCVector();
1450  if(bcVector.size()!=0) {
1451  os << " Periodic BCs:\n";
1452  for(std::size_t i=0;i<bcVector.size();i++)
1453  os << " " << bcVector[i]->getString() << "\n";
1454  os << std::endl;
1455  }
1456 
1457  // print all available fields in meta data
1458  os << " Fields = ";
1459  const stk::mesh::FieldVector & fv = metaData_->get_fields();
1460  for(std::size_t i=0;i<fv.size();i++)
1461  os << "\"" << fv[i]->name() << "\" ";
1462  os << std::endl;
1463 }
1464 
1466 {
1467  std::map<std::string, Teuchos::RCP<shards::CellTopology> >::const_iterator itr;
1468  itr = elementBlockCT_.find(eBlock);
1469 
1470  if(itr==elementBlockCT_.end()) {
1471  std::stringstream ss;
1472  printMetaData(ss);
1473  TEUCHOS_TEST_FOR_EXCEPTION(itr==elementBlockCT_.end(),std::logic_error,
1474  "STK_Interface::getCellTopology: No such element block \"" +eBlock +"\" available.\n\n"
1475  << "STK Meta Data follows: \n" << ss.str());
1476  }
1477 
1478  return itr->second;
1479 }
1480 
1482 {
1483  MPI_Comm newComm;
1484  const int err = MPI_Comm_dup (parallelMach, &newComm);
1485  TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error,
1486  "panzer::STK_Interface: MPI_Comm_dup failed with error \""
1487  << Teuchos::mpiErrorCodeToString (err) << "\".");
1488 
1489  return Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper (newComm,MPI_Comm_free)));
1490 }
1491 
1493 {
1494  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Rebalance not currently supported");
1495 #if 0
1496  // make sure weights have been set (a local operation)
1498 
1499  stk::mesh::Selector selector(getMetaData()->universal_part());
1500  stk::mesh::Selector owned_selector(getMetaData()->locally_owned_part());
1501 
1502  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
1503  out.setOutputToRootOnly(0);
1504  out << "Load balance before: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
1505 
1506  // perform reblance
1507  Teuchos::ParameterList graph;
1508  if(params.begin()!=params.end())
1509  graph.sublist(stk::rebalance::Zoltan::default_parameters_name()) = params;
1510  stk::rebalance::Zoltan zoltan_partition(*mpiComm_->getRawMpiComm(), getDimension(), graph);
1511  stk::rebalance::rebalance(*getBulkData(), owned_selector, &getCoordinatesField(), loadBalField_, zoltan_partition);
1512 
1513  out << "Load balance after: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
1514 
1515  currentLocalId_ = 0;
1516  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
1517 #endif
1518 }
1519 
1522 {
1523  TEUCHOS_ASSERT(this->isInitialized());
1524  return mpiComm_;
1525 }
1526 
1527 } // end namespace panzer_stk
Teuchos::RCP< ElementDescriptor > buildElementDescriptor(stk::mesh::EntityId elmtId, std::vector< stk::mesh::EntityId > &nodes)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
bool isMeshCoordField(const std::string &eBlock, const std::string &fieldName, int &axis) const
void addNodeset(const std::string &name)
Teuchos::RCP< stk::mesh::BulkData > bulkData_
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
std::vector< stk::mesh::Part * > facesPartVec_
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToSolution_
std::string containingBlockId(stk::mesh::Entity elmt) const
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
stk::mesh::EntityRank getFaceRank() const
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.
std::size_t elementLocalId(stk::mesh::Entity elmt) const
stk::mesh::Field< double > SolutionFieldType
std::vector< stk::mesh::Part * > edgesPartVec_
std::map< std::string, double > blockWeights_
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
stk::mesh::EntityRank getNodeRank() const
void addSolutionField(const std::string &fieldName, const std::string &blockId)
std::map< std::string, Teuchos::RCP< shards::CellTopology > > elementBlockCT_
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
std::map< std::string, stk::mesh::Part * > sidesets_
std::vector< std::size_t > entityCounts_
void getElementsSharingNode(stk::mesh::EntityId nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing a single node
void getSubcellIndices(unsigned entityRank, stk::mesh::EntityId elementId, std::vector< stk::mesh::EntityId > &subcellIds) const
stk::mesh::EntityRank getSideRank() const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true)
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
static const std::string nodesString
const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const
void instantiateBulkData(stk::ParallelMachine parallelMach)
bool is_null() const
void getElementBlockNames(std::vector< std::string > &names) const
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) 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_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
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
static const std::string coordsString
stk::mesh::Field< double, stk::mesh::Cartesian > VectorFieldType
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
std::size_t getEntityCounts(unsigned entityRank) const
get the global counts for the entity of specified rank
void writeToExodus(const std::string &filename)
Write this mesh and associated fields to the given output file.
std::map< std::string, stk::mesh::Part * > elementBlocks_
bool isInitialized() const
Has initialize been called on this mesh object?
std::vector< stk::mesh::Part * > nodesPartVec_
static const std::string edgesString
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void printMetaData(std::ostream &os) const
void initializeFieldsInSTK(const std::map< std::pair< std::string, std::string >, SolutionFieldType *> &nameToField, bool setupIO)
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getElementsOrderedByLID() const
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector() const
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< std::vector< stk::mesh::Entity > > orderedElementVector_
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
void addMeshCoordFields(const std::string &blockId, const std::vector< std::string > &coordField, const std::string &dispPrefix)
bool validBlockId(const std::string &blockId) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
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)
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void getSidesetNames(std::vector< std::string > &name) const
std::map< std::string, std::vector< std::string > > meshCoordFields_
static const std::string facesString
Teuchos::RCP< Teuchos::MpiComm< int > > getSafeCommunicator(stk::ParallelMachine parallelMach) const
std::pair< Teuchos::RCP< std::vector< std::pair< std::size_t, std::size_t > > >, Teuchos::RCP< std::vector< unsigned int > > > getPeriodicNodePairing() const
void getNodesetNames(std::vector< std::string > &name) const
void setupExodusFile(const std::string &filename)
Set up an output Exodus file for writing results.
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_
const VectorFieldType & getCoordinatesField() const
void print(std::ostream &os) const
stk::mesh::EntityRank getElementRank() const
stk::mesh::Field< ProcIdData > ProcIdFieldType
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
#define TEUCHOS_ASSERT(assertion_test)
void addCellField(const std::string &fieldName, const std::string &blockId)
RCP< const T > getConst() const
void getMyNodes(const std::string &sideName, const std::string &blockName, std::vector< stk::mesh::Entity > &nodes) 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
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
void getOwnedElementsSharingNode(stk::mesh::Entity node, std::vector< stk::mesh::Entity > &elements, std::vector< int > &relIds) const
stk::mesh::Part * getNodeset(const std::string &name) const
stk::mesh::Field< double > * getSolutionField(const std::string &fieldName, const std::string &blockId) const
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToCellField_
std::map< std::string, stk::mesh::Part * > nodesets_
void getNeighborElements(std::vector< stk::mesh::Entity > &elements) const