Panzer  Version of the Day
Panzer_STK_Interface.hpp
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 #ifndef __Panzer_STK_Interface_hpp__
44 #define __Panzer_STK_Interface_hpp__
45 
46 #include <Teuchos_RCP.hpp>
47 #include <Teuchos_DefaultMpiComm.hpp>
48 
49 #include <stk_mesh/base/Types.hpp>
50 #include <stk_mesh/base/MetaData.hpp>
51 #include <stk_mesh/base/BulkData.hpp>
52 #include <stk_mesh/base/Field.hpp>
53 #include <stk_mesh/base/FieldBase.hpp>
54 #include <stk_mesh/base/MetaData.hpp>
55 #include <stk_mesh/base/CoordinateSystems.hpp>
56 
57 #include "Kokkos_Core.hpp"
58 
59 #include <Shards_CellTopology.hpp>
60 #include <Shards_CellTopologyData.h>
61 
62 #include <PanzerAdaptersSTK_config.hpp>
63 #include <Kokkos_ViewFactory.hpp>
64 
65 #include <unordered_map>
66 
67 #ifdef PANZER_HAVE_IOSS
68 #include <stk_io/StkMeshIoBroker.hpp>
69 #endif
70 
71 #ifdef PANZER_HAVE_PERCEPT
72 namespace percept {
73  class PerceptMesh;
74  class URP_Heterogeneous_3D;
75 }
76 #endif
77 
78 namespace panzer_stk {
79 
80 class PeriodicBC_MatcherBase;
81 
86 public:
87  ElementDescriptor(stk::mesh::EntityId gid,const std::vector<stk::mesh::EntityId> & nodes);
88  virtual ~ElementDescriptor();
89 
90  stk::mesh::EntityId getGID() const { return gid_; }
91  const std::vector<stk::mesh::EntityId> & getNodes() const { return nodes_; }
92 protected:
93  stk::mesh::EntityId gid_;
94  std::vector<stk::mesh::EntityId> nodes_;
95 
97 };
98 
101 Teuchos::RCP<ElementDescriptor>
102 buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes);
103 
105 public:
106  typedef double ProcIdData; // ECC: Not sure why?
107  typedef stk::mesh::Field<double> SolutionFieldType;
108  typedef stk::mesh::Field<double,stk::mesh::Cartesian> VectorFieldType;
109  typedef stk::mesh::Field<ProcIdData> ProcIdFieldType;
110 
111  // some simple exception classes
112  struct ElementBlockException : public std::logic_error
113  { ElementBlockException(const std::string & what) : std::logic_error(what) {} };
114 
115  struct SidesetException : public std::logic_error
116  { SidesetException(const std::string & what) : std::logic_error(what) {} };
117 
118  struct EdgeBlockException : public std::logic_error
119  { EdgeBlockException(const std::string & what) : std::logic_error(what) {} };
120 
121  struct FaceBlockException : public std::logic_error
122  { FaceBlockException(const std::string & what) : std::logic_error(what) {} };
123 
124  STK_Interface();
125 
128  STK_Interface(unsigned dim);
129 
130  STK_Interface(Teuchos::RCP<stk::mesh::MetaData> metaData);
131 
132  // functions called before initialize
134 
137  void addElementBlock(const std::string & name,const CellTopologyData * ctData);
138 
141  void addEdgeBlock(const std::string & name,const CellTopologyData * ctData);
142 
145  void addFaceBlock(const std::string & name,const CellTopologyData * ctData);
146 
149  void addSideset(const std::string & name,const CellTopologyData * ctData);
150 
153  void addNodeset(const std::string & name);
154 
157  void addSolutionField(const std::string & fieldName,const std::string & blockId);
158 
161  void addCellField(const std::string & fieldName,const std::string & blockId);
162 
165  void addEdgeField(const std::string & fieldName,const std::string & blockId);
166 
169  void addFaceField(const std::string & fieldName,const std::string & blockId);
170 
179  void addMeshCoordFields(const std::string & blockId,
180  const std::vector<std::string> & coordField,
181  const std::string & dispPrefix);
182 
187  void addInformationRecords(const std::vector<std::string> & info_records);
188 
190 
202  void initialize(stk::ParallelMachine parallelMach,bool setupIO=true,
203  const bool buildRefinementSupport = false);
204 
210  void instantiateBulkData(stk::ParallelMachine parallelMach);
211 
212  // functions to manage and manipulate bulk data
214 
217  void beginModification();
218 
221  void endModification();
222 
228  void addNode(stk::mesh::EntityId gid, const std::vector<double> & coord);
229 
230  void addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block);
231 
232  void addEdges();
233 
234  void addFaces();
235 
238  void addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset);
239 
242  void addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset);
243 
246  void addEntityToEdgeBlock(stk::mesh::Entity entity,stk::mesh::Part * edgeblock);
249  void addEntitiesToEdgeBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * edgeblock);
250 
253  void addEntityToFaceBlock(stk::mesh::Entity entity,stk::mesh::Part * faceblock);
256  void addEntitiesToFaceBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * faceblock);
257 
258  // Methods to interrogate the mesh topology and structure
260 
264  { return *coordinatesField_; }
265 
269  { return *edgesField_; }
270 
272  { return *facesField_; }
273 
276  const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const;
277 
280  const double * getNodeCoordinates(stk::mesh::Entity node) const;
281 
284  void getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
285  std::vector<stk::mesh::EntityId> & subcellIds) const;
286 
289  void getMyElements(std::vector<stk::mesh::Entity> & elements) const;
290 
293  void getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
294 
298  void getNeighborElements(std::vector<stk::mesh::Entity> & elements) const;
299 
302  void getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
303 
306  void getMyEdges(std::vector<stk::mesh::Entity> & edges) const;
307 
315  void getMyEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const;
316 
325  void getMyEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const;
326 
334  void getAllEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const;
335 
344  void getAllEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const;
345 
348  void getMyFaces(std::vector<stk::mesh::Entity> & faces) const;
349 
357  void getMyFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const;
358 
367  void getMyFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const;
368 
376  void getAllFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const;
377 
386  void getAllFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const;
387 
395  void getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
396 
405  void getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
406 
414  void getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
415 
424  void getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
425 
435  void getMyNodes(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const;
436 
445  stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const;
446 
447  // Utility functions
449 
458  void
459  writeToExodus(const std::string& filename,
460  const bool append = false);
461 
479  void
480  setupExodusFile(const std::string& filename,
481  const bool append = false,
482  const bool append_after_restart_time = false,
483  const double restart_time = 0.0);
484 
497  void
499  double timestep);
500 
520  void
522  const std::string& key,
523  const int& value);
524 
544  void
546  const std::string& key,
547  const double& value);
548 
568  void
570  const std::string& key,
571  const std::vector<int>& value);
572 
592  void
594  const std::string& key,
595  const std::vector<double>& value);
596 
597  // Accessor functions
599 
601  Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;
602 
603  Teuchos::RCP<stk::mesh::BulkData> getBulkData() const { return bulkData_; }
604  Teuchos::RCP<stk::mesh::MetaData> getMetaData() const { return metaData_; }
605 
606 #ifdef PANZER_HAVE_PERCEPT
607  Teuchos::RCP<percept::PerceptMesh> getRefinedMesh() const
609  { TEUCHOS_ASSERT(Teuchos::nonnull(refinedMesh_)); return refinedMesh_; }
610 #endif
611 
612  bool isWritable() const;
613 
614  bool isModifiable() const
615  { if(bulkData_==Teuchos::null) return false;
616  return bulkData_->in_modifiable_state(); }
617 
619  unsigned getDimension() const
620  { return dimension_; }
621 
623  std::size_t getNumElementBlocks() const
624  { return elementBlocks_.size(); }
625 
633  void getElementBlockNames(std::vector<std::string> & names) const;
634 
642  void getSidesetNames(std::vector<std::string> & name) const;
643 
651  void getNodesetNames(std::vector<std::string> & name) const;
652 
660  void getEdgeBlockNames(std::vector<std::string> & names) const;
661 
669  void getFaceBlockNames(std::vector<std::string> & names) const;
670 
672  stk::mesh::Part * getOwnedPart() const
673  { return &getMetaData()->locally_owned_part(); } // I don't like the pointer access here, but it will do for now!
674 
676  stk::mesh::Part * getElementBlockPart(const std::string & name) const
677  {
678  std::map<std::string, stk::mesh::Part*>::const_iterator itr = elementBlocks_.find(name); // Element blocks
679  if(itr==elementBlocks_.end()) return 0;
680  return itr->second;
681  }
682 
684  stk::mesh::Part * getEdgeBlock(const std::string & name) const
685  {
686  std::map<std::string, stk::mesh::Part*>::const_iterator itr = edgeBlocks_.find(name); // edge blocks
687  if(itr==edgeBlocks_.end()) return 0;
688  return itr->second;
689  }
690 
692  stk::mesh::Part * getFaceBlock(const std::string & name) const
693  {
694  std::map<std::string, stk::mesh::Part*>::const_iterator itr = faceBlocks_.find(name); // face blocks
695  if(itr==faceBlocks_.end()) return 0;
696  return itr->second;
697  }
698 
700  std::size_t getNumSidesets() const
701  { return sidesets_.size(); }
702 
703  stk::mesh::Part * getSideset(const std::string & name) const
704  {
705  auto itr = sidesets_.find(name);
706  return (itr != sidesets_.end()) ? itr->second : nullptr;
707  }
708 
710  std::size_t getNumNodesets() const
711  { return nodesets_.size(); }
712 
713  stk::mesh::Part * getNodeset(const std::string & name) const
714  {
715  auto itr = nodesets_.find(name);
716  return (itr != nodesets_.end()) ? itr->second : nullptr;
717  }
718 
720  std::size_t getEntityCounts(unsigned entityRank) const;
721 
723  stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const;
724 
725  // Utilities
727 
729  void getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const;
730 
732  void getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const;
733 
736  void getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
737  std::vector<int> & relIds) const;
738 
741  void getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
742  std::vector<int> & relIds, unsigned int matchType) const;
743 
744 
746  void getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeId,std::vector<stk::mesh::Entity> & elements) const;
747 
749  void buildSubcells();
750 
753  std::size_t elementLocalId(stk::mesh::Entity elmt) const;
754 
757  std::size_t elementLocalId(stk::mesh::EntityId gid) const;
758 
761  inline stk::mesh::EntityId elementGlobalId(std::size_t lid) const
762  { return bulkData_->identifier((*orderedElementVector_)[lid]); }
763 
766  inline stk::mesh::EntityId elementGlobalId(stk::mesh::Entity elmt) const
767  { return bulkData_->identifier(elmt); }
768 
771  bool isEdgeLocal(stk::mesh::Entity edge) const;
772 
775  bool isEdgeLocal(stk::mesh::EntityId gid) const;
776 
779  std::size_t edgeLocalId(stk::mesh::Entity elmt) const;
780 
783  std::size_t edgeLocalId(stk::mesh::EntityId gid) const;
784 
787  inline stk::mesh::EntityId edgeGlobalId(std::size_t lid) const
788  { return bulkData_->identifier((*orderedEdgeVector_)[lid]); }
789 
792  inline stk::mesh::EntityId edgeGlobalId(stk::mesh::Entity edge) const
793  { return bulkData_->identifier(edge); }
794 
797  bool isFaceLocal(stk::mesh::Entity face) const;
798 
801  bool isFaceLocal(stk::mesh::EntityId gid) const;
802 
805  std::size_t faceLocalId(stk::mesh::Entity elmt) const;
806 
809  std::size_t faceLocalId(stk::mesh::EntityId gid) const;
810 
813  inline stk::mesh::EntityId faceGlobalId(std::size_t lid) const
814  { return bulkData_->identifier((*orderedFaceVector_)[lid]); }
815 
818  inline stk::mesh::EntityId faceGlobalId(stk::mesh::Entity face) const
819  { return bulkData_->identifier(face); }
820 
823  inline unsigned entityOwnerRank(stk::mesh::Entity entity) const
824  { return bulkData_->parallel_owner_rank(entity); }
825 
828  inline bool isValid(stk::mesh::Entity entity) const
829  { return bulkData_->is_valid(entity); }
830 
833  std::string containingBlockId(stk::mesh::Entity elmt) const;
834 
839  stk::mesh::Field<double> * getSolutionField(const std::string & fieldName,
840  const std::string & blockId) const;
841 
846  stk::mesh::Field<double> * getCellField(const std::string & fieldName,
847  const std::string & blockId) const;
848 
853  stk::mesh::Field<double> * getEdgeField(const std::string & fieldName,
854  const std::string & blockId) const;
855 
860  stk::mesh::Field<double> * getFaceField(const std::string & fieldName,
861  const std::string & blockId) const;
862 
864 
866  bool isInitialized() const { return initialized_; }
867 
871  Teuchos::RCP<const std::vector<stk::mesh::Entity> > getElementsOrderedByLID() const;
872 
887  template <typename ArrayT>
888  void setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
889  const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
890 
905  template <typename ArrayT>
906  void getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
907  const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const;
908 
923  template <typename ArrayT>
924  void setCellFieldData(const std::string & fieldName,const std::string & blockId,
925  const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
926 
930  Teuchos::RCP<const std::vector<stk::mesh::Entity> > getEdgesOrderedByLID() const;
931 
935  Teuchos::RCP<const std::vector<stk::mesh::Entity> > getFacesOrderedByLID() const;
936 
951  template <typename ArrayT>
952  void setEdgeFieldData(const std::string & fieldName,const std::string & blockId,
953  const std::vector<std::size_t> & localEdgeIds,const ArrayT & edgeValues,double scaleValue=1.0);
954 
969  template <typename ArrayT>
970  void setFaceFieldData(const std::string & fieldName,const std::string & blockId,
971  const std::vector<std::size_t> & localFaceIds,const ArrayT & faceValues,double scaleValue=1.0);
972 
981  template <typename ArrayT>
982  void getElementVertices(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
983 
992  template <typename ArrayT>
993  void getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
994 
1004  template <typename ArrayT>
1005  void getElementVertices(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
1006 
1016  template <typename ArrayT>
1017  void getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1018 
1027  template <typename ArrayT>
1028  void getElementVerticesNoResize(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
1029 
1038  template <typename ArrayT>
1039  void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1040 
1050  template <typename ArrayT>
1051  void getElementVerticesNoResize(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
1052 
1062  template <typename ArrayT>
1063  void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1064 
1065  // const stk::mesh::FEMInterface & getFEMInterface() const
1066  // { return *femPtr_; }
1067 
1068  stk::mesh::EntityRank getElementRank() const { return stk::topology::ELEMENT_RANK; }
1069  stk::mesh::EntityRank getSideRank() const { return metaData_->side_rank(); }
1070  stk::mesh::EntityRank getFaceRank() const { return stk::topology::FACE_RANK; }
1071  stk::mesh::EntityRank getEdgeRank() const { return stk::topology::EDGE_RANK; }
1072  stk::mesh::EntityRank getNodeRank() const { return stk::topology::NODE_RANK; }
1073 
1076  void initializeFromMetaData();
1077 
1080  void buildLocalElementIDs();
1081 
1084  void buildLocalEdgeIDs();
1085 
1088  void buildLocalFaceIDs();
1089 
1092  const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
1094  { return periodicBCs_; }
1095 
1098  std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
1100  { return periodicBCs_; }
1101 
1107  void addPeriodicBC(const Teuchos::RCP<const PeriodicBC_MatcherBase> & bc)
1108  { periodicBCs_.push_back(bc); }
1109 
1115  void addPeriodicBCs(const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bc_vec)
1116  { periodicBCs_.insert(periodicBCs_.end(),bc_vec.begin(),bc_vec.end()); }
1117 
1118  std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
1119  getPeriodicNodePairing() const;
1120 
1123  bool validBlockId(const std::string & blockId) const;
1124 
1127  void print(std::ostream & os) const;
1128 
1131  void printMetaData(std::ostream & os) const;
1132 
1135  Teuchos::RCP<const shards::CellTopology> getCellTopology(const std::string & eBlock) const;
1136 
1141  double getCurrentStateTime() const { return currentStateTime_; }
1142 
1148  double getInitialStateTime() const { return initialStateTime_; }
1149 
1154  void setInitialStateTime(double value) { initialStateTime_ = value; }
1155 
1158  void rebalance(const Teuchos::ParameterList & params);
1159 
1163  void setBlockWeight(const std::string & blockId,double weight)
1164  { blockWeights_[blockId] = weight; }
1165 
1172  void setUseFieldCoordinates(bool useFieldCoordinates)
1173  { useFieldCoordinates_ = useFieldCoordinates; }
1174 
1177  { return useFieldCoordinates_; }
1178 
1180  void setUseLowerCaseForIO(bool useLowerCase)
1181  { useLowerCase_ = useLowerCase; }
1182 
1185  { return useLowerCase_; }
1186 
1197  template <typename ArrayT>
1198  void getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1199 
1200  template <typename ArrayT>
1201  void getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
1202  const std::string & eBlock, ArrayT & vertices) const;
1203 
1213  template <typename ArrayT>
1214  void getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1215 
1216  template <typename ArrayT>
1217  void getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1218 
1224  void refineMesh(const int numberOfLevels, const bool deleteParentElements);
1225 
1226 public: // static operations
1227  static const std::string coordsString;
1228  static const std::string nodesString;
1229  static const std::string edgesString;
1230  static const std::string edgeBlockString;
1231  static const std::string faceBlockString;
1232  static const std::string facesString;
1233 
1234 protected:
1235 
1238  void buildEntityCounts();
1239 
1242  void buildMaxEntityIds();
1243 
1247  void initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
1248  bool setupIO);
1249 
1254  Teuchos::RCP<Teuchos::MpiComm<int> > getSafeCommunicator(stk::ParallelMachine parallelMach) const;
1255 
1262 
1266  bool isMeshCoordField(const std::string & eBlock,const std::string & fieldName,int & axis) const;
1267 
1283  template <typename ArrayT>
1284  void setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
1285  const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues);
1286 
1287  std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > periodicBCs_;
1288 
1289  Teuchos::RCP<stk::mesh::MetaData> metaData_;
1290  Teuchos::RCP<stk::mesh::BulkData> bulkData_;
1291 #ifdef PANZER_HAVE_PERCEPT
1292  Teuchos::RCP<percept::PerceptMesh> refinedMesh_;
1293  Teuchos::RCP<percept::URP_Heterogeneous_3D> breakPattern_;
1294 #endif
1295 
1296  std::map<std::string, stk::mesh::Part*> elementBlocks_; // Element blocks
1297  std::map<std::string, stk::mesh::Part*> sidesets_; // Side sets
1298  std::map<std::string, stk::mesh::Part*> nodesets_; // Node sets
1299  std::map<std::string, stk::mesh::Part*> edgeBlocks_; // Edge blocks
1300  std::map<std::string, stk::mesh::Part*> faceBlocks_; // Face blocks
1301 
1302  std::map<std::string, Teuchos::RCP<shards::CellTopology> > elementBlockCT_;
1303  std::map<std::string, Teuchos::RCP<shards::CellTopology> > edgeBlockCT_;
1304  std::map<std::string, Teuchos::RCP<shards::CellTopology> > faceBlockCT_;
1305 
1306  // for storing/accessing nodes
1307  stk::mesh::Part * nodesPart_;
1308  std::vector<stk::mesh::Part*> nodesPartVec_;
1309  stk::mesh::Part * edgesPart_;
1310  std::vector<stk::mesh::Part*> edgesPartVec_;
1311  stk::mesh::Part * facesPart_;
1312  std::vector<stk::mesh::Part*> facesPartVec_;
1313 
1319 
1320  // maps field names to solution field stk mesh handles
1321  std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToSolution_;
1322  std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToCellField_;
1323  std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToEdgeField_;
1324  std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToFaceField_;
1325 
1326  // use a set to maintain a list of unique information records
1327  std::set<std::string> informationRecords_;
1328 
1329  unsigned dimension_;
1330 
1332 
1333  // how many elements, faces, edges, and nodes are there globally
1334  std::vector<std::size_t> entityCounts_;
1335 
1336  // what is maximum entity ID
1337  std::vector<stk::mesh::EntityId> maxEntityId_;
1338 
1339  unsigned procRank_;
1340  std::size_t currentLocalId_;
1341 
1342  Teuchos::RCP<Teuchos::MpiComm<int> > mpiComm_;
1343 
1344  double initialStateTime_; // the time stamp at the time this object was constructed (default 0.0)
1345  double currentStateTime_; // the time stamp set by the user when writeToExodus is called (default 0.0)
1346 
1347 #ifdef PANZER_HAVE_IOSS
1348  // I/O support
1349  Teuchos::RCP<stk::io::StkMeshIoBroker> meshData_;
1350  int meshIndex_;
1351 
1356  enum class GlobalVariable
1357  {
1358  ADD,
1359  WRITE
1360  }; // end of enum class GlobalVariable
1361 
1378  void
1379  globalToExodus(
1380  const GlobalVariable& flag);
1381 
1385  Teuchos::ParameterList globalData_;
1386 #endif
1387 
1388  // uses lazy evaluation
1389  mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedElementVector_;
1390 
1391  // uses lazy evaluation
1392  mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedEdgeVector_;
1393 
1394  // uses lazy evaluation
1395  mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedFaceVector_;
1396 
1397  // for element block weights
1398  std::map<std::string,double> blockWeights_;
1399 
1400  std::unordered_map<stk::mesh::EntityId,std::size_t> localIDHash_;
1401  std::unordered_map<stk::mesh::EntityId,std::size_t> localEdgeIDHash_;
1402  std::unordered_map<stk::mesh::EntityId,std::size_t> localFaceIDHash_;
1403 
1404  // Store mesh displacement fields by element block. This map
1405  // goes like this meshCoordFields_[eBlock][axis_index] => coordinate FieldName
1406  // goes like this meshDispFields_[eBlock][axis_index] => displacement FieldName
1407  std::map<std::string,std::vector<std::string> > meshCoordFields_; // coordinate fields written by user
1408  std::map<std::string,std::vector<std::string> > meshDispFields_; // displacement fields, output to exodus
1409 
1411 
1413 
1414  // Object describing how to sort a vector of elements using
1415  // local ID as the key, very short lived object
1417  public:
1418  LocalIdCompare(const STK_Interface * mesh) : mesh_(mesh) {}
1419 
1420  // Compares two stk mesh entities based on local ID
1421  bool operator() (stk::mesh::Entity a,stk::mesh::Entity b)
1422  { return mesh_->elementLocalId(a) < mesh_->elementLocalId(b);}
1423 
1424  private:
1426  };
1427 };
1428 
1429 template <typename ArrayT>
1430 void STK_Interface::setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
1431  const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
1432 {
1433  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1434  auto solutionValues_h = Kokkos::create_mirror_view(solutionValues);
1435  Kokkos::deep_copy(solutionValues_h, solutionValues);
1436 
1437  int field_axis = -1;
1438  if(isMeshCoordField(blockId,fieldName,field_axis)) {
1439  setDispFieldData(fieldName,blockId,field_axis,localElementIds,solutionValues_h);
1440  return;
1441  }
1442 
1443  SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1444 
1445  for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1446  std::size_t localId = localElementIds[cell];
1447  stk::mesh::Entity element = elements[localId];
1448 
1449  // loop over nodes set solution values
1450  const size_t num_nodes = bulkData_->num_nodes(element);
1451  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1452  for(std::size_t i=0; i<num_nodes; ++i) {
1453  stk::mesh::Entity node = nodes[i];
1454 
1455  double * solnData = stk::mesh::field_data(*field,node);
1456  // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1457  solnData[0] = scaleValue*solutionValues_h(cell,i);
1458  }
1459  }
1460 }
1461 
1462 template <typename ArrayT>
1463 void STK_Interface::setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
1464  const std::vector<std::size_t> & localElementIds,const ArrayT & dispValues)
1465 {
1466  TEUCHOS_ASSERT(axis>=0); // sanity check
1467 
1468  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1469 
1470  SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1471  const VectorFieldType & coord_field = this->getCoordinatesField();
1472 
1473  for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1474  std::size_t localId = localElementIds[cell];
1475  stk::mesh::Entity element = elements[localId];
1476 
1477  // loop over nodes set solution values
1478  const size_t num_nodes = bulkData_->num_nodes(element);
1479  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1480  for(std::size_t i=0; i<num_nodes; ++i) {
1481  stk::mesh::Entity node = nodes[i];
1482 
1483  double * solnData = stk::mesh::field_data(*field,node);
1484  double * coordData = stk::mesh::field_data(coord_field,node);
1485  // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1486  solnData[0] = dispValues(cell,i)-coordData[axis];
1487  }
1488  }
1489 }
1490 
1491 template <typename ArrayT>
1492 void STK_Interface::getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
1493  const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const
1494 {
1495  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1496 
1497  solutionValues = Kokkos::createDynRankView(solutionValues,
1498  "solutionValues",
1499  localElementIds.size(),
1500  bulkData_->num_nodes(elements[localElementIds[0]]));
1501 
1502  SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1503 
1504  for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1505  std::size_t localId = localElementIds[cell];
1506  stk::mesh::Entity element = elements[localId];
1507 
1508  // loop over nodes set solution values
1509  const size_t num_nodes = bulkData_->num_nodes(element);
1510  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1511  for(std::size_t i=0; i<num_nodes; ++i) {
1512  stk::mesh::Entity node = nodes[i];
1513 
1514  double * solnData = stk::mesh::field_data(*field,node);
1515  // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1516  solutionValues(cell,i) = solnData[0];
1517  }
1518  }
1519 }
1520 
1521 template <typename ArrayT>
1522 void STK_Interface::setCellFieldData(const std::string & fieldName,const std::string & blockId,
1523  const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
1524 {
1525  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1526 
1527  SolutionFieldType * field = this->getCellField(fieldName,blockId);
1528 
1529  auto solutionValues_h = Kokkos::create_mirror_view(solutionValues);
1530  Kokkos::deep_copy(solutionValues_h, solutionValues);
1531 
1532  for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1533  std::size_t localId = localElementIds[cell];
1534  stk::mesh::Entity element = elements[localId];
1535 
1536  double * solnData = stk::mesh::field_data(*field,element);
1537  TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1538  solnData[0] = scaleValue*solutionValues_h.access(cell,0);
1539  }
1540 }
1541 
1542 template <typename ArrayT>
1543 void STK_Interface::setEdgeFieldData(const std::string & fieldName,const std::string & blockId,
1544  const std::vector<std::size_t> & localEdgeIds,const ArrayT & edgeValues,double scaleValue)
1545 {
1546  const std::vector<stk::mesh::Entity> & edges = *(this->getEdgesOrderedByLID());
1547 
1548  SolutionFieldType * field = this->getEdgeField(fieldName,blockId);
1549 
1550  auto edgeValues_h = Kokkos::create_mirror_view(edgeValues);
1551  Kokkos::deep_copy(edgeValues_h, edgeValues);
1552 
1553  for(std::size_t idx=0;idx<localEdgeIds.size();idx++) {
1554  std::size_t localId = localEdgeIds[idx];
1555  stk::mesh::Entity edge = edges[localId];
1556 
1557  double * solnData = stk::mesh::field_data(*field,edge);
1558  TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1559  solnData[0] = scaleValue*edgeValues_h.access(idx,0);
1560  }
1561 }
1562 
1563 template <typename ArrayT>
1564 void STK_Interface::setFaceFieldData(const std::string & fieldName,const std::string & blockId,
1565  const std::vector<std::size_t> & localFaceIds,const ArrayT & faceValues,double scaleValue)
1566 {
1567  const std::vector<stk::mesh::Entity> & faces = *(this->getFacesOrderedByLID());
1568 
1569  SolutionFieldType * field = this->getFaceField(fieldName,blockId);
1570 
1571  auto faceValues_h = Kokkos::create_mirror_view(faceValues);
1572  Kokkos::deep_copy(faceValues_h, faceValues);
1573 
1574  for(std::size_t idx=0;idx<localFaceIds.size();idx++) {
1575  std::size_t localId = localFaceIds[idx];
1576  stk::mesh::Entity face = faces[localId];
1577 
1578  double * solnData = stk::mesh::field_data(*field,face);
1579  TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1580  solnData[0] = scaleValue*faceValues_h.access(idx,0);
1581  }
1582 }
1583 
1584 template <typename ArrayT>
1585 void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
1586 {
1587  if(!useFieldCoordinates_) {
1588  //
1589  // gather from the intrinsic mesh coordinates (non-lagrangian)
1590  //
1591 
1592  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1593 
1594  // convert to a vector of entity objects
1595  std::vector<stk::mesh::Entity> selected_elements;
1596  for(std::size_t cell=0;cell<localElementIds.size();cell++)
1597  selected_elements.push_back(elements[localElementIds[cell]]);
1598 
1599  getElementVertices_FromCoords(selected_elements,vertices);
1600  }
1601  else {
1602  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1603  "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
1604  "without specifying an element block.");
1605  }
1606 }
1607 
1608 template <typename ArrayT>
1609 void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1610 {
1611  if(!useFieldCoordinates_) {
1612  getElementVertices_FromCoords(elements,vertices);
1613  }
1614  else {
1615  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1616  "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
1617  "without specifying an element block.");
1618  }
1619 }
1620 
1621 template <typename ArrayT>
1622 void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1623 {
1624  if(!useFieldCoordinates_) {
1625  getElementVertices_FromCoords(elements,vertices);
1626  }
1627  else {
1628  getElementVertices_FromField(elements,eBlock,vertices);
1629  }
1630 }
1631 
1632 template <typename ArrayT>
1633 void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
1634 {
1635  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1636 
1637  // convert to a vector of entity objects
1638  std::vector<stk::mesh::Entity> selected_elements;
1639  for(std::size_t cell=0;cell<localElementIds.size();cell++)
1640  selected_elements.push_back(elements[localElementIds[cell]]);
1641 
1642  if(!useFieldCoordinates_) {
1643  getElementVertices_FromCoords(selected_elements,vertices);
1644  }
1645  else {
1646  getElementVertices_FromField(selected_elements,eBlock,vertices);
1647  }
1648 }
1649 
1650 template <typename ArrayT>
1651 void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
1652 {
1653  if(!useFieldCoordinates_) {
1654  //
1655  // gather from the intrinsic mesh coordinates (non-lagrangian)
1656  //
1657 
1658  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1659 
1660  // convert to a vector of entity objects
1661  std::vector<stk::mesh::Entity> selected_elements;
1662  for(std::size_t cell=0;cell<localElementIds.size();cell++)
1663  selected_elements.push_back(elements[localElementIds[cell]]);
1664 
1665  getElementVertices_FromCoordsNoResize(selected_elements,vertices);
1666  }
1667  else {
1668  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1669  "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
1670  "without specifying an element block.");
1671  }
1672 }
1673 
1674 template <typename ArrayT>
1675 void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1676 {
1677  if(!useFieldCoordinates_) {
1678  getElementVertices_FromCoordsNoResize(elements,vertices);
1679  }
1680  else {
1681  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1682  "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
1683  "without specifying an element block.");
1684  }
1685 }
1686 
1687 template <typename ArrayT>
1688 void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1689 {
1690  if(!useFieldCoordinates_) {
1691  getElementVertices_FromCoordsNoResize(elements,vertices);
1692  }
1693  else {
1694  getElementVertices_FromFieldNoResize(elements,eBlock,vertices);
1695  }
1696 }
1697 
1698 template <typename ArrayT>
1699 void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
1700 {
1701  const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1702 
1703  // convert to a vector of entity objects
1704  std::vector<stk::mesh::Entity> selected_elements;
1705  for(std::size_t cell=0;cell<localElementIds.size();cell++)
1706  selected_elements.push_back(elements[localElementIds[cell]]);
1707 
1708  if(!useFieldCoordinates_) {
1709  getElementVertices_FromCoordsNoResize(selected_elements,vertices);
1710  }
1711  else {
1712  getElementVertices_FromFieldNoResize(selected_elements,eBlock,vertices);
1713  }
1714 }
1715 
1716 template <typename ArrayT>
1717 void STK_Interface::getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1718 {
1719  // nothing to do! silently return
1720  if(elements.size() == 0) {
1721  vertices = Kokkos::createDynRankView(vertices, "vertices", 0, 0, 0);
1722  return;
1723  }
1724 
1725  //
1726  // gather from the intrinsic mesh coordinates (non-lagrangian)
1727  //
1728 
1729  // get *master* cell toplogy...(belongs to first element)
1730  const auto masterVertexCount
1731  = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1732 
1733  // allocate space
1734  vertices = Kokkos::createDynRankView(vertices, "vertices", elements.size(), masterVertexCount,getDimension());
1735  auto vertices_h = Kokkos::create_mirror_view(vertices);
1736  Kokkos::deep_copy(vertices_h, vertices);
1737 
1738  // loop over each requested element
1739  const auto dim = getDimension();
1740  for(std::size_t cell = 0; cell < elements.size(); cell++) {
1741  const auto element = elements[cell];
1742  TEUCHOS_ASSERT(element != 0);
1743 
1744  const auto vertexCount
1745  = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
1746  TEUCHOS_TEST_FOR_EXCEPTION(vertexCount != masterVertexCount, std::runtime_error,
1747  "In call to STK_Interface::getElementVertices all elements "
1748  "must have the same vertex count!");
1749 
1750  // loop over all element nodes
1751  const size_t num_nodes = bulkData_->num_nodes(element);
1752  auto const* nodes = bulkData_->begin_nodes(element);
1753  TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
1754  "In call to STK_Interface::getElementVertices cardinality of "
1755  "element node relations must be the vertex count!");
1756  for(std::size_t node = 0; node < num_nodes; ++node) {
1757  const double * coord = getNodeCoordinates(nodes[node]);
1758 
1759  // set each dimension of the coordinate
1760  for(unsigned d=0;d<dim;d++)
1761  vertices_h(cell,node,d) = coord[d];
1762  }
1763  }
1764  Kokkos::deep_copy(vertices, vertices_h);
1765 }
1766 
1767 template <typename ArrayT>
1768 void STK_Interface::getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1769 {
1770  // nothing to do! silently return
1771  if(elements.size()==0) {
1772  return;
1773  }
1774 
1775  //
1776  // gather from the intrinsic mesh coordinates (non-lagrangian)
1777  //
1778 
1779  // get *master* cell toplogy...(belongs to first element)
1780  unsigned masterVertexCount
1781  = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1782 
1783  // loop over each requested element
1784  unsigned dim = getDimension();
1785  auto vertices_h = Kokkos::create_mirror_view(vertices);
1786  for(std::size_t cell=0;cell<elements.size();cell++) {
1787  stk::mesh::Entity element = elements[cell];
1788  TEUCHOS_ASSERT(element!=0);
1789 
1790  unsigned vertexCount
1791  = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
1792  TEUCHOS_TEST_FOR_EXCEPTION(vertexCount!=masterVertexCount,std::runtime_error,
1793  "In call to STK_Interface::getElementVertices all elements "
1794  "must have the same vertex count!");
1795 
1796  // loop over all element nodes
1797  const size_t num_nodes = bulkData_->num_nodes(element);
1798  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1799  TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
1800  "In call to STK_Interface::getElementVertices cardinality of "
1801  "element node relations must be the vertex count!");
1802  for(std::size_t node=0; node<num_nodes; ++node) {
1803  const double * coord = getNodeCoordinates(nodes[node]);
1804 
1805  // set each dimension of the coordinate
1806  for(unsigned d=0;d<dim;d++)
1807  vertices_h(cell,node,d) = coord[d];
1808  }
1809  }
1810  Kokkos::deep_copy(vertices, vertices_h);
1811 }
1812 
1813 template <typename ArrayT>
1814 void STK_Interface::getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1815 {
1816  TEUCHOS_ASSERT(useFieldCoordinates_);
1817 
1818  // nothing to do! silently return
1819  if(elements.size()==0) {
1820  vertices = Kokkos::createDynRankView(vertices,"vertices",0,0,0);
1821  return;
1822  }
1823 
1824  // get *master* cell toplogy...(belongs to first element)
1825  unsigned masterVertexCount
1826  = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1827 
1828  // allocate space
1829  vertices = Kokkos::createDynRankView(vertices,"vertices",elements.size(),masterVertexCount,getDimension());
1830  auto vertices_h = Kokkos::create_mirror_view(vertices);
1831  std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
1832  if(itr==meshCoordFields_.end()) {
1833  // no coordinate field set for this element block
1834  TEUCHOS_ASSERT(false);
1835  }
1836 
1837  const std::vector<std::string> & coordField = itr->second;
1838  std::vector<SolutionFieldType*> fields(getDimension());
1839  for(std::size_t d=0;d<fields.size();d++) {
1840  fields[d] = this->getSolutionField(coordField[d],eBlock);
1841  }
1842 
1843  for(std::size_t cell=0;cell<elements.size();cell++) {
1844  stk::mesh::Entity element = elements[cell];
1845 
1846  // loop over nodes set solution values
1847  const size_t num_nodes = bulkData_->num_nodes(element);
1848  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1849  for(std::size_t i=0; i<num_nodes; ++i) {
1850  stk::mesh::Entity node = nodes[i];
1851 
1852  const double * coord = getNodeCoordinates(node);
1853 
1854  for(unsigned d=0;d<getDimension();d++) {
1855  double * solnData = stk::mesh::field_data(*fields[d],node);
1856 
1857  // recall mesh field coordinates are stored as displacements
1858  // from the mesh coordinates, make sure to add them together
1859  vertices_h(cell,i,d) = solnData[0]+coord[d];
1860  }
1861  }
1862  }
1863  Kokkos::deep_copy(vertices, vertices_h);
1864 }
1865 
1866 template <typename ArrayT>
1867 void STK_Interface::getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
1868  const std::string & eBlock, ArrayT & vertices) const
1869 {
1870  TEUCHOS_ASSERT(useFieldCoordinates_);
1871 
1872  // nothing to do! silently return
1873  if(elements.size()==0) {
1874  return;
1875  }
1876 
1877  std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
1878  if(itr==meshCoordFields_.end()) {
1879  // no coordinate field set for this element block
1880  TEUCHOS_ASSERT(false);
1881  }
1882 
1883  const std::vector<std::string> & coordField = itr->second;
1884  std::vector<SolutionFieldType*> fields(getDimension());
1885  for(std::size_t d=0;d<fields.size();d++) {
1886  fields[d] = this->getSolutionField(coordField[d],eBlock);
1887  }
1888 
1889  for(std::size_t cell=0;cell<elements.size();cell++) {
1890  stk::mesh::Entity element = elements[cell];
1891 
1892  // loop over nodes set solution values
1893  const size_t num_nodes = bulkData_->num_nodes(element);
1894  stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1895  for(std::size_t i=0; i<num_nodes; ++i) {
1896  stk::mesh::Entity node = nodes[i];
1897 
1898  const double * coord = getNodeCoordinates(node);
1899 
1900  for(unsigned d=0;d<getDimension();d++) {
1901  double * solnData = stk::mesh::field_data(*fields[d],node);
1902 
1903  // recall mesh field coordinates are stored as displacements
1904  // from the mesh coordinates, make sure to add them together
1905  vertices(cell,i,d) = solnData[0]+coord[d];
1906  }
1907  }
1908  }
1909 }
1910 
1911 }
1912 
1913 #endif
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_
stk::mesh::EntityId faceGlobalId(std::size_t lid) const
bool isMeshCoordField(const std::string &eBlock, const std::string &fieldName, int &axis) const
void addNodeset(const std::string &name)
void getElementVertices(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
Teuchos::RCP< stk::mesh::BulkData > bulkData_
void setEdgeFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localEdgeIds, const ArrayT &edgeValues, double scaleValue=1.0)
void setSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
std::vector< stk::mesh::Part * > facesPartVec_
stk::mesh::EntityId getGID() const
void getElementVertices_FromCoordsNoResize(const std::vector< stk::mesh::Entity > &elements, ArrayT &vertices) const
std::vector< stk::mesh::EntityId > nodes_
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.
void setDispFieldData(const std::string &fieldName, const std::string &blockId, int axis, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues)
void addPeriodicBC(const Teuchos::RCP< const PeriodicBC_MatcherBase > &bc)
std::size_t elementLocalId(stk::mesh::Entity elmt) const
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector()
stk::mesh::Field< double > SolutionFieldType
std::vector< stk::mesh::Part * > edgesPartVec_
const std::vector< stk::mesh::EntityId > & getNodes() const
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_
std::size_t getNumSidesets() const
get the side set count
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
stk::mesh::Part * getOwnedPart() const
Get a pointer to the locally owned part.
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
bool isValid(stk::mesh::Entity entity) 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
const VectorFieldType & getFacesField() const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
std::map< std::string, std::vector< std::string > > meshDispFields_
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack... dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
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.
void setUseFieldCoordinates(bool useFieldCoordinates)
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
void getElementVertices_FromField(const std::vector< stk::mesh::Entity > &elements, const std::string &eBlock, ArrayT &vertices) 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
void setFaceFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localFaceIds, const ArrayT &faceValues, double scaleValue=1.0)
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_
void getElementVertices_FromCoords(const std::vector< stk::mesh::Entity > &elements, ArrayT &vertices) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
stk::mesh::EntityId edgeGlobalId(std::size_t lid) const
void getElementVertices_FromFieldNoResize(const std::vector< stk::mesh::Entity > &elements, const std::string &eBlock, ArrayT &vertices) const
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?
stk::mesh::EntityId edgeGlobalId(stk::mesh::Entity edge) const
std::vector< stk::mesh::Part * > nodesPartVec_
static const std::string edgeBlockString
static const std::string edgesString
std::size_t getNumElementBlocks() const
get the block count
void printMetaData(std::ostream &os) const
void initializeFieldsInSTK(const std::map< std::pair< std::string, std::string >, SolutionFieldType *> &nameToField, bool setupIO)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCs_
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)
void setBlockWeight(const std::string &blockId, double weight)
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 getElementVerticesNoResize(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
void getAllEdges(const std::string &edgeBlockName, std::vector< stk::mesh::Entity > &edges) const
void setInitialStateTime(double value)
std::size_t getNumNodesets() const
get the side set count
Teuchos::RCP< stk::mesh::MetaData > metaData_
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedElementVector_
void getEdgeBlockNames(std::vector< std::string > &names) const
std::map< std::string, Teuchos::RCP< shards::CellTopology > > edgeBlockCT_
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
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)
void setUseLowerCaseForIO(bool useLowerCase)
bool validBlockId(const std::string &blockId) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
void getSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, ArrayT &solutionValues) const
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_
void setCellFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
std::unordered_map< stk::mesh::EntityId, std::size_t > localEdgeIDHash_
ProcIdFieldType * getProcessorIdField()
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
stk::mesh::EntityId elementGlobalId(stk::mesh::Entity elmt) 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_
unsigned entityOwnerRank(stk::mesh::Entity entity) const
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
void addPeriodicBCs(const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &bc_vec)
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)
stk::mesh::EntityId faceGlobalId(stk::mesh::Entity face) const
bool operator()(stk::mesh::Entity a, stk::mesh::Entity b)
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_
const VectorFieldType & getEdgesField() const
std::map< std::string, stk::mesh::Part * > nodesets_
void getNeighborElements(std::vector< stk::mesh::Entity > &elements) const