Panzer_STK_Interface.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//           Panzer: A partial differential equation assembly
//       engine for strongly coupled complex multiphysics systems
//                 Copyright (2011) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
// Eric C. Cyr (eccyr@sandia.gov)
// ***********************************************************************
// @HEADER
 
#ifndef __Panzer_STK_Interface_hpp__
#define __Panzer_STK_Interface_hpp__
 
#include <Teuchos_RCP.hpp>
#include <Teuchos_DefaultMpiComm.hpp>
 
#include <stk_mesh/base/Types.hpp>
#include <stk_mesh/base/MetaData.hpp>
#include <stk_mesh/base/BulkData.hpp>
#include <stk_mesh/base/Field.hpp>
#include <stk_mesh/base/FieldBase.hpp>
#include <stk_mesh/base/MetaData.hpp>
#include <stk_mesh/base/CoordinateSystems.hpp>
 
#include "Kokkos_Core.hpp"
 
#include <Shards_CellTopology.hpp>
#include <Shards_CellTopologyData.h>
 
#include <PanzerAdaptersSTK_config.hpp>
#include <Kokkos_ViewFactory.hpp>
 
#include <unordered_map>
 
#ifdef PANZER_HAVE_IOSS
#include <stk_io/StkMeshIoBroker.hpp>
#endif
 
namespace panzer_stk {
 
class PeriodicBC_MatcherBase;
 
class ElementDescriptor {
public:
   ElementDescriptor(stk::mesh::EntityId gid,const std::vector<stk::mesh::EntityId> & nodes);
   virtual ~ElementDescriptor();
 
   stk::mesh::EntityId getGID() const { return gid_; }
   const std::vector<stk::mesh::EntityId> & getNodes() const { return nodes_; }
protected:
   stk::mesh::EntityId gid_;
   std::vector<stk::mesh::EntityId> nodes_;
 
   ElementDescriptor();
};
 
Teuchos::RCP<ElementDescriptor> 
buildElementDescriptor(stk::mesh::EntityId elmtId,std::vector<stk::mesh::EntityId> & nodes);
 
class STK_Interface {
public:
   typedef double ProcIdData; // ECC: Not sure why?
   typedef stk::mesh::Field<double> SolutionFieldType;
   typedef stk::mesh::Field<double,stk::mesh::Cartesian> VectorFieldType;
   typedef stk::mesh::Field<ProcIdData> ProcIdFieldType;
 
   // some simple exception classes
   struct ElementBlockException : public std::logic_error 
   { ElementBlockException(const std::string & what) : std::logic_error(what) {} };
 
   struct SidesetException : public std::logic_error 
   { SidesetException(const std::string & what) : std::logic_error(what) {} };
 
   STK_Interface();
 
   STK_Interface(unsigned dim);
 
   STK_Interface(Teuchos::RCP<stk::mesh::MetaData> metaData);
 
   // functions called before initialize
 
   void addElementBlock(const std::string & name,const CellTopologyData * ctData);
 
   void addSideset(const std::string & name,const CellTopologyData * ctData);
 
   void addNodeset(const std::string & name);
 
   void addSolutionField(const std::string & fieldName,const std::string & blockId);
 
   void addCellField(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 initialize(stk::ParallelMachine parallelMach,bool setupIO=true);
 
   void instantiateBulkData(stk::ParallelMachine parallelMach);
 
   // functions to manage and manipulate bulk data
  
   void beginModification();
 
   void endModification();
 
   void addNode(stk::mesh::EntityId gid, const std::vector<double> & coord);
 
   void addElement(const Teuchos::RCP<ElementDescriptor> & ed,stk::mesh::Part * block);
 
   void addEdges();
 
   void addFaces();
 
   void addEntityToSideset(stk::mesh::Entity entity,stk::mesh::Part * sideset);
 
   void addEntityToNodeset(stk::mesh::Entity entity,stk::mesh::Part * nodeset);
 
   // Methods to interrogate the mesh topology and structure
 
   const VectorFieldType & getCoordinatesField() const
   { return *coordinatesField_; }
 
   const VectorFieldType & getEdgesField() const
   { return *edgesField_; }
 
   const VectorFieldType & getFacesField() const
   { return *facesField_; }
 
   const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const;
 
   const double * getNodeCoordinates(stk::mesh::Entity node) const;
 
   void getSubcellIndices(unsigned entityRank,stk::mesh::EntityId elementId,
                          std::vector<stk::mesh::EntityId> & subcellIds) const;
 
   void getMyElements(std::vector<stk::mesh::Entity> & elements) const;
 
   void getMyElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
 
   void getNeighborElements(std::vector<stk::mesh::Entity> & elements) const;
 
   void getNeighborElements(const std::string & blockID,std::vector<stk::mesh::Entity> & elements) const;
 
   void getMySides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
 
   void getMySides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
 
   void getAllSides(const std::string & sideName,std::vector<stk::mesh::Entity> & sides) const;
 
   void getAllSides(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & sides) const;
 
   void getMyNodes(const std::string & sideName,const std::string & blockName,std::vector<stk::mesh::Entity> & nodes) const;
 
   stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const;
 
   // Utility functions
 
  void
  writeToExodus(
    const std::string& filename);
 
  void
  setupExodusFile(
    const std::string& filename);
 
  void
  writeToExodus(
    double timestep);
 
  void
  addGlobalToExodus(
    const std::string& key,
    const int&         value);
 
  void
  addGlobalToExodus(
    const std::string& key,
    const double&      value);
 
  void
  addGlobalToExodus(
    const std::string&      key,
    const std::vector<int>& value);
 
  void
  addGlobalToExodus(
    const std::string&         key,
    const std::vector<double>& value);
 
   // Accessor functions
 
   Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;
 
   Teuchos::RCP<stk::mesh::BulkData> getBulkData() const { return bulkData_; }
   Teuchos::RCP<stk::mesh::MetaData> getMetaData() const { return metaData_; }
 
   bool isWritable() const;
 
   bool isModifiable() const
   {  if(bulkData_==Teuchos::null) return false;
      return bulkData_->in_modifiable_state(); }
 
   unsigned getDimension() const
   { return dimension_; }
 
   std::size_t getNumElementBlocks() const
   { return elementBlocks_.size(); }
 
   void getElementBlockNames(std::vector<std::string> & names) const;
 
   void getSidesetNames(std::vector<std::string> & name) const;
 
   void getNodesetNames(std::vector<std::string> & name) const;
 
   stk::mesh::Part * getOwnedPart() const
   { return &getMetaData()->locally_owned_part(); } // I don't like the pointer access here, but it will do for now!
 
   stk::mesh::Part * getElementBlockPart(const std::string & name) const
   { 
      std::map<std::string, stk::mesh::Part*>::const_iterator itr = elementBlocks_.find(name);   // Element blocks
      if(itr==elementBlocks_.end()) return 0;
      return itr->second; 
   }
 
   std::size_t getNumSidesets() const
   { return sidesets_.size(); }
 
   stk::mesh::Part * getSideset(const std::string & name) const
   { return sidesets_.find(name)->second; }
 
   std::size_t getNumNodesets() const
   { return nodesets_.size(); }
 
   stk::mesh::Part * getNodeset(const std::string & name) const
   { return nodesets_.find(name)->second; }
 
   std::size_t getEntityCounts(unsigned entityRank) const;
 
   stk::mesh::EntityId getMaxEntityId(unsigned entityRank) const;
 
   // Utilities
 
   void getElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements) const;
 
   void getNodeIdsForElement(stk::mesh::Entity element,std::vector<stk::mesh::EntityId> & nodeIds) const;
 
   void getOwnedElementsSharingNode(stk::mesh::Entity node,std::vector<stk::mesh::Entity> & elements,
                                    std::vector<int> & relIds) const;
 
   void getOwnedElementsSharingNode(stk::mesh::EntityId nodeId,std::vector<stk::mesh::Entity> & elements,
                                    std::vector<int> & relIds, unsigned int matchType) const;
 
 
   void getElementsSharingNodes(const std::vector<stk::mesh::EntityId> nodeId,std::vector<stk::mesh::Entity> & elements) const;
 
   void buildSubcells();
 
   std::size_t elementLocalId(stk::mesh::Entity elmt) const;
 
   std::size_t elementLocalId(stk::mesh::EntityId gid) const;
 
   inline stk::mesh::EntityId elementGlobalId(std::size_t lid) const
   { return bulkData_->identifier((*orderedElementVector_)[lid]); }
 
   inline stk::mesh::EntityId elementGlobalId(stk::mesh::Entity elmt) const
   { return bulkData_->identifier(elmt); }
 
  inline unsigned entityOwnerRank(stk::mesh::Entity entity) const
  { return bulkData_->parallel_owner_rank(entity); }
 
  inline bool isValid(stk::mesh::Entity entity) const
  { return bulkData_->is_valid(entity); }
 
   std::string containingBlockId(stk::mesh::Entity elmt) const;
 
   stk::mesh::Field<double> * getSolutionField(const std::string & fieldName,
                                               const std::string & blockId) const;
 
   stk::mesh::Field<double> * getCellField(const std::string & fieldName,
                                           const std::string & blockId) const;
 
   ProcIdFieldType * getProcessorIdField() { return processorIdField_; }
 
   bool isInitialized() const { return initialized_; }
 
   Teuchos::RCP<const std::vector<stk::mesh::Entity> > getElementsOrderedByLID() const;
 
   template <typename ArrayT>
   void setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
                             const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
 
   template <typename ArrayT>
   void getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
                             const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const;
 
   template <typename ArrayT>
   void setCellFieldData(const std::string & fieldName,const std::string & blockId,
                         const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue=1.0);
 
   template <typename ArrayT>
   void getElementVertices(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVerticesNoResize(const std::vector<std::size_t> & localIds, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVerticesNoResize(const std::vector<std::size_t> & localIds,const std::string & eBlock, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
 
   // const stk::mesh::FEMInterface & getFEMInterface() const 
   // { return *femPtr_; }
 
   stk::mesh::EntityRank getElementRank() const { return stk::topology::ELEMENT_RANK; }
   stk::mesh::EntityRank getSideRank() const { return metaData_->side_rank(); }
   stk::mesh::EntityRank getFaceRank() const { return stk::topology::FACE_RANK; }
   stk::mesh::EntityRank getEdgeRank() const { return stk::topology::EDGE_RANK; }
   stk::mesh::EntityRank getNodeRank() const { return stk::topology::NODE_RANK; }
 
   void initializeFromMetaData();
 
   void buildLocalElementIDs();
 
   const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
   getPeriodicBCVector() const
   { return periodicBCs_; }
 
   std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > &
   getPeriodicBCVector()
   { return periodicBCs_; }
 
   void addPeriodicBC(const Teuchos::RCP<const PeriodicBC_MatcherBase> & bc)
   { periodicBCs_.push_back(bc); }
 
   void addPeriodicBCs(const std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > & bc_vec)
   { periodicBCs_.insert(periodicBCs_.end(),bc_vec.begin(),bc_vec.end()); }
 
   std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > >
   getPeriodicNodePairing() const;
 
   bool validBlockId(const std::string & blockId) const;
 
   void print(std::ostream & os) const;
 
   void printMetaData(std::ostream & os) const;
 
   Teuchos::RCP<const shards::CellTopology> getCellTopology(const std::string & eBlock) const;
 
   double getCurrentStateTime() const { return currentStateTime_; }
 
   double getInitialStateTime() const { return initialStateTime_; }
 
   void setInitialStateTime(double value) { initialStateTime_ = value; }
 
   void rebalance(const Teuchos::ParameterList & params);
 
   void setBlockWeight(const std::string & blockId,double weight)
   { blockWeights_[blockId] = weight; }
 
   void setUseFieldCoordinates(bool useFieldCoordinates)
   { useFieldCoordinates_ = useFieldCoordinates; }
 
   bool getUseFieldCoordinates() const
   { return useFieldCoordinates_; }
 
   void setUseLowerCaseForIO(bool useLowerCase) 
   { useLowerCase_ = useLowerCase; }
 
   bool getUseLowerCaseForIO() const
   { return useLowerCase_; }
 
   template <typename ArrayT>
   void getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
                                             const std::string & eBlock, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
 
   template <typename ArrayT>
   void getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const;
 
public: // static operations
   static const std::string coordsString;
   static const std::string nodesString;
   static const std::string edgesString;
   static const std::string facesString;
 
protected:
 
   void buildEntityCounts();
 
   void buildMaxEntityIds();
 
   void initializeFieldsInSTK(const std::map<std::pair<std::string,std::string>,SolutionFieldType*> & nameToField,
                              bool setupIO);
 
   Teuchos::RCP<Teuchos::MpiComm<int> > getSafeCommunicator(stk::ParallelMachine parallelMach) const;
 
   void applyElementLoadBalanceWeights();
 
   bool isMeshCoordField(const std::string & eBlock,const std::string & fieldName,int & axis) const;
 
   template <typename ArrayT>
   void setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
                         const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues);
 
   std::vector<Teuchos::RCP<const PeriodicBC_MatcherBase> > periodicBCs_;
 
   Teuchos::RCP<stk::mesh::MetaData> metaData_;
   Teuchos::RCP<stk::mesh::BulkData> bulkData_;
 
   std::map<std::string, stk::mesh::Part*> elementBlocks_;   // Element blocks
   std::map<std::string, stk::mesh::Part*> sidesets_; // Side sets 
   std::map<std::string, stk::mesh::Part*> nodesets_; // Node sets 
   std::map<std::string, Teuchos::RCP<shards::CellTopology> > elementBlockCT_;
 
   // for storing/accessing nodes
   stk::mesh::Part * nodesPart_;
   std::vector<stk::mesh::Part*> nodesPartVec_;
   stk::mesh::Part * edgesPart_;
   std::vector<stk::mesh::Part*> edgesPartVec_;
   stk::mesh::Part * facesPart_;
   std::vector<stk::mesh::Part*> facesPartVec_;
 
   VectorFieldType * coordinatesField_;
   VectorFieldType * edgesField_;
   VectorFieldType * facesField_;
   ProcIdFieldType * processorIdField_;
   SolutionFieldType * loadBalField_;
   
   // maps field names to solution field stk mesh handles
   std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToSolution_;
   std::map<std::pair<std::string,std::string>,SolutionFieldType*> fieldNameToCellField_;
 
   unsigned dimension_;
 
   bool initialized_;
 
   // how many elements, faces, edges, and nodes are there globally
   std::vector<std::size_t> entityCounts_;
 
   // what is maximum entity ID
   std::vector<stk::mesh::EntityId> maxEntityId_;
 
   unsigned procRank_;
   std::size_t currentLocalId_;
 
   Teuchos::RCP<Teuchos::MpiComm<int> > mpiComm_;
 
   double initialStateTime_; // the time stamp at the time this object was constructed (default 0.0)
   double currentStateTime_; // the time stamp set by the user when writeToExodus is called (default 0.0)
 
#ifdef PANZER_HAVE_IOSS
   // I/O support
   Teuchos::RCP<stk::io::StkMeshIoBroker> meshData_;
   int meshIndex_;
 
  enum class GlobalVariable
  {
    ADD,
    WRITE
  }; // end of enum class GlobalVariable
 
  void
  globalToExodus(
    const GlobalVariable& flag);
 
  Teuchos::ParameterList globalData_;
#endif
 
   // uses lazy evaluation
   mutable Teuchos::RCP<std::vector<stk::mesh::Entity> > orderedElementVector_;
 
   // for element block weights
   std::map<std::string,double> blockWeights_;
 
   std::unordered_map<stk::mesh::EntityId,std::size_t> localIDHash_;
 
   // Store mesh displacement fields by element block. This map
   // goes like this meshCoordFields_[eBlock][axis_index] => coordinate FieldName
   // goes like this meshDispFields_[eBlock][axis_index] => displacement FieldName
   std::map<std::string,std::vector<std::string> > meshCoordFields_; // coordinate  fields written by user
   std::map<std::string,std::vector<std::string> > meshDispFields_;  // displacement fields, output to exodus
 
   bool useFieldCoordinates_;
 
   bool useLowerCase_;
 
   // Object describing how to sort a vector of elements using
   // local ID as the key, very short lived object
   class LocalIdCompare {
   public:
     LocalIdCompare(const STK_Interface * mesh) : mesh_(mesh) {}
   
     // Compares two stk mesh entities based on local ID
     bool operator() (stk::mesh::Entity a,stk::mesh::Entity b) 
     { return mesh_->elementLocalId(a) < mesh_->elementLocalId(b);}
   
   private:
     const STK_Interface * mesh_;
   };
};
 
template <typename ArrayT>
void STK_Interface::setSolutionFieldData(const std::string & fieldName,const std::string & blockId,
                                         const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
{
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   int field_axis = -1;
   if(isMeshCoordField(blockId,fieldName,field_axis)) {
     setDispFieldData(fieldName,blockId,field_axis,localElementIds,solutionValues); 
     return;
   }
 
   SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
 
   for(std::size_t cell=0;cell<localElementIds.size();cell++) {
      std::size_t localId = localElementIds[cell];
      stk::mesh::Entity element = elements[localId];
 
      // loop over nodes set solution values
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      for(std::size_t i=0; i<num_nodes; ++i) {
        stk::mesh::Entity node = nodes[i];
 
        double * solnData = stk::mesh::field_data(*field,node);
        // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
        solnData[0] = scaleValue*solutionValues(cell,i);
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::setDispFieldData(const std::string & fieldName,const std::string & blockId,int axis,
                                     const std::vector<std::size_t> & localElementIds,const ArrayT & dispValues)
{
   TEUCHOS_ASSERT(axis>=0); // sanity check
 
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
   const VectorFieldType & coord_field = this->getCoordinatesField();
 
   for(std::size_t cell=0;cell<localElementIds.size();cell++) {
      std::size_t localId = localElementIds[cell];
      stk::mesh::Entity element = elements[localId];
 
      // loop over nodes set solution values
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      for(std::size_t i=0; i<num_nodes; ++i) {
        stk::mesh::Entity node = nodes[i];
 
        double * solnData = stk::mesh::field_data(*field,node);
        double * coordData = stk::mesh::field_data(coord_field,node);
        // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
        solnData[0] = dispValues(cell,i)-coordData[axis];
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::getSolutionFieldData(const std::string & fieldName,const std::string & blockId,
                                         const std::vector<std::size_t> & localElementIds,ArrayT & solutionValues) const
{
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   solutionValues = Kokkos::createDynRankView(solutionValues,
                "solutionValues",
                localElementIds.size(),
                bulkData_->num_nodes(elements[localElementIds[0]]));
 
   SolutionFieldType * field = this->getSolutionField(fieldName,blockId);
 
   for(std::size_t cell=0;cell<localElementIds.size();cell++) {
      std::size_t localId = localElementIds[cell];
      stk::mesh::Entity element = elements[localId];
 
      // loop over nodes set solution values
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      for(std::size_t i=0; i<num_nodes; ++i) {
        stk::mesh::Entity node = nodes[i];
 
        double * solnData = stk::mesh::field_data(*field,node);
        // TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
        solutionValues(cell,i) = solnData[0];
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::setCellFieldData(const std::string & fieldName,const std::string & blockId,
                                     const std::vector<std::size_t> & localElementIds,const ArrayT & solutionValues,double scaleValue)
{
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   SolutionFieldType * field = this->getCellField(fieldName,blockId);
 
   for(std::size_t cell=0;cell<localElementIds.size();cell++) {
      std::size_t localId = localElementIds[cell];
      stk::mesh::Entity element = elements[localId];
 
      double * solnData = stk::mesh::field_data(*field,element);
      TEUCHOS_ASSERT(solnData!=0); // only needed if blockId is not specified
      solnData[0] = scaleValue*solutionValues(cell,0);
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     //
     // gather from the intrinsic mesh coordinates (non-lagrangian)
     //
 
     const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
    
     // convert to a vector of entity objects
     std::vector<stk::mesh::Entity> selected_elements;
     for(std::size_t cell=0;cell<localElementIds.size();cell++) 
       selected_elements.push_back(elements[localElementIds[cell]]);
 
     getElementVertices_FromCoords(selected_elements,vertices);
   }
   else {
     TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
                                "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
                                "without specifying an element block.");
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoords(elements,vertices);
   }
   else {
     TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
                                "STK_Interface::getElementVertices: Cannot call this method when field coordinates are used "
                                "without specifying an element block.");
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoords(elements,vertices);
   }
   else {
     getElementVertices_FromField(elements,eBlock,vertices);
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
{
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   // convert to a vector of entity objects
   std::vector<stk::mesh::Entity> selected_elements;
   for(std::size_t cell=0;cell<localElementIds.size();cell++) 
     selected_elements.push_back(elements[localElementIds[cell]]);
 
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoords(selected_elements,vertices);
   }
   else {
     getElementVertices_FromField(selected_elements,eBlock,vertices);
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     //
     // gather from the intrinsic mesh coordinates (non-lagrangian)
     //
 
     const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
    
     // convert to a vector of entity objects
     std::vector<stk::mesh::Entity> selected_elements;
     for(std::size_t cell=0;cell<localElementIds.size();cell++) 
       selected_elements.push_back(elements[localElementIds[cell]]);
 
     getElementVertices_FromCoordsNoResize(selected_elements,vertices);
   }
   else {
     TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
                                "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
                                "without specifying an element block.");
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoordsNoResize(elements,vertices);
   }
   else {
     TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
                                "STK_Interface::getElementVerticesNoResize: Cannot call this method when field coordinates are used "
                                "without specifying an element block.");
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVerticesNoResize(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
{
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoordsNoResize(elements,vertices);
   }
   else {
     getElementVertices_FromFieldNoResize(elements,eBlock,vertices);
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVerticesNoResize(const std::vector<std::size_t> & localElementIds,const std::string & eBlock, ArrayT & vertices) const
{
   const std::vector<stk::mesh::Entity> & elements = *(this->getElementsOrderedByLID());
 
   // convert to a vector of entity objects
   std::vector<stk::mesh::Entity> selected_elements;
   for(std::size_t cell=0;cell<localElementIds.size();cell++) 
     selected_elements.push_back(elements[localElementIds[cell]]);
 
   if(!useFieldCoordinates_) {
     getElementVertices_FromCoordsNoResize(selected_elements,vertices);
   }
   else {
     getElementVertices_FromFieldNoResize(selected_elements,eBlock,vertices);
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices_FromCoords(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
{
   // nothing to do! silently return
   if(elements.size()==0) {
     vertices = Kokkos::createDynRankView(vertices,"vertices",0,0,0);
      return;
   }
 
   //
   // gather from the intrinsic mesh coordinates (non-lagrangian)
   //
 
   // get *master* cell toplogy...(belongs to first element)
   unsigned masterVertexCount 
     = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
 
   // allocate space
   vertices = Kokkos::createDynRankView(vertices,"vertices",elements.size(),masterVertexCount,getDimension());
 
   // loop over each requested element
   unsigned dim = getDimension();
   for(std::size_t cell=0;cell<elements.size();cell++) {
      stk::mesh::Entity element = elements[cell];
      TEUCHOS_ASSERT(element!=0);
 
      unsigned vertexCount 
        = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
      TEUCHOS_TEST_FOR_EXCEPTION(vertexCount!=masterVertexCount,std::runtime_error,
                         "In call to STK_Interface::getElementVertices all elements "
                         "must have the same vertex count!");
 
      // loop over all element nodes
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
                         "In call to STK_Interface::getElementVertices cardinality of "
                                 "element node relations must be the vertex count!");
      for(std::size_t node=0; node<num_nodes; ++node) {
        const double * coord = getNodeCoordinates(nodes[node]);
 
        // set each dimension of the coordinate
        for(unsigned d=0;d<dim;d++)
          vertices(cell,node,d) = coord[d];
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices_FromCoordsNoResize(const std::vector<stk::mesh::Entity> & elements, ArrayT & vertices) const
{
   // nothing to do! silently return
   if(elements.size()==0) {
      return;
   }
 
   //
   // gather from the intrinsic mesh coordinates (non-lagrangian)
   //
 
   // get *master* cell toplogy...(belongs to first element)
   unsigned masterVertexCount 
     = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
 
   // loop over each requested element
   unsigned dim = getDimension();
   for(std::size_t cell=0;cell<elements.size();cell++) {
      stk::mesh::Entity element = elements[cell];
      TEUCHOS_ASSERT(element!=0);
 
      unsigned vertexCount 
        = stk::mesh::get_cell_topology(bulkData_->bucket(element).topology()).getCellTopologyData()->vertex_count;
      TEUCHOS_TEST_FOR_EXCEPTION(vertexCount!=masterVertexCount,std::runtime_error,
                         "In call to STK_Interface::getElementVertices all elements "
                         "must have the same vertex count!");
 
      // loop over all element nodes
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      TEUCHOS_TEST_FOR_EXCEPTION(num_nodes!=masterVertexCount,std::runtime_error,
                         "In call to STK_Interface::getElementVertices cardinality of "
                                 "element node relations must be the vertex count!");
      for(std::size_t node=0; node<num_nodes; ++node) {
        const double * coord = getNodeCoordinates(nodes[node]);
 
        // set each dimension of the coordinate
        for(unsigned d=0;d<dim;d++)
          vertices(cell,node,d) = coord[d];
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices_FromField(const std::vector<stk::mesh::Entity> & elements,const std::string & eBlock, ArrayT & vertices) const
{
   TEUCHOS_ASSERT(useFieldCoordinates_);
 
   // nothing to do! silently return
   if(elements.size()==0) {
     vertices = Kokkos::createDynRankView(vertices,"vertices",0,0,0);
      return;
   }
 
   // get *master* cell toplogy...(belongs to first element)
   unsigned masterVertexCount 
     = stk::mesh::get_cell_topology(bulkData_->bucket(elements[0]).topology()).getCellTopologyData()->vertex_count;
 
   // allocate space
   vertices = Kokkos::createDynRankView(vertices,"vertices",elements.size(),masterVertexCount,getDimension());
 
   std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
   if(itr==meshCoordFields_.end()) {
     // no coordinate field set for this element block
     TEUCHOS_ASSERT(false);
   }
 
   const std::vector<std::string> & coordField = itr->second;
   std::vector<SolutionFieldType*> fields(getDimension());
   for(std::size_t d=0;d<fields.size();d++) {
     fields[d] = this->getSolutionField(coordField[d],eBlock);
   }
 
   for(std::size_t cell=0;cell<elements.size();cell++) {
      stk::mesh::Entity element = elements[cell];
 
      // loop over nodes set solution values
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      for(std::size_t i=0; i<num_nodes; ++i) {
        stk::mesh::Entity node = nodes[i];
 
        const double * coord = getNodeCoordinates(node);
 
        for(unsigned d=0;d<getDimension();d++) {
          double * solnData = stk::mesh::field_data(*fields[d],node);
 
          // recall mesh field coordinates are stored as displacements
          // from the mesh coordinates, make sure to add them together
          vertices(cell,i,d) = solnData[0]+coord[d];
        }
      }
   }
}
 
template <typename ArrayT>
void STK_Interface::getElementVertices_FromFieldNoResize(const std::vector<stk::mesh::Entity> & elements,
                                                         const std::string & eBlock, ArrayT & vertices) const
{
   TEUCHOS_ASSERT(useFieldCoordinates_);
 
   // nothing to do! silently return
   if(elements.size()==0) {
      return;
   }
 
   std::map<std::string,std::vector<std::string> >::const_iterator itr = meshCoordFields_.find(eBlock);
   if(itr==meshCoordFields_.end()) {
     // no coordinate field set for this element block
     TEUCHOS_ASSERT(false);
   }
 
   const std::vector<std::string> & coordField = itr->second;
   std::vector<SolutionFieldType*> fields(getDimension());
   for(std::size_t d=0;d<fields.size();d++) {
     fields[d] = this->getSolutionField(coordField[d],eBlock);
   }
 
   for(std::size_t cell=0;cell<elements.size();cell++) {
      stk::mesh::Entity element = elements[cell];
 
      // loop over nodes set solution values
      const size_t num_nodes = bulkData_->num_nodes(element);
      stk::mesh::Entity const* nodes = bulkData_->begin_nodes(element);
      for(std::size_t i=0; i<num_nodes; ++i) {
        stk::mesh::Entity node = nodes[i];
 
        const double * coord = getNodeCoordinates(node);
 
        for(unsigned d=0;d<getDimension();d++) {
          double * solnData = stk::mesh::field_data(*fields[d],node);
 
          // recall mesh field coordinates are stored as displacements
          // from the mesh coordinates, make sure to add them together
          vertices(cell,i,d) = solnData[0]+coord[d];
        }
      }
   }
}
 
}
 
#endif