Panzer Version of the Day
Loading...
Searching...
No Matches
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
72namespace percept {
73 class PerceptMesh;
74 class URP_Heterogeneous_3D;
75}
76#endif
77
78namespace panzer_stk {
79
80class PeriodicBC_MatcherBase;
81
86public:
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_; }
92protected:
93 stk::mesh::EntityId gid_;
94 std::vector<stk::mesh::EntityId> nodes_;
95
97};
98
101Teuchos::RCP<ElementDescriptor>
102buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes);
103
105public:
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
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 void addEdgeBlock(const std::string & elemBlockName,
143 const std::string & edgeBlockName,
144 const stk::topology & topology);
145 void addEdgeBlock(const std::string & elemBlockName,
146 const std::string & edgeBlockName,
147 const CellTopologyData * ctData);
148
151 void addFaceBlock(const std::string & elemBlockName,
152 const std::string & faceBlockName,
153 const stk::topology & topology);
154 void addFaceBlock(const std::string & elemBlockName,
155 const std::string & faceBlockName,
156 const CellTopologyData * ctData);
157
160 void addSideset(const std::string & name,const CellTopologyData * ctData);
161
164 void addNodeset(const std::string & name);
165
168 void addSolutionField(const std::string & fieldName,const std::string & blockId);
169
172 void addCellField(const std::string & fieldName,const std::string & blockId);
173
176 void addEdgeField(const std::string & fieldName,const std::string & blockId);
177
180 void addFaceField(const std::string & fieldName,const std::string & blockId);
181
190 void addMeshCoordFields(const std::string & blockId,
191 const std::vector<std::string> & coordField,
192 const std::string & dispPrefix);
193
198 void addInformationRecords(const std::vector<std::string> & info_records);
199
201
213 void initialize(stk::ParallelMachine parallelMach,bool setupIO=true,
214 const bool buildRefinementSupport = false);
215
221 void instantiateBulkData(stk::ParallelMachine parallelMach);
222
223 // functions to manage and manipulate bulk data
225
228 void beginModification();
229
232 void endModification();
233
239 void addNode(stk::mesh::EntityId gid, const std::vector<double> & coord);
240
241 void addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block);
242
243 void addEdges();
244
245 void addFaces();
246
249 void addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset);
250
253 void addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset);
254
257 void addEntityToEdgeBlock(stk::mesh::Entity entity,stk::mesh::Part * edgeblock);
260 void addEntitiesToEdgeBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * edgeblock);
261
264 void addEntityToFaceBlock(stk::mesh::Entity entity,stk::mesh::Part * faceblock);
267 void addEntitiesToFaceBlock(std::vector<stk::mesh::Entity> entities,stk::mesh::Part * faceblock);
268
269 // Methods to interrogate the mesh topology and structure
271
275 { return *coordinatesField_; }
276
280 { return *edgesField_; }
281
283 { return *facesField_; }
284
287 const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const;
288
291 const double * getNodeCoordinates(stk::mesh::Entity node) const;
292
295 void getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
296 std::vector<stk::mesh::EntityId> & subcellIds) const;
297
300 void getMyElements(std::vector<stk::mesh::Entity> & elements) const;
301
304 void getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
305
309 void getNeighborElements(std::vector<stk::mesh::Entity> & elements) const;
310
313 void getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
314
317 void getMyEdges(std::vector<stk::mesh::Entity> & edges) const;
318
326 void getMyEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const;
327
336 void getMyEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const;
337
345 void getAllEdges(const std::string & edgeBlockName,std::vector<stk::mesh::Entity> & edges) const;
346
355 void getAllEdges(const std::string & edgeBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & edges) const;
356
359 void getMyFaces(std::vector<stk::mesh::Entity> & faces) const;
360
368 void getMyFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const;
369
378 void getMyFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const;
379
387 void getAllFaces(const std::string & faceBlockName,std::vector<stk::mesh::Entity> & faces) const;
388
397 void getAllFaces(const std::string & faceBlockName,const std::string & blockName,std::vector<stk::mesh::Entity> & faces) const;
398
406 void getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
407
416 void getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
417
425 void getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
426
435 void getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
436
446 void getMyNodes(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const;
447
456 stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const;
457
458 // Utility functions
460
469 void
470 writeToExodus(const std::string& filename,
471 const bool append = false);
472
490 void
491 setupExodusFile(const std::string& filename,
492 const bool append = false,
493 const bool append_after_restart_time = false,
494 const double restart_time = 0.0);
495
508 void
510 double timestep);
511
531 void
533 const std::string& key,
534 const int& value);
535
555 void
557 const std::string& key,
558 const double& value);
559
579 void
581 const std::string& key,
582 const std::vector<int>& value);
583
603 void
605 const std::string& key,
606 const std::vector<double>& value);
607
608 // Accessor functions
610
612 Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;
613
614 Teuchos::RCP<stk::mesh::BulkData> getBulkData() const { return bulkData_; }
615 Teuchos::RCP<stk::mesh::MetaData> getMetaData() const { return metaData_; }
616
617#ifdef PANZER_HAVE_PERCEPT
619 Teuchos::RCP<percept::PerceptMesh> getRefinedMesh() const
620 { TEUCHOS_ASSERT(Teuchos::nonnull(refinedMesh_)); return refinedMesh_; }
621#endif
622
623 bool isWritable() const;
624
625 bool isModifiable() const
626 { if(bulkData_==Teuchos::null) return false;
627 return bulkData_->in_modifiable_state(); }
628
630 unsigned getDimension() const
631 { return dimension_; }
632
634 std::size_t getNumElementBlocks() const
635 { return elementBlocks_.size(); }
636
644 void getElementBlockNames(std::vector<std::string> & names) const;
645
653 void getSidesetNames(std::vector<std::string> & name) const;
654
662 void getNodesetNames(std::vector<std::string> & name) const;
663
671 void getEdgeBlockNames(std::vector<std::string> & names) const;
672
680 void getFaceBlockNames(std::vector<std::string> & names) const;
681
683 stk::mesh::Part * getOwnedPart() const
684 { return &getMetaData()->locally_owned_part(); } // I don't like the pointer access here, but it will do for now!
685
687 stk::mesh::Part * getElementBlockPart(const std::string & name) const
688 {
689 std::map<std::string, stk::mesh::Part*>::const_iterator itr = elementBlocks_.find(name); // Element blocks
690 if(itr==elementBlocks_.end()) return 0;
691 return itr->second;
692 }
693
695 stk::mesh::Part * getEdgeBlock(const std::string & name) const
696 {
697 std::map<std::string, stk::mesh::Part*>::const_iterator itr = edgeBlocks_.find(name); // edge blocks
698 if(itr==edgeBlocks_.end()) return 0;
699 return itr->second;
700 }
701
703 stk::mesh::Part * getFaceBlock(const std::string & name) const
704 {
705 std::map<std::string, stk::mesh::Part*>::const_iterator itr = faceBlocks_.find(name); // face blocks
706 if(itr==faceBlocks_.end()) return 0;
707 return itr->second;
708 }
709
711 std::size_t getNumSidesets() const
712 { return sidesets_.size(); }
713
714 stk::mesh::Part * getSideset(const std::string & name) const
715 {
716 auto itr = sidesets_.find(name);
717 return (itr != sidesets_.end()) ? itr->second : nullptr;
718 }
719
721 std::size_t getNumNodesets() const
722 { return nodesets_.size(); }
723
724 stk::mesh::Part * getNodeset(const std::string & name) const
725 {
726 auto itr = nodesets_.find(name);
727 return (itr != nodesets_.end()) ? itr->second : nullptr;
728 }
729
731 std::size_t getEntityCounts(unsigned entityRank) const;
732
734 stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const;
735
736 // Utilities
738
740 void getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const;
741
743 void getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const;
744
747 void getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
748 std::vector<int> & relIds) const;
749
752 void getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
753 std::vector<int> & relIds, unsigned int matchType) const;
754
755
757 void getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeId,std::vector<stk::mesh::Entity> & elements) const;
758
760 void buildSubcells();
761
764 std::size_t elementLocalId(stk::mesh::Entity elmt) const;
765
768 std::size_t elementLocalId(stk::mesh::EntityId gid) const;
769
772 inline stk::mesh::EntityId elementGlobalId(std::size_t lid) const
773 { return bulkData_->identifier((*orderedElementVector_)[lid]); }
774
777 inline stk::mesh::EntityId elementGlobalId(stk::mesh::Entity elmt) const
778 { return bulkData_->identifier(elmt); }
779
782 bool isEdgeLocal(stk::mesh::Entity edge) const;
783
786 bool isEdgeLocal(stk::mesh::EntityId gid) const;
787
790 std::size_t edgeLocalId(stk::mesh::Entity elmt) const;
791
794 std::size_t edgeLocalId(stk::mesh::EntityId gid) const;
795
798 inline stk::mesh::EntityId edgeGlobalId(std::size_t lid) const
799 { return bulkData_->identifier((*orderedEdgeVector_)[lid]); }
800
803 inline stk::mesh::EntityId edgeGlobalId(stk::mesh::Entity edge) const
804 { return bulkData_->identifier(edge); }
805
808 bool isFaceLocal(stk::mesh::Entity face) const;
809
812 bool isFaceLocal(stk::mesh::EntityId gid) const;
813
816 std::size_t faceLocalId(stk::mesh::Entity elmt) const;
817
820 std::size_t faceLocalId(stk::mesh::EntityId gid) const;
821
824 inline stk::mesh::EntityId faceGlobalId(std::size_t lid) const
825 { return bulkData_->identifier((*orderedFaceVector_)[lid]); }
826
829 inline stk::mesh::EntityId faceGlobalId(stk::mesh::Entity face) const
830 { return bulkData_->identifier(face); }
831
834 inline unsigned entityOwnerRank(stk::mesh::Entity entity) const
835 { return bulkData_->parallel_owner_rank(entity); }
836
839 inline bool isValid(stk::mesh::Entity entity) const
840 { return bulkData_->is_valid(entity); }
841
844 std::string containingBlockId(stk::mesh::Entity elmt) const;
845
850 stk::mesh::Field<double> * getSolutionField(const std::string & fieldName,
851 const std::string & blockId) const;
852
857 stk::mesh::Field<double> * getCellField(const std::string & fieldName,
858 const std::string & blockId) const;
859
864 stk::mesh::Field<double> * getEdgeField(const std::string & fieldName,
865 const std::string & blockId) const;
866
871 stk::mesh::Field<double> * getFaceField(const std::string & fieldName,
872 const std::string & blockId) const;
873
875
877 bool isInitialized() const { return initialized_; }
878
882 Teuchos::RCP<const std::vector<stk::mesh::Entity> > getElementsOrderedByLID() const;
883
898 template <typename ArrayT>
899 void setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
900 const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
901
916 template <typename ArrayT>
917 void getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
918 const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const;
919
934 template <typename ArrayT>
935 void setCellFieldData(const std::string & fieldName,const std::string & blockId,
936 const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
937
941 Teuchos::RCP<const std::vector<stk::mesh::Entity> > getEdgesOrderedByLID() const;
942
946 Teuchos::RCP<const std::vector<stk::mesh::Entity> > getFacesOrderedByLID() const;
947
962 template <typename ArrayT>
963 void setEdgeFieldData(const std::string & fieldName,const std::string & blockId,
964 const std::vector<std::size_t> & localEdgeIds,const ArrayT & edgeValues,double scaleValue=1.0);
965
980 template <typename ArrayT>
981 void setFaceFieldData(const std::string & fieldName,const std::string & blockId,
982 const std::vector<std::size_t> & localFaceIds,const ArrayT & faceValues,double scaleValue=1.0);
983
992 template <typename ArrayT>
993 void getElementVertices(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
994
1003 template <typename ArrayT>
1004 void getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1005
1015 template <typename ArrayT>
1016 void getElementVertices(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
1017
1027 template <typename ArrayT>
1028 void getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1029
1038 template <typename ArrayT>
1039 void getElementVerticesNoResize(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
1040
1049 template <typename ArrayT>
1050 void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1051
1061 template <typename ArrayT>
1062 void getElementVerticesNoResize(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
1063
1073 template <typename ArrayT>
1074 void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1075
1076 // const stk::mesh::FEMInterface & getFEMInterface() const
1077 // { return *femPtr_; }
1078
1079 stk::mesh::EntityRank getElementRank() const { return stk::topology::ELEMENT_RANK; }
1080 stk::mesh::EntityRank getSideRank() const { return metaData_->side_rank(); }
1081 stk::mesh::EntityRank getFaceRank() const { return stk::topology::FACE_RANK; }
1082 stk::mesh::EntityRank getEdgeRank() const { return stk::topology::EDGE_RANK; }
1083 stk::mesh::EntityRank getNodeRank() const { return stk::topology::NODE_RANK; }
1084
1088
1091 void buildLocalElementIDs();
1092
1095 void buildLocalEdgeIDs();
1096
1099 void buildLocalFaceIDs();
1100
1103 const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
1105 { return periodicBCs_; }
1106
1109 std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
1111 { return periodicBCs_; }
1112
1115 const bool & useBoundingBoxSearch() const
1116 { return useBBoxSearch_; }
1117
1119 void setBoundingBoxSearchFlag(const bool & searchFlag)
1120 { useBBoxSearch_ = searchFlag; return; }
1121
1127 void addPeriodicBC(const Teuchos::RCP<const PeriodicBC_MatcherBase> & bc)
1128 { periodicBCs_.push_back(bc); }
1129
1135 void addPeriodicBCs(const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bc_vec)
1136 { periodicBCs_.insert(periodicBCs_.end(),bc_vec.begin(),bc_vec.end()); }
1137
1140 std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
1141 getPeriodicNodePairing() const;
1142
1145 bool validBlockId(const std::string & blockId) const;
1146
1149 void print(std::ostream & os) const;
1150
1153 void printMetaData(std::ostream & os) const;
1154
1157 Teuchos::RCP<const shards::CellTopology> getCellTopology(const std::string & eBlock) const;
1158
1163 double getCurrentStateTime() const { return currentStateTime_; }
1164
1170 double getInitialStateTime() const { return initialStateTime_; }
1171
1176 void setInitialStateTime(double value) { initialStateTime_ = value; }
1177
1180 void rebalance(const Teuchos::ParameterList & params);
1181
1185 void setBlockWeight(const std::string & blockId,double weight)
1186 { blockWeights_[blockId] = weight; }
1187
1194 void setUseFieldCoordinates(bool useFieldCoordinates)
1195 { useFieldCoordinates_ = useFieldCoordinates; }
1196
1199 { return useFieldCoordinates_; }
1200
1202 void setUseLowerCaseForIO(bool useLowerCase)
1203 { useLowerCase_ = useLowerCase; }
1204
1207 { return useLowerCase_; }
1208
1219 template <typename ArrayT>
1220 void getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
1221
1222 template <typename ArrayT>
1223 void getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
1224 const std::string & eBlock, ArrayT & vertices) const;
1225
1235 template <typename ArrayT>
1236 void getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1237
1238 template <typename ArrayT>
1239 void getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
1240
1246 void refineMesh(const int numberOfLevels, const bool deleteParentElements);
1247
1248public: // static operations
1249 static const std::string coordsString;
1250 static const std::string nodesString;
1251 static const std::string edgesString;
1252 static const std::string edgeBlockString;
1253 static const std::string faceBlockString;
1254 static const std::string facesString;
1255
1256protected:
1257
1260 void buildEntityCounts();
1261
1264 void buildMaxEntityIds();
1265
1269 void initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
1270 bool setupIO);
1271
1276 Teuchos::RCP<Teuchos::MpiComm<int> > getSafeCommunicator(stk::ParallelMachine parallelMach) const;
1277
1284
1288 bool isMeshCoordField(const std::string & eBlock,const std::string & fieldName,int & axis) const;
1289
1305 template <typename ArrayT>
1306 void setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
1307 const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues);
1308
1309 std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > periodicBCs_;
1310 bool useBBoxSearch_ = false; // TODO swap this to change default periodic BC search (see also PeriodicBC_Parser.cpp)
1311
1312 Teuchos::RCP<stk::mesh::MetaData> metaData_;
1313 Teuchos::RCP<stk::mesh::BulkData> bulkData_;
1314#ifdef PANZER_HAVE_PERCEPT
1315 Teuchos::RCP<percept::PerceptMesh> refinedMesh_;
1316 Teuchos::RCP<percept::URP_Heterogeneous_3D> breakPattern_;
1317#endif
1318
1319 std::map<std::string, stk::mesh::Part*> elementBlocks_; // Element blocks
1320 std::map<std::string, stk::mesh::Part*> sidesets_; // Side sets
1321 std::map<std::string, stk::mesh::Part*> nodesets_; // Node sets
1322 std::map<std::string, stk::mesh::Part*> edgeBlocks_; // Edge blocks
1323 std::map<std::string, stk::mesh::Part*> faceBlocks_; // Face blocks
1324
1325 std::map<std::string, Teuchos::RCP<shards::CellTopology> > elementBlockCT_;
1326
1327 // for storing/accessing nodes
1328 stk::mesh::Part * nodesPart_;
1329 std::vector<stk::mesh::Part*> nodesPartVec_;
1330 stk::mesh::Part * edgesPart_;
1331 std::vector<stk::mesh::Part*> edgesPartVec_;
1332 stk::mesh::Part * facesPart_;
1333 std::vector<stk::mesh::Part*> facesPartVec_;
1334
1340
1341 // maps field names to solution field stk mesh handles
1342 std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToSolution_;
1343 std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToCellField_;
1344 std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToEdgeField_;
1345 std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToFaceField_;
1346
1347 // use a set to maintain a list of unique information records
1348 std::set<std::string> informationRecords_;
1349
1350 unsigned dimension_;
1351
1353
1354 // how many elements, faces, edges, and nodes are there globally
1355 std::vector<std::size_t> entityCounts_;
1356
1357 // what is maximum entity ID
1358 std::vector<stk::mesh::EntityId> maxEntityId_;
1359
1360 unsigned procRank_;
1361 std::size_t currentLocalId_;
1362
1363 Teuchos::RCP<Teuchos::MpiComm<int> > mpiComm_;
1364
1365 double initialStateTime_; // the time stamp at the time this object was constructed (default 0.0)
1366 double currentStateTime_; // the time stamp set by the user when writeToExodus is called (default 0.0)
1367
1368#ifdef PANZER_HAVE_IOSS
1369 // I/O support
1370 Teuchos::RCP<stk::io::StkMeshIoBroker> meshData_;
1371 int meshIndex_;
1372
1377 enum class GlobalVariable
1378 {
1379 ADD,
1380 WRITE
1381 }; // end of enum class GlobalVariable
1382
1399 void
1400 globalToExodus(
1401 const GlobalVariable& flag);
1402
1406 Teuchos::ParameterList globalData_;
1407#endif
1408
1409 // uses lazy evaluation
1410 mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedElementVector_;
1411
1412 // uses lazy evaluation
1413 mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedEdgeVector_;
1414
1415 // uses lazy evaluation
1416 mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedFaceVector_;
1417
1418 // for element block weights
1419 std::map<std::string,double> blockWeights_;
1420
1421 std::unordered_map<stk::mesh::EntityId,std::size_t> localIDHash_;
1422 std::unordered_map<stk::mesh::EntityId,std::size_t> localEdgeIDHash_;
1423 std::unordered_map<stk::mesh::EntityId,std::size_t> localFaceIDHash_;
1424
1425 // Store mesh displacement fields by element block. This map
1426 // goes like this meshCoordFields_[eBlock][axis_index] => coordinate FieldName
1427 // goes like this meshDispFields_[eBlock][axis_index] => displacement FieldName
1428 std::map<std::string,std::vector<std::string> > meshCoordFields_; // coordinate fields written by user
1429 std::map<std::string,std::vector<std::string> > meshDispFields_; // displacement fields, output to exodus
1430
1432
1434
1435 // Object describing how to sort a vector of elements using
1436 // local ID as the key, very short lived object
1438 public:
1439 LocalIdCompare(const STK_Interface * mesh) : mesh_(mesh) {}
1440
1441 // Compares two stk mesh entities based on local ID
1442 bool operator() (stk::mesh::Entity a,stk::mesh::Entity b)
1443 { return mesh_->elementLocalId(a) < mesh_->elementLocalId(b);}
1444
1445 private:
1447 };
1448};
1449
1450template <typename ArrayT>
1451void STK_Interface::setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
1452 const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
1453{
1454 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1455 auto solutionValues_h = Kokkos::create_mirror_view(solutionValues);
1456 Kokkos::deep_copy(solutionValues_h, solutionValues);
1457
1458 int field_axis = -1;
1459 if(isMeshCoordField(blockId,fieldName,field_axis)) {
1460 setDispFieldData(fieldName,blockId,field_axis,localElementIds,solutionValues_h);
1461 return;
1462 }
1463
1464 SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1465
1466 for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1467 std::size_t localId = localElementIds[cell];
1468 stk::mesh::Entity element = elements[localId];
1469
1470 // loop over nodes set solution values
1471 const size_t num_nodes = bulkData_->num_nodes(element);
1472 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1473 for(std::size_t i=0; i<num_nodes; ++i) {
1474 stk::mesh::Entity node = nodes[i];
1475
1476 double * solnData = stk::mesh::field_data(*field,node);
1477 // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1478 solnData[0] = scaleValue*solutionValues_h(cell,i);
1479 }
1480 }
1481}
1482
1483template <typename ArrayT>
1484void STK_Interface::setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
1485 const std::vector<std::size_t> & localElementIds,const ArrayT & dispValues)
1486{
1487 TEUCHOS_ASSERT(axis>=0); // sanity check
1488
1489 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1490
1491 SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1492 const VectorFieldType & coord_field = this->getCoordinatesField();
1493
1494 for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1495 std::size_t localId = localElementIds[cell];
1496 stk::mesh::Entity element = elements[localId];
1497
1498 // loop over nodes set solution values
1499 const size_t num_nodes = bulkData_->num_nodes(element);
1500 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1501 for(std::size_t i=0; i<num_nodes; ++i) {
1502 stk::mesh::Entity node = nodes[i];
1503
1504 double * solnData = stk::mesh::field_data(*field,node);
1505 double * coordData = stk::mesh::field_data(coord_field,node);
1506 // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1507 solnData[0] = dispValues(cell,i)-coordData[axis];
1508 }
1509 }
1510}
1511
1512template <typename ArrayT>
1513void STK_Interface::getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
1514 const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const
1515{
1516 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1517
1518 solutionValues = Kokkos::createDynRankView(solutionValues,
1519 "solutionValues",
1520 localElementIds.size(),
1521 bulkData_->num_nodes(elements[localElementIds[0]]));
1522
1523 SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
1524
1525 for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1526 std::size_t localId = localElementIds[cell];
1527 stk::mesh::Entity element = elements[localId];
1528
1529 // loop over nodes set solution values
1530 const size_t num_nodes = bulkData_->num_nodes(element);
1531 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1532 for(std::size_t i=0; i<num_nodes; ++i) {
1533 stk::mesh::Entity node = nodes[i];
1534
1535 double * solnData = stk::mesh::field_data(*field,node);
1536 // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1537 solutionValues(cell,i) = solnData[0];
1538 }
1539 }
1540}
1541
1542template <typename ArrayT>
1543void STK_Interface::setCellFieldData(const std::string & fieldName,const std::string & blockId,
1544 const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
1545{
1546 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1547
1548 SolutionFieldType * field = this->getCellField(fieldName,blockId);
1549
1550 auto solutionValues_h = Kokkos::create_mirror_view(solutionValues);
1551 Kokkos::deep_copy(solutionValues_h, solutionValues);
1552
1553 for(std::size_t cell=0;cell<localElementIds.size();cell++) {
1554 std::size_t localId = localElementIds[cell];
1555 stk::mesh::Entity element = elements[localId];
1556
1557 double * solnData = stk::mesh::field_data(*field,element);
1558 TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1559 solnData[0] = scaleValue*solutionValues_h.access(cell,0);
1560 }
1561}
1562
1563template <typename ArrayT>
1564void STK_Interface::setEdgeFieldData(const std::string & fieldName,const std::string & blockId,
1565 const std::vector<std::size_t> & localEdgeIds,const ArrayT & edgeValues,double scaleValue)
1566{
1567 const std::vector<stk::mesh::Entity> & edges = *(this->getEdgesOrderedByLID());
1568
1569 SolutionFieldType * field = this->getEdgeField(fieldName,blockId);
1570
1571 auto edgeValues_h = Kokkos::create_mirror_view(edgeValues);
1572 Kokkos::deep_copy(edgeValues_h, edgeValues);
1573
1574 for(std::size_t idx=0;idx<localEdgeIds.size();idx++) {
1575 std::size_t localId = localEdgeIds[idx];
1576 stk::mesh::Entity edge = edges[localId];
1577
1578 double * solnData = stk::mesh::field_data(*field,edge);
1579 TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1580 solnData[0] = scaleValue*edgeValues_h.access(idx,0);
1581 }
1582}
1583
1584template <typename ArrayT>
1585void STK_Interface::setFaceFieldData(const std::string & fieldName,const std::string & blockId,
1586 const std::vector<std::size_t> & localFaceIds,const ArrayT & faceValues,double scaleValue)
1587{
1588 const std::vector<stk::mesh::Entity> & faces = *(this->getFacesOrderedByLID());
1589
1590 SolutionFieldType * field = this->getFaceField(fieldName,blockId);
1591
1592 auto faceValues_h = Kokkos::create_mirror_view(faceValues);
1593 Kokkos::deep_copy(faceValues_h, faceValues);
1594
1595 for(std::size_t idx=0;idx<localFaceIds.size();idx++) {
1596 std::size_t localId = localFaceIds[idx];
1597 stk::mesh::Entity face = faces[localId];
1598
1599 double * solnData = stk::mesh::field_data(*field,face);
1600 TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
1601 solnData[0] = scaleValue*faceValues_h.access(idx,0);
1602 }
1603}
1604
1605template <typename ArrayT>
1606void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
1607{
1609 //
1610 // gather from the intrinsic mesh coordinates (non-lagrangian)
1611 //
1612
1613 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1614
1615 // convert to a vector of entity objects
1616 std::vector<stk::mesh::Entity> selected_elements;
1617 for(std::size_t cell=0;cell<localElementIds.size();cell++)
1618 selected_elements.push_back(elements[localElementIds[cell]]);
1619
1620 getElementVertices_FromCoords(selected_elements,vertices);
1621 }
1622 else {
1623 TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1624 "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
1625 "without specifying an element block.");
1626 }
1627}
1628
1629template <typename ArrayT>
1630void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1631{
1633 getElementVertices_FromCoords(elements,vertices);
1634 }
1635 else {
1636 TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1637 "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
1638 "without specifying an element block.");
1639 }
1640}
1641
1642template <typename ArrayT>
1643void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1644{
1646 getElementVertices_FromCoords(elements,vertices);
1647 }
1648 else {
1649 getElementVertices_FromField(elements,eBlock,vertices);
1650 }
1651}
1652
1653template <typename ArrayT>
1654void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
1655{
1656 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1657
1658 // convert to a vector of entity objects
1659 std::vector<stk::mesh::Entity> selected_elements;
1660 for(std::size_t cell=0;cell<localElementIds.size();cell++)
1661 selected_elements.push_back(elements[localElementIds[cell]]);
1662
1664 getElementVertices_FromCoords(selected_elements,vertices);
1665 }
1666 else {
1667 getElementVertices_FromField(selected_elements,eBlock,vertices);
1668 }
1669}
1670
1671template <typename ArrayT>
1672void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
1673{
1675 //
1676 // gather from the intrinsic mesh coordinates (non-lagrangian)
1677 //
1678
1679 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1680
1681 // convert to a vector of entity objects
1682 std::vector<stk::mesh::Entity> selected_elements;
1683 for(std::size_t cell=0;cell<localElementIds.size();cell++)
1684 selected_elements.push_back(elements[localElementIds[cell]]);
1685
1686 getElementVertices_FromCoordsNoResize(selected_elements,vertices);
1687 }
1688 else {
1689 TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1690 "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
1691 "without specifying an element block.");
1692 }
1693}
1694
1695template <typename ArrayT>
1696void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1697{
1699 getElementVertices_FromCoordsNoResize(elements,vertices);
1700 }
1701 else {
1702 TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
1703 "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
1704 "without specifying an element block.");
1705 }
1706}
1707
1708template <typename ArrayT>
1709void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1710{
1712 getElementVertices_FromCoordsNoResize(elements,vertices);
1713 }
1714 else {
1715 getElementVertices_FromFieldNoResize(elements,eBlock,vertices);
1716 }
1717}
1718
1719template <typename ArrayT>
1720void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
1721{
1722 const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
1723
1724 // convert to a vector of entity objects
1725 std::vector<stk::mesh::Entity> selected_elements;
1726 for(std::size_t cell=0;cell<localElementIds.size();cell++)
1727 selected_elements.push_back(elements[localElementIds[cell]]);
1728
1730 getElementVertices_FromCoordsNoResize(selected_elements,vertices);
1731 }
1732 else {
1733 getElementVertices_FromFieldNoResize(selected_elements,eBlock,vertices);
1734 }
1735}
1736
1737template <typename ArrayT>
1738void STK_Interface::getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1739{
1740 // nothing to do! silently return
1741 if(elements.size() == 0) {
1742 vertices = Kokkos::createDynRankView(vertices, "vertices", 0, 0, 0);
1743 return;
1744 }
1745
1746 //
1747 // gather from the intrinsic mesh coordinates (non-lagrangian)
1748 //
1749
1750 // get *master* cell toplogy...(belongs to first element)
1751 const auto masterVertexCount
1752 = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1753
1754 // allocate space
1755 vertices = Kokkos::createDynRankView(vertices, "vertices", elements.size(), masterVertexCount,getDimension());
1756 auto vertices_h = Kokkos::create_mirror_view(vertices);
1757 Kokkos::deep_copy(vertices_h, vertices);
1758
1759 // loop over each requested element
1760 const auto dim = getDimension();
1761 for(std::size_t cell = 0; cell < elements.size(); cell++) {
1762 const auto element = elements[cell];
1763 TEUCHOS_ASSERT(element != 0);
1764
1765 const auto vertexCount
1766 = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
1767 TEUCHOS_TEST_FOR_EXCEPTION(vertexCount != masterVertexCount, std::runtime_error,
1768 "In call to STK_Interface::getElementVertices all elements "
1769 "must have the same vertex count!");
1770
1771 // loop over all element nodes
1772 const size_t num_nodes = bulkData_->num_nodes(element);
1773 auto const* nodes = bulkData_->begin_nodes(element);
1774 TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
1775 "In call to STK_Interface::getElementVertices cardinality of "
1776 "element node relations must be the vertex count!");
1777 for(std::size_t node = 0; node < num_nodes; ++node) {
1778 const double * coord = getNodeCoordinates(nodes[node]);
1779
1780 // set each dimension of the coordinate
1781 for(unsigned d=0;d<dim;d++)
1782 vertices_h(cell,node,d) = coord[d];
1783 }
1784 }
1785 Kokkos::deep_copy(vertices, vertices_h);
1786}
1787
1788template <typename ArrayT>
1789void STK_Interface::getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
1790{
1791 // nothing to do! silently return
1792 if(elements.size()==0) {
1793 return;
1794 }
1795
1796 //
1797 // gather from the intrinsic mesh coordinates (non-lagrangian)
1798 //
1799
1800 // get *master* cell toplogy...(belongs to first element)
1801 unsigned masterVertexCount
1802 = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1803
1804 // loop over each requested element
1805 unsigned dim = getDimension();
1806 auto vertices_h = Kokkos::create_mirror_view(vertices);
1807 for(std::size_t cell=0;cell<elements.size();cell++) {
1808 stk::mesh::Entity element = elements[cell];
1809 TEUCHOS_ASSERT(element!=0);
1810
1811 unsigned vertexCount
1812 = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
1813 TEUCHOS_TEST_FOR_EXCEPTION(vertexCount!=masterVertexCount,std::runtime_error,
1814 "In call to STK_Interface::getElementVertices all elements "
1815 "must have the same vertex count!");
1816
1817 // loop over all element nodes
1818 const size_t num_nodes = bulkData_->num_nodes(element);
1819 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1820 TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
1821 "In call to STK_Interface::getElementVertices cardinality of "
1822 "element node relations must be the vertex count!");
1823 for(std::size_t node=0; node<num_nodes; ++node) {
1824 const double * coord = getNodeCoordinates(nodes[node]);
1825
1826 // set each dimension of the coordinate
1827 for(unsigned d=0;d<dim;d++)
1828 vertices_h(cell,node,d) = coord[d];
1829 }
1830 }
1831 Kokkos::deep_copy(vertices, vertices_h);
1832}
1833
1834template <typename ArrayT>
1835void STK_Interface::getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
1836{
1837 TEUCHOS_ASSERT(useFieldCoordinates_);
1838
1839 // nothing to do! silently return
1840 if(elements.size()==0) {
1841 vertices = Kokkos::createDynRankView(vertices,"vertices",0,0,0);
1842 return;
1843 }
1844
1845 // get *master* cell toplogy...(belongs to first element)
1846 unsigned masterVertexCount
1847 = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
1848
1849 // allocate space
1850 vertices = Kokkos::createDynRankView(vertices,"vertices",elements.size(),masterVertexCount,getDimension());
1851 auto vertices_h = Kokkos::create_mirror_view(vertices);
1852 std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
1853 if(itr==meshCoordFields_.end()) {
1854 // no coordinate field set for this element block
1855 TEUCHOS_ASSERT(false);
1856 }
1857
1858 const std::vector<std::string> & coordField = itr->second;
1859 std::vector<SolutionFieldType*> fields(getDimension());
1860 for(std::size_t d=0;d<fields.size();d++) {
1861 fields[d] = this->getSolutionField(coordField[d],eBlock);
1862 }
1863
1864 for(std::size_t cell=0;cell<elements.size();cell++) {
1865 stk::mesh::Entity element = elements[cell];
1866
1867 // loop over nodes set solution values
1868 const size_t num_nodes = bulkData_->num_nodes(element);
1869 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1870 for(std::size_t i=0; i<num_nodes; ++i) {
1871 stk::mesh::Entity node = nodes[i];
1872
1873 const double * coord = getNodeCoordinates(node);
1874
1875 for(unsigned d=0;d<getDimension();d++) {
1876 double * solnData = stk::mesh::field_data(*fields[d],node);
1877
1878 // recall mesh field coordinates are stored as displacements
1879 // from the mesh coordinates, make sure to add them together
1880 vertices_h(cell,i,d) = solnData[0]+coord[d];
1881 }
1882 }
1883 }
1884 Kokkos::deep_copy(vertices, vertices_h);
1885}
1886
1887template <typename ArrayT>
1888void STK_Interface::getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
1889 const std::string & eBlock, ArrayT & vertices) const
1890{
1891 TEUCHOS_ASSERT(useFieldCoordinates_);
1892
1893 // nothing to do! silently return
1894 if(elements.size()==0) {
1895 return;
1896 }
1897
1898 std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
1899 if(itr==meshCoordFields_.end()) {
1900 // no coordinate field set for this element block
1901 TEUCHOS_ASSERT(false);
1902 }
1903
1904 const std::vector<std::string> & coordField = itr->second;
1905 std::vector<SolutionFieldType*> fields(getDimension());
1906 for(std::size_t d=0;d<fields.size();d++) {
1907 fields[d] = this->getSolutionField(coordField[d],eBlock);
1908 }
1909
1910 for(std::size_t cell=0;cell<elements.size();cell++) {
1911 stk::mesh::Entity element = elements[cell];
1912
1913 // loop over nodes set solution values
1914 const size_t num_nodes = bulkData_->num_nodes(element);
1915 stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
1916 for(std::size_t i=0; i<num_nodes; ++i) {
1917 stk::mesh::Entity node = nodes[i];
1918
1919 const double * coord = getNodeCoordinates(node);
1920
1921 for(unsigned d=0;d<getDimension();d++) {
1922 double * solnData = stk::mesh::field_data(*fields[d],node);
1923
1924 // recall mesh field coordinates are stored as displacements
1925 // from the mesh coordinates, make sure to add them together
1926 vertices(cell,i,d) = solnData[0]+coord[d];
1927 }
1928 }
1929 }
1930}
1931
1932}
1933
1934#endif
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
const std::vector< stk::mesh::EntityId > & getNodes() const
std::vector< stk::mesh::EntityId > nodes_
stk::mesh::EntityId getGID() const
bool operator()(stk::mesh::Entity a, stk::mesh::Entity b)
void rebalance(const Teuchos::ParameterList &params)
const VectorFieldType & getEdgesField() const
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
Teuchos::RCP< stk::mesh::MetaData > metaData_
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedFaceVector_
std::pair< Teuchos::RCP< std::vector< std::pair< std::size_t, std::size_t > > >, Teuchos::RCP< std::vector< unsigned int > > > getPeriodicNodePairing() const
void setFaceFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localFaceIds, const ArrayT &faceValues, double scaleValue=1.0)
bool isValid(stk::mesh::Entity entity) const
const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector() const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
stk::mesh::Field< double > SolutionFieldType
void getMyFaces(std::vector< stk::mesh::Entity > &faces) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
static const std::string edgesString
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.
static const std::string edgeBlockString
void getElementVertices(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
void getElementVertices_FromCoords(const std::vector< stk::mesh::Entity > &elements, ArrayT &vertices) const
std::size_t getNumSidesets() const
get the side set count
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedElementVector_
std::map< std::string, stk::mesh::Part * > sidesets_
void addPeriodicBC(const Teuchos::RCP< const PeriodicBC_MatcherBase > &bc)
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::unordered_map< stk::mesh::EntityId, std::size_t > localFaceIDHash_
void getElementsSharingNodes(const std::vector< stk::mesh::EntityId > nodeId, std::vector< stk::mesh::Entity > &elements) const
get a set of elements sharing multiple nodes
std::size_t getEntityCounts(unsigned entityRank) const
get the global counts for the entity of specified rank
void getFaceBlockNames(std::vector< std::string > &names) const
stk::mesh::Field< double, stk::mesh::Cartesian > VectorFieldType
bool isFaceLocal(stk::mesh::Entity face) const
void addEntityToEdgeBlock(stk::mesh::Entity entity, stk::mesh::Part *edgeblock)
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getElementsOrderedByLID() const
bool isInitialized() const
Has initialize been called on this mesh object?
std::map< std::string, stk::mesh::Part * > faceBlocks_
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
stk::mesh::EntityRank getNodeRank() const
void addInformationRecords(const std::vector< std::string > &info_records)
std::string containingBlockId(stk::mesh::Entity elmt) const
void getEdgeBlockNames(std::vector< std::string > &names) const
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
void getElementVertices_FromField(const std::vector< stk::mesh::Entity > &elements, const std::string &eBlock, ArrayT &vertices) const
const VectorFieldType & getFacesField() const
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
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getEdgesOrderedByLID() const
void buildSubcells()
force the mesh to build subcells: edges and faces
void setCellFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const
get max entity ID of type entityRank
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > & getPeriodicBCVector()
stk::mesh::EntityRank getElementRank() const
void getOwnedElementsSharingNode(stk::mesh::Entity node, std::vector< stk::mesh::Entity > &elements, std::vector< int > &relIds) const
const double * getNodeCoordinates(stk::mesh::EntityId nodeId) 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 setBlockWeight(const std::string &blockId, double weight)
void addEdgeField(const std::string &fieldName, const std::string &blockId)
stk::mesh::EntityId faceGlobalId(stk::mesh::Entity face) const
std::unordered_map< stk::mesh::EntityId, std::size_t > localEdgeIDHash_
void addEntityToFaceBlock(stk::mesh::Entity entity, stk::mesh::Part *faceblock)
void getNeighborElements(std::vector< stk::mesh::Entity > &elements) const
void getMyNodes(const std::string &sideName, const std::string &blockName, std::vector< stk::mesh::Entity > &nodes) const
std::vector< stk::mesh::Part * > edgesPartVec_
const bool & useBoundingBoxSearch() const
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 addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void getElementVerticesNoResize(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
bool isEdgeLocal(stk::mesh::Entity edge) const
void addPeriodicBCs(const std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &bc_vec)
void refineMesh(const int numberOfLevels, const bool deleteParentElements)
std::map< std::string, double > blockWeights_
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToEdgeField_
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToFaceField_
void setSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void addCellField(const std::string &fieldName, const std::string &blockId)
void setInitialStateTime(double value)
stk::mesh::Field< double > * getSolutionField(const std::string &fieldName, const std::string &blockId) const
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToSolution_
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void getElementVertices_FromFieldNoResize(const std::vector< stk::mesh::Entity > &elements, const std::string &eBlock, ArrayT &vertices) const
void getMySides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void addMeshCoordFields(const std::string &blockId, const std::vector< std::string > &coordField, const std::string &dispPrefix)
stk::mesh::Field< double > * getEdgeField(const std::string &fieldName, const std::string &blockId) const
std::vector< stk::mesh::EntityId > maxEntityId_
void getElementBlockNames(std::vector< std::string > &names) const
Teuchos::RCP< stk::mesh::BulkData > bulkData_
void addNodeset(const std::string &name)
void getSidesetNames(std::vector< std::string > &name) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
stk::mesh::EntityRank getFaceRank() const
void getMyEdges(std::vector< stk::mesh::Entity > &edges) const
std::size_t getNumNodesets() const
get the side set count
std::set< std::string > informationRecords_
void getAllEdges(const std::string &edgeBlockName, std::vector< stk::mesh::Entity > &edges) const
std::vector< stk::mesh::Part * > nodesPartVec_
void initializeFieldsInSTK(const std::map< std::pair< std::string, std::string >, SolutionFieldType * > &nameToField, bool setupIO)
static const std::string nodesString
const VectorFieldType & getCoordinatesField() const
void setUseFieldCoordinates(bool useFieldCoordinates)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Field< double > * getFaceField(const std::string &fieldName, const std::string &blockId) const
static const std::string facesString
void writeToExodus(const std::string &filename, const bool append=false)
Write this mesh and associated fields to the given output file.
unsigned getDimension() const
get the dimension
bool isMeshCoordField(const std::string &eBlock, const std::string &fieldName, int &axis) const
std::map< std::string, stk::mesh::Part * > nodesets_
void getElementVertices_FromCoordsNoResize(const std::vector< stk::mesh::Entity > &elements, ArrayT &vertices) const
std::map< std::string, std::vector< std::string > > meshCoordFields_
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
std::size_t faceLocalId(stk::mesh::Entity elmt) const
void getAllFaces(const std::string &faceBlockName, std::vector< stk::mesh::Entity > &faces) const
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.
Teuchos::RCP< std::vector< stk::mesh::Entity > > orderedEdgeVector_
stk::mesh::Part * getNodeset(const std::string &name) const
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
stk::mesh::EntityId edgeGlobalId(std::size_t lid) const
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void getNodesetNames(std::vector< std::string > &name) const
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
stk::mesh::Field< ProcIdData > ProcIdFieldType
std::unordered_map< stk::mesh::EntityId, std::size_t > localIDHash_
void addSolutionField(const std::string &fieldName, const std::string &blockId)
void getSubcellIndices(unsigned entityRank, stk::mesh::EntityId elementId, std::vector< stk::mesh::EntityId > &subcellIds) const
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCs_
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
Teuchos::RCP< const std::vector< stk::mesh::Entity > > getFacesOrderedByLID() const
void setDispFieldData(const std::string &fieldName, const std::string &blockId, int axis, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues)
stk::mesh::EntityId elementGlobalId(stk::mesh::Entity elmt) const
std::map< std::string, Teuchos::RCP< shards::CellTopology > > elementBlockCT_
stk::mesh::EntityRank getSideRank() const
Teuchos::RCP< Teuchos::MpiComm< int > > mpiComm_
bool validBlockId(const std::string &blockId) const
stk::mesh::EntityId edgeGlobalId(stk::mesh::Entity edge) const
void addFaceField(const std::string &fieldName, const std::string &blockId)
std::map< std::string, stk::mesh::Part * > edgeBlocks_
void print(std::ostream &os) const
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
stk::mesh::Part * getOwnedPart() const
Get a pointer to the locally owned part.
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
ProcIdFieldType * getProcessorIdField()
void addFaceBlock(const std::string &elemBlockName, const std::string &faceBlockName, const stk::topology &topology)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
std::map< std::string, std::vector< std::string > > meshDispFields_
static const std::string faceBlockString
std::vector< stk::mesh::Part * > facesPartVec_
void printMetaData(std::ostream &os) const
void instantiateBulkData(stk::ParallelMachine parallelMach)
stk::mesh::EntityId faceGlobalId(std::size_t lid) const
std::map< std::pair< std::string, std::string >, SolutionFieldType * > fieldNameToCellField_
std::size_t edgeLocalId(stk::mesh::Entity elmt) const
std::map< std::string, stk::mesh::Part * > elementBlocks_
static const std::string coordsString
void getSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, ArrayT &solutionValues) const
void setUseLowerCaseForIO(bool useLowerCase)
Teuchos::RCP< Teuchos::MpiComm< int > > getSafeCommunicator(stk::ParallelMachine parallelMach) const
std::size_t getNumElementBlocks() const
get the block count
void setBoundingBoxSearchFlag(const bool &searchFlag)
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
Teuchos::RCP< ElementDescriptor > buildElementDescriptor(stk::mesh::EntityId elmtId, std::vector< stk::mesh::EntityId > &nodes)