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