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