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