Grid.H

Go to the documentation of this file.

#ifndef __GRID_H__
#define __GRID_H__

#include "PCPPUtil.H"
#include "PCPPCommUtil.H"
#include "Stencil.H"

#ifdef USE_OMP
#include <omp.h>
#endif

using namespace pcpp::util;
using namespace fixtures;

namespace simulation{
  namespace grid {

    enum TopoType { CARTESIAN=0, UNIRECT, RECTILINEAR, CURVILINEAR, NUMTOPOTYPE };
    static const int NUMGRIDTYPES = NUMTOPOTYPE;
    static const char *topoNames[] = {"Cartesian","Uniform-rectangular","Rectilinear","Curvilinear","Unknown"};
    int ResolveTopoName(const std::string &inName);
    std::string ResolveTopoType(int inType);



    struct subregion {
      std::string regionName;
      int normalDirection;
      pcpp::IndexIntervalType regionInterval;
      pcpp::IndexIntervalType regionPartitionInterval;
      pcpp::IndexIntervalType regionPartitionBufferInterval;
      std::vector<size_t> partitionBufferIndices;
      // Threading support
      std::vector<pcpp::IndexIntervalType> threadPartitionIntervals;
      std::vector<pcpp::IndexIntervalType> threadPartitionBufferIntervals;
      std::vector<std::vector<size_t> >    threadPartitionBufferIndices;
    };


    int InitSubRegionFromString(const std::string &inString,const std::vector<size_t> &gridSizes,
                                subregion &inSubRegion);

    pcpp::IndexIntervalType ResolveBufferInterval(const pcpp::IndexIntervalType &partitionInterval,
                                                  const pcpp::IndexIntervalType &partitionBufferInterval,
                                                  const pcpp::IndexIntervalType &partitionSubInterval);
      
    class halo {
      
    public:

      struct messagebuffers {
        int numComponents;
        int componentSize;
        std::vector<size_t> numPointsSend;
        std::vector<size_t> numPointsRecv;
        std::vector<double *> sendBuffers;
        std::vector<double *> recvBuffers;
        std::vector<int *> sendFlagBuffers;
        std::vector<int *> recvFlagBuffers;
      };

    public:

      halo() : doFill(false), numThreads(1) {};

      halo(const pcpp::IndexIntervalType &inGlobalExtent,
           const pcpp::IndexIntervalType &inPartitionExtent) :
        doFill(false), globalGridExtent(inGlobalExtent), 
        globalPartitionExtent(inPartitionExtent), numThreads(1)
      {};
      
      void SetGridInterval(const pcpp::IndexIntervalType &inExtent) { globalGridExtent = inExtent; };
      void SetPartitionInterval(const pcpp::IndexIntervalType &inExtent) { globalPartitionExtent = inExtent; };

      bool SetFill(bool yesno){ doFill = yesno; return(doFill); };
      bool Fill() {return(doFill); };

      void SetLocalBufferSizes(const std::vector<size_t> &inBufferSizes);

      void SetLocalPartitionExtent(const pcpp::IndexIntervalType &inExtent){ localPartitionExtent = inExtent; };

      std::vector<pcpp::IndexIntervalType> CreateRemoteHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                   const pcpp::IndexIntervalType &partitionExtent,
                                                                   std::vector<size_t> &haloSizes);
      std::vector<pcpp::IndexIntervalType> CreateRemoteHaloExtents(std::vector<size_t> &haloSizes);
      
      std::vector<pcpp::IndexIntervalType> CreateLocalHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                  const pcpp::IndexIntervalType &partitionExtent,
                                                                  std::vector<size_t> &haloSizes);
      std::vector<pcpp::IndexIntervalType> CreateLocalHaloExtents(std::vector<size_t> &haloSizes);
      
      int CreateSimpleSendIndices();
      int CreateSimpleRecvIndices();


      // -------- Threaded interface --------------
      int CreateThreadSendIndices(int threadId);
      int CreateThreadRecvIndices(int threadId);
      int UnpackSimpleRecvBuffers(int threadId);
      //      int PackSimpleSendBuffers(const state::base &inX, int threadId);
      void SetRemoteHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents);
      void SetThreadExtent(int myThreadId,pcpp::IndexIntervalType threadInterval);
      void SetNumThreads(int numThreadsIn); 
      int NumThreads() { return (numThreads); }; 
      // ----------------------------------------
  
      void SetLocalHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents);
      int DestroyMessageBuffer(int messageId);
      int CreateMessageBuffers(int numComponents,int componentSize=8);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const double *sourceBuffer);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const int *sourceBuffer);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const double *sourceBuffer,
                             const int threadId);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const int *sourceBuffer,
                             const int threadId);
      int SendMessage(const int messageId,const std::vector<int> &neighborRanks,
                      fixtures::CommunicatorType &inComm);
      int ReceiveMessage(const int messageId,const std::vector<int> &neighborRanks,
                         fixtures::CommunicatorType &inComm);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,double *targetBuffer);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,double *targetBuffer,
                               const int threadId);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,int *targetBuffer);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,int *targetBuffer,
                               const int threadId);
        
      void SetNeighbors(const std::vector<bool> &inNeighbors){
        haveNeighbor = inNeighbors;
      };


      inline pcpp::IndexIntervalType &GlobalExtent() { return(globalGridExtent); };
      inline pcpp::IndexIntervalType &PartitionExtent() { return (globalPartitionExtent); };
      inline std::vector<pcpp::IndexIntervalType> &RemoteHaloExtents() { return(remoteHaloExtents); };
      inline std::vector<pcpp::IndexIntervalType> &LocalHaloExtents() { return(localHaloExtents); };
      const std::vector<pcpp::IndexIntervalType> &RemoteHaloBufferExtents() const { return(remoteHaloBufferExtents); };
      const std::vector<pcpp::IndexIntervalType> &LocalHaloBufferExtents() const { return(localHaloBufferExtents); };

      std::vector<messagebuffers> &CommunicationBuffers() { return(communicationBuffers); };
      const std::vector<std::vector<size_t> > &SendIndices() const {return(sendIndices);};
      const std::vector<std::vector<size_t> > &RecvIndices() const {return(recvIndices);};
      inline std::vector<int> &PeriodicDirs() { return(periodicDirs); };
      void SetPeriodicDirs(const std::vector<int> &inPeriodic)
      { periodicDirs = inPeriodic; };
      void ReportSimpleBufferContents(std::ostream &outStream);

      void DestroyAll();

      ~halo(){
        DestroyAll();
      }
      
    protected:
      
      bool haveRecvData;
      bool haveSendData;
      bool doFill;
      int numThreads;
      size_t numPointsBuffer;

      std::vector<int> periodicDirs;

      // These are the extents of the partition and
      // halos, resp.
      pcpp::IndexIntervalType globalGridExtent;
      pcpp::IndexIntervalType globalPartitionExtent;
      pcpp::IndexIntervalType localPartitionExtent;
      std::vector<size_t> localBufferSizes;

      std::vector<bool> haveNeighbor;
      std::vector<pcpp::IndexIntervalType> remoteHaloExtents;
      std::vector<pcpp::IndexIntervalType> localHaloExtents;
      std::vector<pcpp::IndexIntervalType> localHaloBufferExtents;
      std::vector<pcpp::IndexIntervalType> remoteHaloBufferExtents;
      std::vector<pcpp::IndexIntervalType> sendExtents;

      std::vector<pcpp::IndexIntervalType> threadExtents;
      std::vector<pcpp::IndexIntervalType> threadBufferExtents;

      // For each send required, these indicate the index ranges
      // and a flat array of indices into the state data
      std::vector<pcpp::RemoteCollisionType> sendCollisions;
      std::vector<std::vector<size_t> > sendIndices;
      std::vector<std::vector<std::vector<size_t> > > threadSendIndices;
      std::vector<std::vector<std::vector<size_t> > > threadSendBufferIndices;
      
      // These indicate the receive buffer intervals
      std::vector<pcpp::RemoteCollisionType>            recvCollisions;
      std::vector<std::vector<size_t> >                 recvIndices;
      std::vector<std::pair<int,std::vector<size_t> > > haloRecvIndices;
      // receive intervals for threads
      std::vector<std::vector<std::vector<size_t> > >   threadRecvIndices;
      std::vector<std::vector<std::vector<size_t> > >   threadRecvBufferIndices;
      std::vector<std::vector<std::vector<size_t> > >   threadHaloBufferIndices;


      
      
      // Communication buffers are densely stuffed and
      // packed/unpacked with special functions
      std::vector<messagebuffers> communicationBuffers;
      
    };


    struct partition_info {
      pcpp::IndexIntervalType globalInterval;
      pcpp::IndexIntervalType localInterval;
    };
  
    struct patch_info {
      std::string Name;
      int normalDir;
      pcpp::IndexIntervalType patchInterval;
    };

  
    struct intervalset
    {
      pcpp::IndexIntervalType bufferInterval;
      pcpp::IndexIntervalType opInterval;
      std::vector<pcpp::IndexIntervalType> threadBufferIntervals;
      std::vector<pcpp::IndexIntervalType> threadOpIntervals;
    };

    class base {
    
    protected:
      fixtures::ConfigurationType gridConfig;
      size_t   numPoints;
      int      numDim;

    public:
      virtual int ConfigureGrid(fixtures::ConfigurationType &inConfig,std::ostream &messageStream)
      { 
        /*
          advect1d:Dimension = 3
          advect1d:GridType = UniRect
          advect1d:PhysicalDomain = 0 1 0 1 0 1    # Cube domain [0,1]^3
          advect1d:NumPoints = 24000000 3 3        # vertices in each dir
          advect1d:GenerateGrid = Yes              # do grid generation

        */

        std::string configName(inConfig.GetValue("ConfigName"));
        fixtures::ConfigurationType::ParamType gridParam;
        std::ostringstream paramNameOut;
        paramNameOut << "Dimension " << "GridType " 
                     << "PhysicalExtent " << "NumPoints "
                     << "GenerateGrid " << "PhysicalDomain "
                     << "DecompDirs " << "ThreadDecompDirs";
        std::string configParamNames(paramNameOut.str());
        std::istringstream paramNameIn(configParamNames);
        std::string paramName;
        while(paramNameIn >> paramName){
          std::string configKey(ConfigKey(configName,paramName));
          if(inConfig.IsSet(configKey)){
            gridParam.first = paramName;
            gridParam.second = inConfig.GetValue(configKey);
            gridConfig.push_back(gridParam);
          }
        }
        return(0);
      };

      virtual int Configure(fixtures::ConfigurationType &inConfig)
      { 
        /*
          advect1d:Dimension = 3
          advect1d:GridType = UniRect
          advect1d:PhysicalDomain = 0 1 0 1 0 1    # Cube domain [0,1]^3
          advect1d:NumPoints = 24000000 3 3        # vertices in each dir
          advect1d:GenerateGrid = Yes              # do grid generation

        */

        std::string configName(inConfig.GetValue("ConfigName"));
        fixtures::ConfigurationType::ParamType gridParam;
        std::ostringstream paramNameOut;
        paramNameOut << "Dimension " << "GridType " 
                     << "PhysicalDomain " << "NumPoints "
                     << "GenerateGrid";
        std::string configParamNames(paramNameOut.str());
        std::istringstream paramNameIn(configParamNames);
        std::string paramName;
        while(paramNameIn >> paramName){
          std::string configKey(ConfigKey(configName,paramName));
          if(inConfig.IsSet(configKey)){
            gridParam.first = paramName;
            gridParam.second = inConfig.GetValue(configKey);
            gridConfig.push_back(gridParam);
          }
        }
        return(0);
      };

      virtual std::string ReportConfiguration()
      {
        std::ostringstream Ostr;
        Ostr << gridConfig << std::endl;
        return(Ostr.str());
      };
      int NumDim() const { return(numDim); };
      size_t NumPoints() const { return(numPoints); };
      int Dimension() const { return(numDim); };
    };
  
    class parallel_base : public base {
    
    public:
    
      //    using base::ConfigureGrid;

      virtual int ConfigureGrid(fixtures::ConfigurationType &inConfig,std::ostream &messageStream) {
        //        gridCommunicator.Initialize(inCommunicator);
        return(base::ConfigureGrid(inConfig,messageStream));
      };
      virtual int Configure(fixtures::ConfigurationType &inConfig){
        //        gridCommunicator.Initialize(inCommunicator);
        return(base::Configure(inConfig));
      };

      virtual std::string ReportConfiguration()
      {
        //       std::ostringstream Ostr;
        //       Ostr << base::ReportConfiguration() 
        //         << "grid::parallel_base: ReportConfig" << std::endl;
        //       return(Ostr.str());
        return(base::ReportConfiguration());
      };
      
      fixtures::CommunicatorType &Communicator() const 
      { return(const_cast<fixtures::CommunicatorType &>(gridCommunicator)); };
      
      size_t NumGlobalPoints() { return(numGlobalPoints); };
      size_t NumLocalPoints() { return(numLocalPoints); };
    protected:    
      fixtures::CommunicatorType gridCommunicator;
      size_t numGlobalPoints;
      size_t numLocalPoints;
    };
  

    struct pbsintervals {
      pbsintervals() : numDim(0) {};
      int numDim;
      std::vector<size_t> gridSizes;
      std::vector<size_t> bufferSizes;
      std::vector<size_t> partitionSizes;
      pcpp::IndexIntervalType partitionInterval;
      pcpp::IndexIntervalType partitionBufferInterval;
      //#ifdef USE_OMP
      std::vector<pcpp::IndexIntervalType> threadBufferIntervals;
      std::vector<pcpp::IndexIntervalType> threadPartitionBufferIntervals;
      std::vector<pcpp::IndexIntervalType> threadPartitionIntervals;
      std::vector<int> decompSizes;
      //#endif
    };

    class parallel_blockstructured : public parallel_base
    {
      
    public:

      typedef halo HaloType;

    public:

      int SetupThreads(std::vector<int> &numThreadPartitions);

      parallel_blockstructured() : parallel_base(), haloPtr(&gridHalo), 
                                   stencilSet(NULL), stencilConnectivity(NULL),
                                   xyzMessageId(-1), metricMessageId(-1),
                                   gridType(UNIRECT)
      {};

      
      int NumberOfThreads() const 
      { return(partitionIntervals.threadPartitionBufferIntervals.size()); };

      void SetNumThreads(int numThreads);

      std::vector<pcpp::IndexIntervalType> &ThreadPartitionIntervals()
      { return(partitionIntervals.threadPartitionIntervals); };

      const std::vector<pcpp::IndexIntervalType> &ThreadPartitionIntervals() const
      { return(partitionIntervals.threadPartitionIntervals); };

      std::vector<pcpp::IndexIntervalType> &ThreadPartitionBufferIntervals()
      { return(partitionIntervals.threadPartitionBufferIntervals); };

      const std::vector<pcpp::IndexIntervalType> &ThreadPartitionBufferIntervals() const
      { return(partitionIntervals.threadPartitionBufferIntervals); };

      std::vector<pcpp::IndexIntervalType> &ThreadBufferIntervals()
      { return(partitionIntervals.threadBufferIntervals); };

      const std::vector<pcpp::IndexIntervalType> &ThreadBufferIntervals() const
      { return(partitionIntervals.threadBufferIntervals); };

      int Dimension() const { return (partitionIntervals.numDim); };

      size_t BufferSize() const { return(numPointsBuffer); };

      size_t PartitionSize() const { return(numPointsPartition); };
      
      pbsintervals &Intervals(){ return(partitionIntervals); };
      const pbsintervals &Intervals() const { return(partitionIntervals); };

      void SetDecompSizes(const std::vector<int> &inDecomp)
      {
        partitionIntervals.decompSizes = inDecomp;
      };

      const std::vector<int> &DecompSizes() const
      {
        return(partitionIntervals.decompSizes);
      };

      std::vector<int> &DecompSizes()
      {
        return(partitionIntervals.decompSizes);
      };
      
      void SetGridSizes(const std::vector<size_t> &inSize)
      {
        partitionIntervals.gridSizes = inSize;
        partitionIntervals.numDim = inSize.size();
        numDim = partitionIntervals.numDim;
      };

      const std::vector<size_t> &BufferSizes() const
      {
        return(partitionIntervals.bufferSizes);
      };

      std::vector<size_t> &BufferSizes()
      {
        return(partitionIntervals.bufferSizes);
      };

      const std::vector<size_t> &GridSizes() const
      {
        return(partitionIntervals.gridSizes);
      };

      std::vector<size_t> &GridSizes() 
      { 
        return(partitionIntervals.gridSizes);
      };


      void SetPartitionSizes(const std::vector<size_t> &inPartSize)
      {
        partitionIntervals.partitionSizes = inPartSize;
        numDim = inPartSize.size();
        partitionIntervals.numDim = numDim;
      };

      const std::vector<size_t> &PartitionSizes() const
      {
        return(partitionIntervals.partitionSizes);
      };

      void SetDimensionExtensions(const std::vector<size_t> &inExtensions)
      {
        dimExtensions = inExtensions;
      };

      void SetPartitionInterval(const pcpp::IndexIntervalType &inInterval)
      { partitionIntervals.partitionInterval = inInterval; };

      pcpp::IndexIntervalType &PartitionInterval(){ return(partitionIntervals.partitionInterval); };
      const pcpp::IndexIntervalType &PartitionInterval() const { return(partitionIntervals.partitionInterval); };

      std::vector<double> &CoordinateData(){ return(gridCoordinates); };
      const std::vector<double> &CoordinateData() const { return(gridCoordinates); };

      const double *CoordinateData(int inDim) const
      {
        assert(inDim < numDim);
        return(&gridCoordinates[inDim*numPointsBuffer]);
      };
      double *CoordinateData(int inDim)
      {
        assert(inDim < numDim);
        return(&gridCoordinates[inDim*numPointsBuffer]);
      };

      int ExchangeCoordinates();

      pcpp::IndexIntervalType &PartitionBufferInterval(){ return(partitionIntervals.partitionBufferInterval); };
      const pcpp::IndexIntervalType &PartitionBufferInterval() const { return(partitionIntervals.partitionBufferInterval); };

      int Finalize(bool allocateCoordinateData = false);

      void AllocateCoordinateData()
      {
        size_t numAllocate = numPointsBuffer*numDim;
        if(numAllocate > 0){
          if(gridCoordinates.size() != numAllocate)
            gridCoordinates.resize(numAllocate);
        }
      };

      const std::vector<double> &PhysicalExtent() const {return (physicalExtent); };
      std::vector<double> &PhysicalExtent() {return (physicalExtent); };
      void SetPhysicalExtent(const std::vector<double> &inExtent){
        physicalExtent = inExtent;
      };
      const std::vector<double> &DX() const {return(dX);};
      std::vector<double> &DX() {return(dX);};
      void SetGridSpacings(const std::vector<double> &inSpacing){
        dX = inSpacing;
      };
      void SetHalo(halo &inHalo) { haloPtr = &inHalo;}
      HaloType &Halo() { return(*haloPtr); };
      int GenerateCoordinates(std::ostream &);
      int GenerateCoordinates(int (*)(const std::vector<size_t> &,std::vector<double> &),
                              std::ostream &);
      int GenerateGrid(int (*)(const std::vector<int> &,
                               const std::vector<double> &,
                               const std::vector<size_t> &,
                               const std::vector<size_t> &,
                               std::vector<double> &),
                       std::ostream &);
      
      const std::vector<subregion> &SubRegions() const
      { return(subRegions); };
      
      std::vector<subregion> &SubRegions()
      { return(subRegions); };
      
      int GetRegionIndex(const std::string &regionName) const {
        int regionIndex = 0;
        std::vector<subregion>::const_iterator regIt = subRegions.begin();
        while(regIt != subRegions.end()){
          if(regIt++->regionName == regionName)
            return(regionIndex);
          regionIndex++;
        }
        return(-1);
      };

      // Dumb interface - deprecated
      //    pcpp::IndexIntervalType &LocalPartitionInterval(){ return(localPartInterval); };
      //    const pcpp::IndexIntervalType &LocalPartitionInterval() const { return(localPartInterval); };

      std::vector<bool> &PeriodicDirs()
      { return(periodicDirs); };
      const std::vector<bool> &PeriodicDirs() const
      { return(periodicDirs); };
      std::vector<int> &DecompDirs()
      { return(decompDirs); };
      const std::vector<int> &DecompDirs() const
      { return(decompDirs); };
      void SetDecompDirs(const std::vector<int> &inDecomp)
      { decompDirs = inDecomp; };
      void SetPeriodicDirs(const std::vector<bool> &inPeriodic)
      { periodicDirs = inPeriodic; };
      void SetPeriodicLengths(const std::vector<double> &inLengths)
      { periodicLengths = inLengths; };
      std::vector<int> &ThreadDecompDirs()
      { return(threadDecompDirs); };
      const std::vector<int> &ThreadDecompDirs() const 
      { return(threadDecompDirs); };
      void SetThreadDecompDirs(const std::vector<int> &inThread)
      { threadDecompDirs = inThread; };
      void SetType(int inType)
      {
        gridType = inType; 
      };
      int Type() const
      { return(gridType); };
      void SetScaleFactor(const std::vector<double> &inScale)
      { scaleFactor = inScale; };
      int ParallelSetup(fixtures::CommunicatorType &inComm,pcpp::ParallelTopologyInfoType &cartInfo,
                        int haloDepth,std::ostream &messageStream);

      std::vector<int>    &FlagData() 
      {
        size_t numFlags = BufferSize();
        if(flagData.size() != numFlags)
          flagData.resize(numFlags,0);
        return(flagData);
      };

      const std::vector<int>    &FlagData() const
      {
        return(flagData);
      };

      void DestroyFlags() { std::vector<int>().swap(flagData); };
      size_t FlagMask(int *iMask,int flagValue = 0,int flagBit = 1) const; 
      std::vector<double> &Metric() { return(gridMetric); };
      std::vector<double> Metric() const { return(gridMetric); };
      std::vector<double> &Jacobian() {return(gridJacobian); };
      std::vector<double> Jacobian() const { return(gridJacobian); };      
      std::vector<double> &JacobianMatrix() { return(jacobianMatrix); };
      std::vector<double> JacobianMatrix() const { return(jacobianMatrix); };

      int ComputeUniformRectangularMetrics(std::ostream &messageStream);
      int ComputeRectilinearMetrics(std::ostream &messageStream);
      int ComputeCurvilinearMetrics2D(std::ostream &messageStream);
      int ComputeCurvilinearMetrics3D(std::ostream &messageStream);
      int CurvilinearMetrics(std::ostream &messageStream);
      int ComputeMetrics(std::ostream &messageStream);
      int ComputeMetricIdentities(std::vector<double> &identityData,std::ostream &messageStream);
      int SetupDifferentialOperator(plascom2::operators::stencilset *inStencilSet,int *inStencilConnectivity);
      int ExchangeNodalData(int numComponents,double *dataBuffer);
      int ExchangeJacobianMatrix(std::ostream &messageStream);
      int GradIJK(int numComponents,double *sourceBuffer,
                  double *targetBuffer,std::ostream &messageStream,bool exchange);
      int ComputeJacobianMatrix(std::ostream &messageStream);
      int ComputeJacobianMatrix2(std::ostream &messageStream);
      std::vector<double> &GridGenerationParameters()
      { return(gridGenParams); };
      const std::vector<double> &GridGenerationParameters() const
      { return(gridGenParams); };
      void SetGridGenerationParameters(const std::vector<double> &inParams)
      { gridGenParams = inParams; };
      std::vector<int> &GridGenerationOptions()
      { return(gridGenOptions); };
      const std::vector<int> &GridGenerationOptions() const
      { return(gridGenOptions); };
      void SetGridGenerationOptions(const std::vector<int> &inOptions)
      { gridGenOptions = inOptions; };

    protected:

      int gridType;
      size_t numPointsBuffer;
      size_t numPointsPartition;
      pbsintervals partitionIntervals;
      std::vector<size_t> dimExtensions;
      pcpp::IndexIntervalType globalBufferInterval;
      std::vector<double> gridCoordinates;
      TopoType gridTopoType;
      std::vector<double> gridMetric;
      std::vector<double> gridJacobian;
      std::vector<double> jacobianMatrix;
      std::vector<double> jacobianMatrix2;
      std::vector<double> gridGenParams;
      std::vector<int>    gridGenOptions;
      std::vector<int>    flagData;

      std::vector<double> physicalExtent;
      std::vector<double> dX;
      std::vector<double> scaleFactor;
      plascom2::operators::stencilset *stencilSet;
      int *stencilConnectivity;

      halo gridHalo;
      halo *haloPtr;
      std::vector<bool> periodicDirs;
      std::vector<double> periodicLengths;
      std::vector<int>  decompDirs;
      std::vector<int>  threadDecompDirs;
      int xyzMessageId;
      int metricMessageId;

      std::vector<subregion> subRegions;

      int errorState;
      
    };



    class unstructured : public base 
    {
    };

    template<typename GridType>
    std::vector<double> GenerateUniformGrid(GridType &inGrid,const std::vector<double> &physicalExtent)
    {
      int numDim = inGrid.Dimension();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      pcpp::IndexIntervalType &partInterval(inGrid.PartitionInterval());
      pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      std::vector<double> dX(numDim,0.0);
      size_t numPointsBuffer = inGrid.BufferSize();
      for(int iDim = 0;iDim < numDim;iDim++)
        dX[iDim] = (physicalExtent[2*iDim+1] - physicalExtent[2*iDim])/static_cast<double>(gridSizes[iDim]-1);
      inGrid.AllocateCoordinateData();
      std::vector<double> &gridCoordinates(inGrid.CoordinateData());
    
      if(numDim == 3){
      
        size_t bufferPlane = bufferSizes[0]*bufferSizes[1];
        size_t bufferLine  = bufferSizes[0];
      
        size_t xStart = partInterval[0].first;
        size_t xEnd   = partInterval[0].second;
        size_t yStart = partInterval[1].first;
        size_t yEnd   = partInterval[1].second;
        size_t zStart = partInterval[2].first;
        size_t zEnd   = partInterval[2].second;
      
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        size_t kStart = partitionBufferInterval[2].first;
        size_t kEnd   = partitionBufferInterval[2].second;
      
        size_t kIndex = kStart;
        size_t zIndex = zStart;
      
        for(kIndex = kStart,zIndex = zStart;kIndex <= kEnd;kIndex++,zIndex++){
          size_t jIndex = jStart;
          size_t yIndex = yStart;
          size_t bzIndex = kIndex * (bufferPlane);
          double zCoordinate = zIndex*dX[2] + physicalExtent[4];
          for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
            size_t iIndex = iStart;
            size_t xIndex = xStart;
            size_t byzIndex = bzIndex + jIndex*bufferLine; 
            double yCoordinate = yIndex*dX[1] + physicalExtent[2];
            for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
              size_t bufferIndex = byzIndex + iIndex;
              double xCoordinate = xIndex*dX[0] + physicalExtent[0];
              gridCoordinates[bufferIndex] = xCoordinate;
              gridCoordinates[bufferIndex+numPointsBuffer] = yCoordinate;
              gridCoordinates[bufferIndex+2*numPointsBuffer] = zCoordinate;
            }
          }
        }
      } else if (numDim == 2) {
      
        size_t bufferLine  = bufferSizes[0];
      
        size_t xStart = partInterval[0].first;
        size_t xEnd   = partInterval[0].second;
        size_t yStart = partInterval[1].first;
        size_t yEnd   = partInterval[1].second;
      
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
      
        size_t jIndex = jStart;
        size_t yIndex = yStart;
        for(jIndex = jStart,yIndex = yStart;jIndex <= jEnd;jIndex++,yIndex++){
          size_t iIndex = iStart;
          size_t xIndex = xStart;
          size_t byIndex = jIndex*bufferLine; 
          double yCoordinate = yIndex*dX[1] + physicalExtent[2];
          for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
            size_t bufferIndex = byIndex + iIndex;
            double xCoordinate = xIndex*dX[0] + physicalExtent[0];
            gridCoordinates[bufferIndex] = xCoordinate;
            gridCoordinates[bufferIndex+numPointsBuffer] = yCoordinate;
          }
        }
      } else if (numDim == 1) {
      
      
        size_t xStart = partInterval[0].first;
        size_t xEnd   = partInterval[0].second;
      
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
      
        size_t iIndex = iStart;
        size_t xIndex = xStart;
        for(iIndex = iStart,xIndex = xStart;iIndex <= iEnd;iIndex++,xIndex++){
          gridCoordinates[iIndex] = xIndex*dX[0] + physicalExtent[0];
        }
      }
      return(dX);
    }

  };
};
#endif