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 
64 #ifdef PANZER_HAVE_IOSS
65 #include <Ionit_Initializer.h>
66 #include <stk_io/IossBridge.hpp>
67 #endif
68 
70 
71 #include <set>
72 
73 using Teuchos::RCP;
74 using Teuchos::rcp;
75 
76 namespace panzer_stk {
77 
79 ElementDescriptor::ElementDescriptor(stk::mesh::EntityId gid,const std::vector<stk::mesh::EntityId> & nodes)
80  : gid_(gid), nodes_(nodes) {}
82 
86 buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes)
87 {
88  return Teuchos::rcp(new ElementDescriptor(elmtId,nodes));
89 }
90 
91 const std::string STK_Interface::coordsString = "coordinates";
92 const std::string STK_Interface::nodesString = "nodes";
93 const std::string STK_Interface::edgesString = "edges";
94 const std::string STK_Interface::facesString = "faces";
95 
97  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
98 {
99  metaData_ = rcp(new stk::mesh::MetaData());
100 }
101 
103  : dimension_(0), initialized_(false), currentLocalId_(0), initialStateTime_(0.0), currentStateTime_(0.0), useFieldCoordinates_(false)
104 {
105  metaData_ = metaData;
106 }
107 
109  : dimension_(dim), initialized_(false), currentLocalId_(0), useFieldCoordinates_(false)
110 {
111  std::vector<std::string> entity_rank_names = stk::mesh::entity_rank_names();
112  entity_rank_names.push_back("FAMILY_TREE");
113 
114  metaData_ = rcp(new stk::mesh::MetaData(dimension_,entity_rank_names));
115 
117 }
118 
119 void STK_Interface::addSideset(const std::string & name,const CellTopologyData * ctData)
120 {
123 
124  stk::mesh::Part * sideset = metaData_->get_part(name);
125  if(sideset==NULL)
126  sideset = &metaData_->declare_part_with_topology(name,
127  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
128  sidesets_.insert(std::make_pair(name,sideset));
129 }
130 
131 void STK_Interface::addNodeset(const std::string & name)
132 {
135 
136  stk::mesh::Part * nodeset = metaData_->get_part(name);
137  if(nodeset==NULL) {
138  const CellTopologyData * ctData = shards::getCellTopologyData<shards::Node>();
139  nodeset = &metaData_->declare_part_with_topology(name,
140  stk::mesh::get_topology(shards::CellTopology(ctData), dimension_));
141  }
142  nodesets_.insert(std::make_pair(name,nodeset));
143 }
144 
145 void STK_Interface::addSolutionField(const std::string & fieldName,const std::string & blockId)
146 {
149  "Unknown element block \"" << blockId << "\"");
150  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
151 
152  // add & declare field if not already added...currently assuming linears
153  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
154  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::NODE_RANK, fieldName);
155  if(field==0)
156  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::NODE_RANK, fieldName);
158  }
159 }
160 
161 void STK_Interface::addCellField(const std::string & fieldName,const std::string & blockId)
162 {
165  "Unknown element block \"" << blockId << "\"");
166  std::pair<std::string,std::string> key = std::make_pair(fieldName,blockId);
167 
168  // add & declare field if not already added...currently assuming linears
169  if(fieldNameToCellField_.find(key)==fieldNameToCellField_.end()) {
170  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, fieldName);
171  if(field==0)
172  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, fieldName);
174  }
175 }
176 
177 void STK_Interface::addMeshCoordFields(const std::string & blockId,
178  const std::vector<std::string> & coordNames,
179  const std::string & dispPrefix)
180 {
182  TEUCHOS_ASSERT(dimension_==coordNames.size());
185  "Unknown element block \"" << blockId << "\"");
186 
187  // we only allow one alternative coordinate field
188  TEUCHOS_TEST_FOR_EXCEPTION(meshCoordFields_.find(blockId)!=meshCoordFields_.end(),std::invalid_argument,
189  "STK_Interface::addMeshCoordFields: Can't set more than one set of coordinate "
190  "fields for element block \""+blockId+"\".");
191 
192  // Note that there is a distinction between the key which is used for lookups
193  // and the field that lives on the mesh, which is used for printing the displacement.
194 
195  // just copy the coordinate names
196  meshCoordFields_[blockId] = coordNames;
197 
198  // must fill in the displacement fields
199  std::vector<std::string> & dispFields = meshDispFields_[blockId];
200  dispFields.resize(dimension_);
201 
202  for(unsigned i=0;i<dimension_;i++) {
203  std::pair<std::string,std::string> key = std::make_pair(coordNames[i],blockId);
204  std::string dispName = dispPrefix+coordNames[i];
205 
206  dispFields[i] = dispName; // record this field as a
207  // displacement field
208 
209  // add & declare field if not already added...currently assuming linears
210  if(fieldNameToSolution_.find(key)==fieldNameToSolution_.end()) {
211 
212  SolutionFieldType * field = metaData_->get_field<SolutionFieldType>(stk::topology::NODE_RANK, dispName);
213  if(field==0) {
214  field = &metaData_->declare_field<SolutionFieldType>(stk::topology::NODE_RANK, dispName);
215  }
217  }
218  }
219 }
220 
221 void STK_Interface::initialize(stk::ParallelMachine parallelMach,bool setupIO)
222 {
224  TEUCHOS_ASSERT(dimension_!=0); // no zero dimensional meshes!
225 
226  if(mpiComm_==Teuchos::null)
227  mpiComm_ = getSafeCommunicator(parallelMach);
228 
229  procRank_ = stk::parallel_machine_rank(*mpiComm_->getRawMpiComm());
230 
231  // associating the field with a part: universal part!
232  stk::mesh::put_field( *coordinatesField_ , metaData_->universal_part(), getDimension());
233  stk::mesh::put_field( *edgesField_ , metaData_->universal_part(), getDimension());
234  if (dimension_ > 2)
235  stk::mesh::put_field( *facesField_ , metaData_->universal_part(), getDimension());
236  stk::mesh::put_field( *processorIdField_ , metaData_->universal_part());
237  stk::mesh::put_field( *loadBalField_ , metaData_->universal_part());
238 
241 
242 #ifdef PANZER_HAVE_IOSS
243  if(setupIO) {
244  // setup Exodus file IO
246 
247  // add element blocks
248  {
249  std::map<std::string, stk::mesh::Part*>::iterator itr;
250  for(itr=elementBlocks_.begin();itr!=elementBlocks_.end();++itr)
251  if(!stk::io::is_part_io_part(*itr->second))
252  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
253  }
254 
255  // add side sets
256  {
257  std::map<std::string, stk::mesh::Part*>::iterator itr;
258  for(itr=sidesets_.begin();itr!=sidesets_.end();++itr)
259  if(!stk::io::is_part_io_part(*itr->second))
260  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
261  }
262 
263  // add node sets
264  {
265  std::map<std::string, stk::mesh::Part*>::iterator itr;
266  for(itr=nodesets_.begin();itr!=nodesets_.end();++itr)
267  if(!stk::io::is_part_io_part(*itr->second))
268  stk::io::put_io_part_attribute(*itr->second); // this can only be called once per part
269  }
270 
271  // add nodes
272  if(!stk::io::is_part_io_part(*nodesPart_))
273  stk::io::put_io_part_attribute(*nodesPart_);
274 
275  stk::io::set_field_role(*coordinatesField_, Ioss::Field::MESH);
276  stk::io::set_field_role(*edgesField_, Ioss::Field::MESH);
277  if (dimension_ > 2)
278  stk::io::set_field_role(*facesField_, Ioss::Field::MESH);
279  stk::io::set_field_role(*processorIdField_, Ioss::Field::TRANSIENT);
280  // stk::io::set_field_role(*loadBalField_, Ioss::Field::TRANSIENT);
281  }
282 #endif
283 
284  metaData_->commit();
285  if(bulkData_==Teuchos::null)
286  instantiateBulkData(*mpiComm_->getRawMpiComm());
287 
288  initialized_ = true;
289 }
290 
291 void STK_Interface::initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
292  bool setupIO)
293 {
294  std::set<SolutionFieldType*> uniqueFields;
295  std::map<std::pair<std::string,std::string>,SolutionFieldType*>::const_iterator fieldIter;
296  for(fieldIter=nameToField.begin();fieldIter!=nameToField.end();++fieldIter)
297  uniqueFields.insert(fieldIter->second); // this makes setting up IO easier!
298 
299  {
300  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
301  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
302  stk::mesh::put_field(*(*uniqueFieldIter),metaData_->universal_part());
303  }
304 
305 #ifdef PANZER_HAVE_IOSS
306  if(setupIO) {
307  // add solution fields
308  std::set<SolutionFieldType*>::const_iterator uniqueFieldIter;
309  for(uniqueFieldIter=uniqueFields.begin();uniqueFieldIter!=uniqueFields.end();++uniqueFieldIter)
310  stk::io::set_field_role(*(*uniqueFieldIter), Ioss::Field::TRANSIENT);
311  }
312 #endif
313 }
314 
315 void STK_Interface::instantiateBulkData(stk::ParallelMachine parallelMach)
316 {
317  TEUCHOS_ASSERT(bulkData_==Teuchos::null);
318  if(mpiComm_==Teuchos::null)
319  mpiComm_ = getSafeCommunicator(parallelMach);
320 
321  bulkData_ = rcp(new stk::mesh::BulkData(*metaData_, *mpiComm_->getRawMpiComm()));
322 }
323 
325 {
326  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
327  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"beginModification\"");
328 
329  bulkData_->modification_begin();
330 }
331 
333 {
334  TEUCHOS_TEST_FOR_EXCEPTION(bulkData_==Teuchos::null,std::logic_error,
335  "STK_Interface: Must call \"initialized\" or \"instantiateBulkData\" before \"endModification\"");
336 
337  // TODO: Resolving sharing this way comes at a cost in performance. The STK
338  // team has decided that users need to declare their own sharing. We should
339  // find where shared entities are being created in Panzer and declare it.
340  // Once this is done, the extra code below can be deleted.
341 
342  stk::CommAll comm(bulkData_->parallel());
343 
344  for (int phase=0;phase<2;++phase) {
345  for (int i=0;i<bulkData_->parallel_size();++i) {
346  if ( i != bulkData_->parallel_rank() ) {
347  const stk::mesh::BucketVector& buckets = bulkData_->buckets(stk::topology::NODE_RANK);
348  for (size_t j=0;j<buckets.size();++j) {
349  const stk::mesh::Bucket& bucket = *buckets[j];
350  if ( bucket.owned() ) {
351  for (size_t k=0;k<bucket.size();++k) {
352  stk::mesh::EntityKey key = bulkData_->entity_key(bucket[k]);
353  comm.send_buffer(i).pack<stk::mesh::EntityKey>(key);
354  }
355  }
356  }
357  }
358  }
359 
360  if (phase == 0 ) {
361  comm.allocate_buffers( bulkData_->parallel_size()/4 );
362  }
363  else {
364  comm.communicate();
365  }
366  }
367 
368  for (int i=0;i<bulkData_->parallel_size();++i) {
369  if ( i != bulkData_->parallel_rank() ) {
370  while(comm.recv_buffer(i).remaining()) {
371  stk::mesh::EntityKey key;
372  comm.recv_buffer(i).unpack<stk::mesh::EntityKey>(key);
373  stk::mesh::Entity node = bulkData_->get_entity(key);
374  if ( bulkData_->is_valid(node) ) {
375  bulkData_->add_node_sharing(node, i);
376  }
377  }
378  }
379  }
380 
381 
382  bulkData_->modification_end();
383 
386 }
387 
388 void STK_Interface::addNode(stk::mesh::EntityId gid, const std::vector<double> & coord)
389 {
390  TEUCHOS_TEST_FOR_EXCEPTION(not isModifiable(),std::logic_error,
391  "STK_Interface::addNode: STK_Interface must be modifiable to add a node");
392  TEUCHOS_TEST_FOR_EXCEPTION(coord.size()!=getDimension(),std::logic_error,
393  "STK_Interface::addNode: number of coordinates in vector must mation dimension");
394  TEUCHOS_TEST_FOR_EXCEPTION(gid==0,std::logic_error,
395  "STK_Interface::addNode: STK has STUPID restriction of no zero GIDs, pick something else");
396  stk::mesh::EntityRank nodeRank = getNodeRank();
397 
398  stk::mesh::Entity node = bulkData_->declare_entity(nodeRank,gid,nodesPartVec_);
399 
400  // set coordinate vector
401  double * fieldCoords = stk::mesh::field_data(*coordinatesField_,node);
402  for(std::size_t i=0;i<coord.size();++i)
403  fieldCoords[i] = coord[i];
404 }
405 
406 void STK_Interface::addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset)
407 {
408  std::vector<stk::mesh::Part*> sidesetV;
409  sidesetV.push_back(sideset);
410 
411  bulkData_->change_entity_parts(entity,sidesetV);
412 }
413 
414 void STK_Interface::addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset)
415 {
416  std::vector<stk::mesh::Part*> nodesetV;
417  nodesetV.push_back(nodeset);
418 
419  bulkData_->change_entity_parts(entity,nodesetV);
420 }
421 
422 void STK_Interface::addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block)
423 {
424  std::vector<stk::mesh::Part*> blockVec;
425  blockVec.push_back(block);
426 
427  stk::mesh::EntityRank elementRank = getElementRank();
428  stk::mesh::EntityRank nodeRank = getNodeRank();
429  stk::mesh::Entity element = bulkData_->declare_entity(elementRank,ed->getGID(),blockVec);
430 
431  // build relations that give the mesh structure
432  const std::vector<stk::mesh::EntityId> & nodes = ed->getNodes();
433  for(std::size_t i=0;i<nodes.size();++i) {
434  // add element->node relation
435  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodes[i]);
436  TEUCHOS_ASSERT(bulkData_->is_valid(node));
437  bulkData_->declare_relation(element,node,i);
438  }
439 
440  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
441  procId[0] = Teuchos::as<ProcIdData>(procRank_);
442 }
443 
444 
446 {
447  // loop over elements
448  stk::mesh::EntityRank edgeRank = getEdgeRank();
449  std::vector<stk::mesh::Entity> localElmts;
450  getMyElements(localElmts);
451  std::vector<stk::mesh::Entity>::const_iterator itr;
452  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
453  stk::mesh::Entity element = (*itr);
454  stk::mesh::EntityId gid = bulkData_->identifier(element);
455  std::vector<stk::mesh::EntityId> subcellIds;
456  getSubcellIndices(edgeRank,gid,subcellIds);
457 
458  for(std::size_t i=0;i<subcellIds.size();++i) {
459  stk::mesh::Entity edge = bulkData_->get_entity(edgeRank,subcellIds[i]);
460  stk::mesh::Entity const* relations = bulkData_->begin_nodes(edge);
461 
462  double * node_coord_1 = stk::mesh::field_data(*coordinatesField_,relations[0]);
463  double * node_coord_2 = stk::mesh::field_data(*coordinatesField_,relations[1]);
464 
465  // set coordinate vector
466  double * edgeCoords = stk::mesh::field_data(*edgesField_,edge);
467  for(std::size_t i=0;i<getDimension();++i)
468  edgeCoords[i] = (node_coord_1[i]+node_coord_2[i])/2.0;
469  }
470  }
471 }
472 
474 {
475  // loop over elements
476  stk::mesh::EntityRank faceRank = getFaceRank();
477  std::vector<stk::mesh::Entity> localElmts;
478  getMyElements(localElmts);
479  std::vector<stk::mesh::Entity>::const_iterator itr;
480  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
481  stk::mesh::Entity element = (*itr);
482  stk::mesh::EntityId gid = elementGlobalId(element);
483  std::vector<stk::mesh::EntityId> subcellIds;
484  getSubcellIndices(faceRank,gid,subcellIds);
485 
486  for(std::size_t i=0;i<subcellIds.size();++i) {
487  stk::mesh::Entity face = bulkData_->get_entity(faceRank,subcellIds[i]);
488  stk::mesh::Entity const* relations = bulkData_->begin_nodes(face);
489  const size_t num_relations = bulkData_->num_nodes(face);
490 
491  // set coordinate vector
492  double * faceCoords = stk::mesh::field_data(*facesField_,face);
493  for(std::size_t i=0;i<getDimension();++i){
494  faceCoords[i] = 0.0;
495  for(std::size_t j=0;j<num_relations;++j)
496  faceCoords[i] += stk::mesh::field_data(*coordinatesField_,relations[j])[i];
497  faceCoords[i] /= double(num_relations);
498  }
499  }
500  }
501 }
502 
503 void STK_Interface::writeToExodus(const std::string & filename)
504 {
505  PANZER_FUNC_TIME_MONITOR("STK_Interface::writeToExodus(filename)");
506 
507  #ifdef PANZER_HAVE_IOSS
508  TEUCHOS_ASSERT(mpiComm_!=Teuchos::null);
509  stk::ParallelMachine comm = *mpiComm_->getRawMpiComm();
510 
511  stk::io::StkMeshIoBroker meshData(comm);
512  meshData.set_bulk_data(bulkData_);
513  const int meshIndex = meshData.create_output_mesh(filename, stk::io::WRITE_RESULTS);
514  const stk::mesh::FieldVector &fields = metaData_->get_fields();
515  for (size_t i=0; i < fields.size(); i++) {
516  try {
517  meshData.add_field(meshIndex, *fields[i]);
518  }
519  catch (std::runtime_error const&) { }
520  }
521 
522  meshData.process_output_request(meshIndex, 0.0);
523  #else
524  TEUCHOS_ASSERT(false);
525  #endif
526 }
527 
528 stk::mesh::Entity STK_Interface::findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
529 {
530  const size_t num_rels = bulkData_->num_connectivity(src, tgt_rank);
531  stk::mesh::Entity const* relations = bulkData_->begin(src, tgt_rank);
532  stk::mesh::ConnectivityOrdinal const* ordinals = bulkData_->begin_ordinals(src, tgt_rank);
533  for (size_t i = 0; i < num_rels; ++i) {
534  if (ordinals[i] == static_cast<stk::mesh::ConnectivityOrdinal>(rel_id)) {
535  return relations[i];
536  }
537  }
538 
539  return stk::mesh::Entity();
540 }
541 
542 void STK_Interface::setupTransientExodusFile(const std::string & filename)
543 {
544  PANZER_FUNC_TIME_MONITOR("STK_Interface::setupTransientExodusFile(filename)");
545 
546  #ifdef PANZER_HAVE_IOSS
547  TEUCHOS_ASSERT(mpiComm_!=Teuchos::null);
548  stk::ParallelMachine comm = *mpiComm_->getRawMpiComm();
549 
550  meshData_ = Teuchos::rcp(new stk::io::StkMeshIoBroker(comm));
551  meshData_->set_bulk_data(bulkData_);
552  meshIndex_ = meshData_->create_output_mesh(filename, stk::io::WRITE_RESULTS);
553  const stk::mesh::FieldVector &fields = metaData_->get_fields();
554  for (size_t i=0; i < fields.size(); i++) {
555  try {
556  meshData_->add_field(meshIndex_, *fields[i]);
557  }
558  catch (std::runtime_error const&) { }
559  }
560  #else
561  TEUCHOS_ASSERT(false);
562  #endif
563 }
564 
565 void STK_Interface::writeToExodus(double timestep)
566 {
567  PANZER_FUNC_TIME_MONITOR("STK_Interface::writeToExodus(timestep)");
568 
569  #ifdef PANZER_HAVE_IOSS
570  if(meshData_!=Teuchos::null) {
571  currentStateTime_ = timestep;
572  meshData_->process_output_request(meshIndex_, timestep);
573  }
574  else {
575  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
576  out.setOutputToRootOnly(0);
577  out << "WARNING: Exodus I/O has been disabled or not setup properly, not writing to Exodus" << std::endl;
578  }
579  #else
580  TEUCHOS_ASSERT(false);
581  #endif
582 }
583 
585 {
586  #ifdef PANZER_HAVE_IOSS
587  return true;
588  #else
589  return false;
590  #endif
591 }
592 
593 void STK_Interface::getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const
594 {
595  stk::mesh::EntityRank nodeRank = getNodeRank();
596 
597  // get all relations for node
598  stk::mesh::Entity node = bulkData_->get_entity(nodeRank,nodeId);
599  const size_t numElements = bulkData_->num_elements(node);
600  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
601 
602  // extract elements sharing nodes
603  elements.insert(elements.end(), relations, relations + numElements);
604 }
605 
606 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
607  std::vector<int> & relIds) const
608 {
609  // get all relations for node
610  const size_t numElements = bulkData_->num_elements(node);
611  stk::mesh::Entity const* relations = bulkData_->begin_elements(node);
612  stk::mesh::ConnectivityOrdinal const* rel_ids = bulkData_->begin_element_ordinals(node);
613 
614  // extract elements sharing nodes
615  for (size_t i = 0; i < numElements; ++i) {
616  stk::mesh::Entity element = relations[i];
617 
618  // if owned by this processor
619  if(bulkData_->parallel_owner_rank(element) == static_cast<int>(procRank_)) {
620  elements.push_back(element);
621  relIds.push_back(rel_ids[i]);
622  }
623  }
624 }
625 
626 void STK_Interface::getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
627  std::vector<int> & relIds, unsigned int matchType) const
628 {
629  stk::mesh::EntityRank rank;
630  if(matchType == 0)
631  rank = getNodeRank();
632  else if(matchType == 1)
633  rank = getEdgeRank();
634  else if(matchType == 2)
635  rank = getFaceRank();
636  else
637  TEUCHOS_ASSERT(false);
638 
639  stk::mesh::Entity node = bulkData_->get_entity(rank,nodeId);
640 
641  getOwnedElementsSharingNode(node,elements,relIds);
642 }
643 
644 void STK_Interface::getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeIds,std::vector<stk::mesh::Entity> & elements) const
645 {
646  std::vector<stk::mesh::Entity> current;
647 
648  getElementsSharingNode(nodeIds[0],current); // fill it with elements touching first node
649  std::sort(current.begin(),current.end()); // sort for intersection on the pointer
650 
651  // find intersection with remaining nodes
652  for(std::size_t n=1;n<nodeIds.size();++n) {
653  // get elements associated with next node
654  std::vector<stk::mesh::Entity> nextNode;
655  getElementsSharingNode(nodeIds[n],nextNode); // fill it with elements touching first node
656  std::sort(nextNode.begin(),nextNode.end()); // sort for intersection on the pointer ID
657 
658  // intersect next node elements with current element list
659  std::vector<stk::mesh::Entity> intersection(std::min(nextNode.size(),current.size()));
660  std::vector<stk::mesh::Entity>::const_iterator endItr
661  = std::set_intersection(current.begin(),current.end(),
662  nextNode.begin(),nextNode.end(),
663  intersection.begin());
664  std::size_t newLength = endItr-intersection.begin();
665  intersection.resize(newLength);
666 
667  // store intersection
668  current.clear();
669  current = intersection;
670  }
671 
672  // return the elements computed
673  elements = current;
674 }
675 
676 void STK_Interface::getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const
677 {
678  stk::mesh::Entity const* nodeRel = getBulkData()->begin_nodes(element);
679  const size_t numNodes = getBulkData()->num_nodes(element);
680 
681  nodeIds.reserve(numNodes);
682  for(size_t i = 0; i < numNodes; ++i) {
683  nodeIds.push_back(elementGlobalId(nodeRel[i]));
684  }
685 }
686 
688 {
689  entityCounts_.clear();
690  stk::mesh::comm_mesh_counts(*bulkData_,entityCounts_);
691 }
692 
694 {
695  // developed to mirror "comm_mesh_counts" in stk_mesh/base/Comm.cpp
696 
697  const unsigned entityRankCount = metaData_->entity_rank_count();
698  const size_t commCount = 10; // entityRankCount
699 
700  TEUCHOS_ASSERT(entityRankCount<10);
701 
702  // stk::ParallelMachine mach = bulkData_->parallel();
703  stk::ParallelMachine mach = *mpiComm_->getRawMpiComm();
704 
705  std::vector<stk::mesh::EntityId> local(commCount,0);
706 
707  // determine maximum ID for this processor for each entity type
708  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
709  for(stk::mesh::EntityRank i=stk::topology::NODE_RANK; i<entityRankCount; ++i) {
710  std::vector<stk::mesh::Entity> entities;
711 
712  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(i),entities);
713 
714  // determine maximum ID for this processor
715  std::vector<stk::mesh::Entity>::const_iterator itr;
716  for(itr=entities.begin();itr!=entities.end();++itr) {
717  stk::mesh::EntityId id = bulkData_->identifier(*itr);
718  if(id>local[i])
719  local[i] = id;
720  }
721  }
722 
723  // get largest IDs across processors
724  stk::all_reduce(mach,stk::ReduceMax<10>(&local[0]));
725  maxEntityId_.assign(local.begin(),local.begin()+entityRankCount+1);
726 }
727 
728 std::size_t STK_Interface::getEntityCounts(unsigned entityRank) const
729 {
730  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=entityCounts_.size(),std::logic_error,
731  "STK_Interface::getEntityCounts: Entity counts do not include rank: " << entityRank);
732 
733  return entityCounts_[entityRank];
734 }
735 
736 stk::mesh::EntityId STK_Interface::getMaxEntityId(unsigned entityRank) const
737 {
738  TEUCHOS_TEST_FOR_EXCEPTION(entityRank>=maxEntityId_.size(),std::logic_error,
739  "STK_Interface::getMaxEntityId: Max entity ids do not include rank: " << entityRank);
740 
741  return maxEntityId_[entityRank];
742 }
743 
745 {
746  stk::mesh::PartVector emptyPartVector;
747  stk::mesh::create_adjacent_entities(*bulkData_,emptyPartVector);
748 
751 
752  addEdges();
753  if (dimension_ > 2)
754  addFaces();
755 }
756 
757 const double * STK_Interface::getNodeCoordinates(stk::mesh::EntityId nodeId) const
758 {
759  stk::mesh::Entity node = bulkData_->get_entity(getNodeRank(),nodeId);
760  return stk::mesh::field_data(*coordinatesField_,node);
761 }
762 
763 const double * STK_Interface::getNodeCoordinates(stk::mesh::Entity node) const
764 {
765  return stk::mesh::field_data(*coordinatesField_,node);
766 }
767 
768 void STK_Interface::getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
769  std::vector<stk::mesh::EntityId> & subcellIds) const
770 {
771  stk::mesh::EntityRank elementRank = getElementRank();
772  stk::mesh::Entity cell = bulkData_->get_entity(elementRank,elementId);
773 
774  TEUCHOS_TEST_FOR_EXCEPTION(!bulkData_->is_valid(cell),std::logic_error,
775  "STK_Interface::getSubcellIndices: could not find element requested (GID = " << elementId << ")");
776 
777  const size_t numSubcells = bulkData_->num_connectivity(cell, static_cast<stk::mesh::EntityRank>(entityRank));
778  stk::mesh::Entity const* subcells = bulkData_->begin(cell, static_cast<stk::mesh::EntityRank>(entityRank));
779  subcellIds.clear();
780  subcellIds.resize(numSubcells,0);
781 
782  // loop over relations and fill subcell vector
783  for(size_t i = 0; i < numSubcells; ++i) {
784  stk::mesh::Entity subcell = subcells[i];
785  subcellIds[i] = bulkData_->identifier(subcell);
786  }
787 }
788 
789 void STK_Interface::getMyElements(std::vector<stk::mesh::Entity> & elements) const
790 {
791  // setup local ownership
792  stk::mesh::Selector ownedPart = metaData_->locally_owned_part();
793 
794  // grab elements
795  stk::mesh::EntityRank elementRank = getElementRank();
796  stk::mesh::get_selected_entities(ownedPart,bulkData_->buckets(elementRank),elements);
797 }
798 
799 void STK_Interface::getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
800 {
801  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
802 
803  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
804 
805  // setup local ownership
806  // stk::mesh::Selector block = *elementBlock;
807  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & (*elementBlock);
808 
809  // grab elements
810  stk::mesh::EntityRank elementRank = getElementRank();
811  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(elementRank),elements);
812 }
813 
814 void STK_Interface::getNeighborElements(std::vector<stk::mesh::Entity> & elements) const
815 {
816  // setup local ownership
817  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part());
818 
819  // grab elements
820  stk::mesh::EntityRank elementRank = getElementRank();
821  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
822 }
823 
824 void STK_Interface::getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const
825 {
826  stk::mesh::Part * elementBlock = getElementBlockPart(blockID);
827 
828  TEUCHOS_TEST_FOR_EXCEPTION(elementBlock==0,std::logic_error,"Could not find element block \"" << blockID << "\"");
829 
830  // setup local ownership
831  stk::mesh::Selector neighborBlock = (!metaData_->locally_owned_part()) & (*elementBlock);
832 
833  // grab elements
834  stk::mesh::EntityRank elementRank = getElementRank();
835  stk::mesh::get_selected_entities(neighborBlock,bulkData_->buckets(elementRank),elements);
836 }
837 
838 void STK_Interface::getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
839 {
840  stk::mesh::Part * sidePart = getSideset(sideName);
841  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
842  "Unknown side set \"" << sideName << "\"");
843 
844  stk::mesh::Selector side = *sidePart;
845  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & side;
846 
847  // grab elements
848  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
849 }
850 
851 void STK_Interface::getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
852 {
853  stk::mesh::Part * sidePart = getSideset(sideName);
854  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
856  "Unknown side set \"" << sideName << "\"");
858  "Unknown element block \"" << blockName << "\"");
859 
860  stk::mesh::Selector side = *sidePart;
861  stk::mesh::Selector block = *elmtPart;
862  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & side;
863 
864  // grab elements
865  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getSideRank()),sides);
866 }
867 
868 void STK_Interface::getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const
869 {
870  stk::mesh::Part * sidePart = getSideset(sideName);
871  TEUCHOS_TEST_FOR_EXCEPTION(sidePart==0,std::logic_error,
872  "Unknown side set \"" << sideName << "\"");
873 
874  stk::mesh::Selector side = *sidePart;
875 
876  // grab elements
877  stk::mesh::get_selected_entities(side,bulkData_->buckets(getSideRank()),sides);
878 }
879 
880 void STK_Interface::getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const
881 {
882  stk::mesh::Part * sidePart = getSideset(sideName);
883  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
885  "Unknown side set \"" << sideName << "\"");
887  "Unknown element block \"" << blockName << "\"");
888 
889  stk::mesh::Selector side = *sidePart;
890  stk::mesh::Selector block = *elmtPart;
891  stk::mesh::Selector sideBlock = block & side;
892 
893  // grab elements
894  stk::mesh::get_selected_entities(sideBlock,bulkData_->buckets(getSideRank()),sides);
895 }
896 
897 void STK_Interface::getMyNodes(const std::string & nodesetName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const
898 {
899  stk::mesh::Part * nodePart = getNodeset(nodesetName);
900  stk::mesh::Part * elmtPart = getElementBlockPart(blockName);
902  "Unknown node set \"" << nodesetName << "\"");
904  "Unknown element block \"" << blockName << "\"");
905 
906  stk::mesh::Selector nodeset = *nodePart;
907  stk::mesh::Selector block = *elmtPart;
908  stk::mesh::Selector ownedBlock = metaData_->locally_owned_part() & block & nodeset;
909 
910  // grab elements
911  stk::mesh::get_selected_entities(ownedBlock,bulkData_->buckets(getNodeRank()),nodes);
912 }
913 
914 void STK_Interface::getElementBlockNames(std::vector<std::string> & names) const
915 {
916  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
917 
918  names.clear();
919 
920  // fill vector with automagically ordered string values
921  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr; // Element blocks
922  for(blkItr=elementBlocks_.begin();blkItr!=elementBlocks_.end();++blkItr)
923  names.push_back(blkItr->first);
924 }
925 
926 void STK_Interface::getSidesetNames(std::vector<std::string> & names) const
927 {
928  // TEUCHOS_ASSERT(initialized_); // all blocks must have been added
929 
930  names.clear();
931 
932  // fill vector with automagically ordered string values
933  std::map<std::string, stk::mesh::Part*>::const_iterator sideItr; // Element blocks
934  for(sideItr=sidesets_.begin();sideItr!=sidesets_.end();++sideItr)
935  names.push_back(sideItr->first);
936 }
937 
938 void STK_Interface::getNodesetNames(std::vector<std::string> & names) const
939 {
940  names.clear();
941 
942  // fill vector with automagically ordered string values
943  std::map<std::string, stk::mesh::Part*>::const_iterator nodeItr; // Element blocks
944  for(nodeItr=nodesets_.begin();nodeItr!=nodesets_.end();++nodeItr)
945  names.push_back(nodeItr->first);
946 }
947 
948 std::size_t STK_Interface::elementLocalId(stk::mesh::Entity elmt) const
949 {
950  return elementLocalId(bulkData_->identifier(elmt));
951  // const std::size_t * fieldCoords = stk::mesh::field_data(*localIdField_,*elmt);
952  // return fieldCoords[0];
953 }
954 
955 std::size_t STK_Interface::elementLocalId(stk::mesh::EntityId gid) const
956 {
957  // stk::mesh::EntityRank elementRank = getElementRank();
958  // stk::mesh::Entity elmt = bulkData_->get_entity(elementRank,gid);
959  // TEUCHOS_ASSERT(elmt->owner_rank()==procRank_);
960  // return elementLocalId(elmt);
961  std::unordered_map<stk::mesh::EntityId,std::size_t>::const_iterator itr = localIDHash_.find(gid);
962  TEUCHOS_ASSERT(itr!=localIDHash_.end());
963  return itr->second;
964 }
965 
966 
967 std::string STK_Interface::containingBlockId(stk::mesh::Entity elmt)
968 {
969  std::map<std::string,stk::mesh::Part*>::const_iterator itr;
970  for(itr=elementBlocks_.begin();itr!=elementBlocks_.end();++itr)
971  if(bulkData_->bucket(elmt).member(*itr->second))
972  return itr->first;
973  return "";
974 }
975 
976 stk::mesh::Field<double> * STK_Interface::getSolutionField(const std::string & fieldName,
977  const std::string & blockId) const
978 {
979  // look up field in map
980  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
981  iter = fieldNameToSolution_.find(std::make_pair(fieldName,blockId));
982 
983  // check to make sure field was actually found
984  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToSolution_.end(),std::runtime_error,
985  "Solution field name \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
986 
987  return iter->second;
988 }
989 
990 stk::mesh::Field<double> * STK_Interface::getCellField(const std::string & fieldName,
991  const std::string & blockId) const
992 {
993  // look up field in map
994  std::map<std::pair<std::string,std::string>, SolutionFieldType*>::const_iterator
995  iter = fieldNameToCellField_.find(std::make_pair(fieldName,blockId));
996 
997  // check to make sure field was actually found
998  TEUCHOS_TEST_FOR_EXCEPTION(iter==fieldNameToCellField_.end(),std::runtime_error,
999  "Cell field named \"" << fieldName << "\" in block ID \"" << blockId << "\" was not found");
1000 
1001  return iter->second;
1002 }
1003 
1005 {
1006  using Teuchos::RCP;
1007  using Teuchos::rcp;
1008 
1009  if(orderedElementVector_==Teuchos::null) {
1010  // safe because essentially this is a call to modify a mutable object
1011  const_cast<STK_Interface*>(this)->buildLocalElementIDs();
1012  }
1013 
1015 }
1016 
1017 void STK_Interface::addElementBlock(const std::string & name,const CellTopologyData * ctData)
1018 {
1020 
1021  stk::mesh::Part * block = metaData_->get_part(name);
1022  if(block==0) {
1023  block = &metaData_->declare_part_with_topology(name, stk::mesh::get_topology(shards::CellTopology(ctData)));
1024  }
1025 
1026  // construct cell topology object for this block
1028  = Teuchos::rcp(new shards::CellTopology(ctData));
1029 
1030  // add element block part and cell topology
1031  elementBlocks_.insert(std::make_pair(name,block));
1032  elementBlockCT_.insert(std::make_pair(name,ct));
1033 }
1034 
1036 {
1037  dimension_ = metaData_->spatial_dimension();
1038 
1039  // declare coordinates and node parts
1040  coordinatesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::NODE_RANK, coordsString);
1041  edgesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::EDGE_RANK, edgesString);
1042  if (dimension_ > 2)
1043  facesField_ = &metaData_->declare_field<VectorFieldType>(stk::topology::FACE_RANK, facesString);
1044  processorIdField_ = &metaData_->declare_field<ProcIdFieldType>(stk::topology::ELEMENT_RANK, "PROC_ID");
1045  loadBalField_ = &metaData_->declare_field<SolutionFieldType>(stk::topology::ELEMENT_RANK, "LOAD_BAL");
1046 
1047  // stk::mesh::put_field( *coordinatesField_ , getNodeRank() , metaData_->universal_part() );
1048 
1049  nodesPart_ = &metaData_->declare_part(nodesString,getNodeRank());
1050  nodesPartVec_.push_back(nodesPart_);
1051  edgesPart_ = &metaData_->declare_part(edgesString,getEdgeRank());
1052  edgesPartVec_.push_back(edgesPart_);
1053  if (dimension_ > 2) {
1054  facesPart_ = &metaData_->declare_part(facesString,getFaceRank());
1055  facesPartVec_.push_back(facesPart_);
1056  }
1057 }
1058 
1060 {
1061  currentLocalId_ = 0;
1062 
1063  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
1064 
1065  // might be better (faster) to do this by buckets
1066  std::vector<stk::mesh::Entity> elements;
1067  getMyElements(elements);
1068 
1069  for(std::size_t index=0;index<elements.size();++index) {
1070  stk::mesh::Entity element = elements[index];
1071 
1072  // set processor rank
1073  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1074  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1075 
1076  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1077 
1078  currentLocalId_++;
1079  }
1080 
1081  // copy elements into the ordered element vector
1082  orderedElementVector_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>(elements));
1083 
1084  elements.clear();
1085  getNeighborElements(elements);
1086 
1087  for(std::size_t index=0;index<elements.size();++index) {
1088  stk::mesh::Entity element = elements[index];
1089 
1090  // set processor rank
1091  ProcIdData * procId = stk::mesh::field_data(*processorIdField_,element);
1092  procId[0] = Teuchos::as<ProcIdData>(procRank_);
1093 
1094  localIDHash_[bulkData_->identifier(element)] = currentLocalId_;
1095 
1096  currentLocalId_++;
1097  }
1098 
1099  orderedElementVector_->insert(orderedElementVector_->end(),elements.begin(),elements.end());
1100 }
1101 
1103 {
1104  std::vector<std::string> names;
1105  getElementBlockNames(names);
1106 
1107  for(std::size_t b=0;b<names.size();b++) {
1108  // find user specified block weight, otherwise use 1.0
1109  std::map<std::string,double>::const_iterator bw_itr = blockWeights_.find(names[b]);
1110  double blockWeight = (bw_itr!=blockWeights_.end()) ? bw_itr->second : 1.0;
1111 
1112  std::vector<stk::mesh::Entity> elements;
1113  getMyElements(names[b],elements);
1114 
1115  for(std::size_t index=0;index<elements.size();++index) {
1116  // set local element ID
1117  double * loadBal = stk::mesh::field_data(*loadBalField_,elements[index]);
1118  loadBal[0] = blockWeight;
1119  }
1120  }
1121 }
1122 
1123 bool
1124 STK_Interface::isMeshCoordField(const std::string & eBlock,
1125  const std::string & fieldName,
1126  int & axis) const
1127 {
1128  std::map<std::string,std::vector<std::string> >::const_iterator blkItr = meshCoordFields_.find(eBlock);
1129  if(blkItr==meshCoordFields_.end()) {
1130  return false;
1131  }
1132 
1133  axis = 0;
1134  for(axis=0;axis<Teuchos::as<int>(blkItr->second.size());axis++) {
1135  if(blkItr->second[axis]==fieldName)
1136  break; // found axis, break
1137  }
1138 
1139  if(axis>=Teuchos::as<int>(blkItr->second.size()))
1140  return false;
1141 
1142  return true;
1143 }
1144 
1145 std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
1147 {
1149  Teuchos::RCP<std::vector<unsigned int > > type_vec = rcp(new std::vector<unsigned int>);
1150  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & matchers = getPeriodicBCVector();
1151 
1152  // build up the vectors by looping over the matched pair
1153  for(std::size_t m=0;m<matchers.size();m++){
1154  vec = matchers[m]->getMatchedPair(*this,vec);
1155  unsigned int type;
1156  if(matchers[m]->getType() == "coord")
1157  type = 0;
1158  else if(matchers[m]->getType() == "edge")
1159  type = 1;
1160  else if(matchers[m]->getType() == "face")
1161  type = 2;
1162  else
1163  TEUCHOS_ASSERT(false);
1164  type_vec->insert(type_vec->begin(),vec->size()-type_vec->size(),type);
1165  }
1166 
1167  return std::make_pair(vec,type_vec);
1168 }
1169 
1170 bool STK_Interface::validBlockId(const std::string & blockId) const
1171 {
1172  std::map<std::string, stk::mesh::Part*>::const_iterator blkItr = elementBlocks_.find(blockId);
1173 
1174  return blkItr!=elementBlocks_.end();
1175 }
1176 
1177 void STK_Interface::print(std::ostream & os) const
1178 {
1179  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1180 
1181  getElementBlockNames(blockNames);
1182  getSidesetNames(sidesetNames);
1183  getNodesetNames(nodesetNames);
1184 
1185  os << "STK Mesh data:\n";
1186  os << " Spatial dim = " << getDimension() << "\n";
1187  if(getDimension()==2)
1188  os << " Entity counts (Nodes, Edges, Cells) = ( "
1189  << getEntityCounts(getNodeRank()) << ", "
1190  << getEntityCounts(getEdgeRank()) << ", "
1191  << getEntityCounts(getElementRank()) << " )\n";
1192  else if(getDimension()==3)
1193  os << " Entity counts (Nodes, Edges, Faces, Cells) = ( "
1194  << getEntityCounts(getNodeRank()) << ", "
1195  << getEntityCounts(getEdgeRank()) << ", "
1196  << getEntityCounts(getSideRank()) << ", "
1197  << getEntityCounts(getElementRank()) << " )\n";
1198  else
1199  os << " Entity counts (Nodes, Cells) = ( "
1200  << getEntityCounts(getNodeRank()) << ", "
1201  << getEntityCounts(getElementRank()) << " )\n";
1202 
1203  os << " Element blocks = ";
1204  for(std::size_t i=0;i<blockNames.size();i++)
1205  os << "\"" << blockNames[i] << "\" ";
1206  os << "\n";
1207  os << " Sidesets = ";
1208  for(std::size_t i=0;i<sidesetNames.size();i++)
1209  os << "\"" << sidesetNames[i] << "\" ";
1210  os << "\n";
1211  os << " Nodesets = ";
1212  for(std::size_t i=0;i<nodesetNames.size();i++)
1213  os << "\"" << nodesetNames[i] << "\" ";
1214  os << std::endl;
1215 
1216  // print out periodic boundary conditions
1217  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1218  = getPeriodicBCVector();
1219  if(bcVector.size()!=0) {
1220  os << " Periodic BCs:\n";
1221  for(std::size_t i=0;i<bcVector.size();i++)
1222  os << " " << bcVector[i]->getString() << "\n";
1223  os << std::endl;
1224  }
1225 }
1226 
1227 void STK_Interface::printMetaData(std::ostream & os) const
1228 {
1229  std::vector<std::string> blockNames, sidesetNames, nodesetNames;
1230 
1231  getElementBlockNames(blockNames);
1232  getSidesetNames(sidesetNames);
1233  getNodesetNames(nodesetNames);
1234 
1235  os << "STK Meta data:\n";
1236  os << " Element blocks = ";
1237  for(std::size_t i=0;i<blockNames.size();i++)
1238  os << "\"" << blockNames[i] << "\" ";
1239  os << "\n";
1240  os << " Sidesets = ";
1241  for(std::size_t i=0;i<sidesetNames.size();i++)
1242  os << "\"" << sidesetNames[i] << "\" ";
1243  os << "\n";
1244  os << " Nodesets = ";
1245  for(std::size_t i=0;i<nodesetNames.size();i++)
1246  os << "\"" << nodesetNames[i] << "\" ";
1247  os << std::endl;
1248 
1249  // print out periodic boundary conditions
1250  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bcVector
1251  = getPeriodicBCVector();
1252  if(bcVector.size()!=0) {
1253  os << " Periodic BCs:\n";
1254  for(std::size_t i=0;i<bcVector.size();i++)
1255  os << " " << bcVector[i]->getString() << "\n";
1256  os << std::endl;
1257  }
1258 
1259  // print all available fields in meta data
1260  os << " Fields = ";
1261  const stk::mesh::FieldVector & fv = metaData_->get_fields();
1262  for(std::size_t i=0;i<fv.size();i++)
1263  os << "\"" << fv[i]->name() << "\" ";
1264  os << std::endl;
1265 }
1266 
1268 {
1269  std::map<std::string, Teuchos::RCP<shards::CellTopology> >::const_iterator itr;
1270  itr = elementBlockCT_.find(eBlock);
1271 
1272  if(itr==elementBlockCT_.end()) {
1273  std::stringstream ss;
1274  printMetaData(ss);
1275  TEUCHOS_TEST_FOR_EXCEPTION(itr==elementBlockCT_.end(),std::logic_error,
1276  "STK_Interface::getCellTopology: No such element block \"" +eBlock +"\" available.\n\n"
1277  << "STK Meta Data follows: \n" << ss.str());
1278  }
1279 
1280  return itr->second;
1281 }
1282 
1284 {
1285  MPI_Comm newComm;
1286  const int err = MPI_Comm_dup (parallelMach, &newComm);
1287  TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error,
1288  "panzer::STK_Interface: MPI_Comm_dup failed with error \""
1289  << Teuchos::mpiErrorCodeToString (err) << "\".");
1290 
1291  return Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper (newComm,MPI_Comm_free)));
1292 }
1293 
1295 {
1296  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Rebalance not currently supported");
1297 #if 0
1298  // make sure weights have been set (a local operation)
1300 
1301  stk::mesh::Selector selector(getMetaData()->universal_part());
1302  stk::mesh::Selector owned_selector(getMetaData()->locally_owned_part());
1303 
1304  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
1305  out.setOutputToRootOnly(0);
1306  out << "Load balance before: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
1307 
1308  // perform reblance
1309  Teuchos::ParameterList graph;
1310  if(params.begin()!=params.end())
1311  graph.sublist(stk::rebalance::Zoltan::default_parameters_name()) = params;
1312  stk::rebalance::Zoltan zoltan_partition(*mpiComm_->getRawMpiComm(), getDimension(), graph);
1313  stk::rebalance::rebalance(*getBulkData(), owned_selector, &getCoordinatesField(), loadBalField_, zoltan_partition);
1314 
1315  out << "Load balance after: " << stk::rebalance::check_balance(*getBulkData(), loadBalField_, getElementRank(), &selector) << std::endl;
1316 
1317  currentLocalId_ = 0;
1318  orderedElementVector_ = Teuchos::null; // forces rebuild of ordered lists
1319 #endif
1320 }
1321 
1322 } // 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_
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
stk::mesh::EntityRank getFaceRank() const
std::string containingBlockId(stk::mesh::Entity elmt)
ConstIterator begin() const
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_
#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)
void getElementBlockNames(std::vector< std::string > &names) const
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
ConstIterator end() 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)
void setupTransientExodusFile(const std::string &filename)
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)
void printMetaData(std::ostream &os) const
PHX::MDField< const ScalarT, Cell, IP > field
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_
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)
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
Teuchos::RCP< const Teuchos::Comm< int > > comm
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
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