/*=========================================================================

Copyright (c) 2007, Los Alamos National Security, LLC

All rights reserved.

Copyright 2007. Los Alamos National Security, LLC.
This software was produced under U.S. Government contract DE-AC52-06NA25396
for Los Alamos National Laboratory (LANL), which is operated by
Los Alamos National Security, LLC for the U.S. Department of Energy.
The U.S. Government has rights to use, reproduce, and distribute this software.
NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY,
EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
If software is modified to produce derivative works, such modified software
should be clearly marked, so as not to confuse it with the version available
from LANL.

Additionally, redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following conditions
are met:
-   Redistributions of source code must retain the above copyright notice,
    this list of conditions and the following disclaimer.
-   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.
-   Neither the name of Los Alamos National Security, LLC, Los Alamos National
    Laboratory, LANL, the U.S. Government, nor the names of its contributors
    may be used to endorse or promote products derived from this software
    without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS
"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 LOS ALAMOS NATIONAL SECURITY, LLC OR
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.

=========================================================================*/

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>

#include <stddef.h>

#include <sys/types.h>

#ifdef _WIN32
#include "winDirent.h"
#else
#include <dirent.h>
#endif

#include "Partition.h"
#include "ParticleDistribute.h"

#ifdef USE_VTK_COSMO
#include "vtkStdString.h"
#include "vtkSetGet.h"
#endif

#include <cstring>
using namespace std;

/////////////////////////////////////////////////////////////////////////
//
// Particle data space is partitioned for the number of processors
// which currently is a factor of two but is easily extended.  Particles
// are read in from files where each processor reads one file into a buffer,
// extracts the particles which really belong on the processor (ALIVE) and
// those in a buffer region around the edge (DEAD).  The buffer is then
// passed round robin to every other processor so that all particles are
// examined by all processors.  All dead particles are tagged with the
// neighbor zone (26 neighbors in 3D) so that later halos can be associated
// with zones.
//
/////////////////////////////////////////////////////////////////////////

ParticleDistribute::ParticleDistribute()
{
  // Get the number of processors running this problem and rank
  this->numProc = Partition::getNumProc();
  this->myProc = Partition::getMyProc();

  // Get the number of processors in each dimension
  Partition::getDecompSize(this->layoutSize);

  // Get my position within the Cartesian topology
  Partition::getMyPosition(this->layoutPos);

  // Get neighbors of this processor including the wraparound
  Partition::getNeighbors(this->neighbor);

  this->numberOfAliveParticles = 0;
  this->massConvertFactor = 1.0;
  this->distConvertFactor = 1.0;
}

ParticleDistribute::~ParticleDistribute()
{
}

/////////////////////////////////////////////////////////////////////////
//
// Set parameters for particle distribution
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::setParameters(
                        const string& baseName,
                        POSVEL_T rL,
                        string dataType)
{
  // Base file name which will have processor id appended for actual files
  this->baseFile = baseName;

  // Physical total space and amount of physical space to use for dead particles
  this->boxSize = rL;

  // RECORD format is the binary .cosmo of one particle with all information
  if (dataType == "RECORD")
    this->inputType = RECORD;

  // BLOCK format is Gadget format with a header and x,y,z locations for
  // all particles, then x,y,z velocities for all particles, and all tags
  else if (dataType == "BLOCK")
    this->inputType = BLOCK;

#ifndef USE_VTK_COSMO
  if (this->myProc == MASTER) {
    cout << endl << "------------------------------------" << endl;
    cout << "boxSize:  " << this->boxSize << endl;
  }
#endif
}

/////////////////////////////////////////////////////////////////////////
//
// Set parameters for particle unit conversion
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::setConvertParameters(
                        POSVEL_T massFactor,
                        POSVEL_T distFactor)
{
  this->massConvertFactor = massFactor;
  this->distConvertFactor = distFactor;
}

/////////////////////////////////////////////////////////////////////////
//
// Set box sizes for determining if a particle is in the alive or dead
// region of this processor.  Data space is a DIMENSION torus.
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::initialize()
{
#ifndef USE_VTK_COSMO
#ifdef DEBUG
  if (this->myProc == MASTER)
    cout << "Decomposition: [" << this->layoutSize[0] << ":"
         << this->layoutSize[1] << ":" << this->layoutSize[2] << "]" << endl;
#endif
#endif

  // Set subextents on particle locations for this processor
  POSVEL_T boxStep[DIMENSION];
  for (int dim = 0; dim < DIMENSION; dim++) {
    boxStep[dim] = this->boxSize / this->layoutSize[dim];

    // Alive particles
    this->minAlive[dim] = this->layoutPos[dim] * boxStep[dim];
    this->maxAlive[dim] = this->minAlive[dim] + boxStep[dim];
    if (this->maxAlive[dim] > this->boxSize)
      this->maxAlive[dim] = this->boxSize;
  }
}


void ParticleDistribute::setParticles(vector<POSVEL_T>* xLoc,
                                      vector<POSVEL_T>* yLoc,
                                      vector<POSVEL_T>* zLoc,
                                      vector<POSVEL_T>* xVel,
                                      vector<POSVEL_T>* yVel,
                                      vector<POSVEL_T>* zVel,
                                      vector<POSVEL_T>* mass,
                                      vector<ID_T>* id)
{
  this->xx = xLoc;
  this->yy = yLoc;
  this->zz = zLoc;
  this->vx = xVel;
  this->vy = yVel;
  this->vz = zVel;
  this->ms = mass;
  this->tag = id;
}

#ifndef USE_SERIAL_COSMO
/////////////////////////////////////////////////////////////////////////
//
// Each processor reads 0 or more files, a buffer at a time, and shares
// the particles using an all-to-all with every other processor
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readParticlesAllToAll(int reserveQ)
{
  // Find how many input files there are and deal them between the processors
  // Calculates the max number of files per processor and max number of
  // particles per file so that buffering can be done
  // For round robin sharing determine where to send and receive buffers from
  partitionInputFiles(true);

  // Compute the total number of particles in the problem
  // Compute the maximum number of particles in any one file to set buffer size
  findFileParticleCount();

  // If there is only one input file we don't have to do MPI messaging
  // because each processor will read that same file and extract only
  // the particles in range
  if (this->numberOfFiles == 1) {
    if (this->inputType == RECORD) {
      readFromRecordFile();
    } else {
      readFromBlockFile();
    }
  } else {

  // MPI buffer size might limit the number of particles read from a file
  // and passed round robin
  // Largest file will have a number of buffer chunks to send if it is too large
  // Every processor must send that number of chunks even if its own file
  // does not have that much information

  if (ENFORCE_MAX_READ == true && this->maxParticles > MAX_READ) {
    this->maxRead = MAX_READ;
    this->maxReadsPerFile = (this->maxParticles / this->maxRead) + 1;
  } else {
    this->maxRead = this->maxParticles;
    this->maxReadsPerFile = 1;
  }

  MPI_Datatype CosmoParticleType;

  {
    MPI_Datatype type[3] = { MPI_FLOAT, MPI_INT, MPI_UB };
    int blocklen[3] = { COSMO_FLOAT, COSMO_INT, 1 };
    MPI_Aint disp[3] = { offsetof(CosmoParticle, floatData),
                         offsetof(CosmoParticle, intData),
                         sizeof(CosmoParticle) };
    MPI_Type_struct(3, blocklen, disp, type, &CosmoParticleType);
    MPI_Type_commit(&CosmoParticleType);
  }

  int numProc = Partition::getNumProc();
  std::vector<int> pCounts(numProc), pRecvCounts(numProc),
                   pDisp(numProc),   pRecvDisp(numProc);
  std::vector< std::vector<CosmoParticle> > pByProc(numProc);
  std::vector<CosmoParticle> pBuffer, pRecvBuffer;

  // Allocate space to hold buffer information for reading of files
  // Mass is constant use that float to store the tag
  // Number of particles is the first integer in the buffer

  // Allocate space for the data read from the file
  POSVEL_T *fBlock = 0;
  POSVEL_T *lBlock = 0;
  POSVEL_T *vBlock = 0;
  ID_T* iBlock = 0;

  // RECORD format reads one particle at a time
  if (this->inputType == RECORD) {
    fBlock = new POSVEL_T[COSMO_FLOAT];
    iBlock = new ID_T[COSMO_INT];
  }

  // BLOCK format reads all particles at one time for triples
  else if (this->inputType == BLOCK) {
    lBlock = new POSVEL_T[this->maxRead * DIMENSION];
    vBlock = new POSVEL_T[this->maxRead * DIMENSION];
    iBlock = new ID_T[this->maxRead];
  }

  // Reserve particle storage to minimize reallocation
  int reserveSize = (int) (this->maxFiles * this->maxParticles * DEAD_FACTOR);

  // If multiple processors are reading the same file we can reduce size
  reserveSize /= this->processorsPerFile;

  if(reserveQ) {
#ifndef USE_VTK_COSMO
    cout << "readParticlesRoundRobin reserving vectors" << endl;
#endif
    this->xx->reserve(reserveSize);
    this->yy->reserve(reserveSize);
    this->zz->reserve(reserveSize);
    this->vx->reserve(reserveSize);
    this->vy->reserve(reserveSize);
    this->vz->reserve(reserveSize);
    this->ms->reserve(reserveSize);
    this->tag->reserve(reserveSize);
  }

  // Running total and index into particle data on this processor
  this->particleCount = 0;

  // Using the input files assigned to this processor, read the input
  // and push all-to-all to every other processor
  // this->maxFiles is the maximum number to read on any processor
  // Some processors may have no files to read but must still participate
  // in the round robin distribution

  for (int file = 0; file < this->maxFiles; file++) {

    // Open file to read the data if any for this processor
    ifstream* inStream = 0;
    int firstParticle = 0;
    int numberOfParticles = 0;
    int remainingParticles = 0;

    if ((int)this->inFiles.size() > file) {
      inStream = new ifstream(this->inFiles[file].c_str(), ios::in|ios::binary);

#ifndef USE_VTK_COSMO
      cout << "Rank " << this->myProc << " open file " << inFiles[file]
           << " with " << this->fileParticles[file] << " particles" << endl;
#endif

      // Number of particles read at one time depends on MPI buffer size
      numberOfParticles = this->fileParticles[file];
      if (numberOfParticles > this->maxRead)
        numberOfParticles = this->maxRead;

      // If a file is too large to be passed as an MPI message divide it up
      remainingParticles = this->fileParticles[file];

    } else {
#ifndef USE_VTK_COSMO
      cout << "Rank " << this->myProc << " no file to open " << endl;
#endif
    }

    for (int piece = 0; piece < this->maxReadsPerFile; piece++) {
      // Reset the comm buffers
      for (int i = 0; i < numProc; ++i) {
        pByProc[i].clear();
      }

      pBuffer.clear();
      pRecvBuffer.clear();

      // Processor has a file to read and share via round robin with others
      if (file < (int)this->inFiles.size()) {
        if (this->inputType == RECORD) {
          readFromRecordFile(inStream, firstParticle, numberOfParticles,
                             fBlock, iBlock, pByProc);
        } else {
          readFromBlockFile(inStream, firstParticle, numberOfParticles,
                           this->fileParticles[file],
                           lBlock, vBlock, iBlock, pByProc);
        }
        firstParticle += numberOfParticles;
        remainingParticles -= numberOfParticles;
        if (remainingParticles <= 0)
          numberOfParticles = 0;
        else if (remainingParticles < numberOfParticles)
          numberOfParticles = remainingParticles;
      }

      // Record all sizes into a single buffer and send this to all ranks
      long totalToSend = 0;
      for (int i = 0; i < numProc; ++i) {
        int sz = (int) pByProc[i].size();
        pCounts[i] = sz;
        totalToSend += sz;
      }

      MPI_Alltoall(&pCounts[0], 1, MPI_INT,
                   &pRecvCounts[0], 1, MPI_INT,
                   Partition::getComm());

      // pRecvCounts now holds the number of particles that this rank should
      // get from every other rank
      long totalToRecv = 0;
      for (int i = 0; i < numProc; ++i) {
        totalToRecv += pRecvCounts[i];
      }

      // Allocate and pack the buffer with all particles to send
      pBuffer.reserve(totalToSend);
      for (int i = 0; i < numProc; ++i)
      for (int j = 0; j < (int) pByProc[i].size(); ++j) {
        pBuffer.push_back(pByProc[i][j]);
      }

      // Calculate displacements
      pDisp[0] = pRecvDisp[0] = 0;
      for (int i = 1; i < numProc; ++i) {
        pDisp[i] = pDisp[i-1] + pCounts[i-1];
        pRecvDisp[i] = pRecvDisp[i-1] + pRecvCounts[i-1];
      }

      // Send all particles to their new homes
      pRecvBuffer.resize(totalToRecv);
      MPI_Alltoallv(&pBuffer[0], &pCounts[0], &pDisp[0], CosmoParticleType,
                    &pRecvBuffer[0], &pRecvCounts[0], &pRecvDisp[0], CosmoParticleType,
                    Partition::getComm());

      // We now have all of our particles, put them in our local arrays
      for (long i = 0; i < totalToRecv; ++i) {
        POSVEL_T loc[DIMENSION], vel[DIMENSION], mass;
        ID_T id;

        int j = 0;
        for (int dim = 0; dim < DIMENSION; dim++) {
          loc[dim] = pRecvBuffer[i].floatData[j++];
        }
        for (int dim = 0; dim < DIMENSION; dim++) {
          vel[dim] = pRecvBuffer[i].floatData[j++];
        }
        id = pRecvBuffer[i].intData[0];

        // Is the particle ALIVE on this processor
        if ((loc[0] >= minAlive[0] && loc[0] < maxAlive[0]) &&
            (loc[1] >= minAlive[1] && loc[1] < maxAlive[1]) &&
            (loc[2] >= minAlive[2] && loc[2] < maxAlive[2])) {

          this->xx->push_back(loc[0]);
          this->yy->push_back(loc[1]);
          this->zz->push_back(loc[2]);
          this->vx->push_back(vel[0]);
          this->vy->push_back(vel[1]);
          this->vz->push_back(vel[2]);
          this->ms->push_back(mass);
          this->tag->push_back(id);

          this->numberOfAliveParticles++;
          this->particleCount++;
        }
      }
    }

    // Can delete the read buffers as soon as last file is read because
    // information has been transferred into the double buffer1
    if (file == (this->maxFiles - 1)) {
      if (this->inputType == RECORD) {
        delete [] fBlock;
        delete [] iBlock;
      } else if (this->inputType == BLOCK) {
        delete [] lBlock;
        delete [] vBlock;
        delete [] iBlock;
      }
    }

    if ((int)this->inFiles.size() > file)
      inStream->close();
  }

  // Count the particles across processors
  long totalAliveParticles = 0;
  MPI_Allreduce((void*) &this->numberOfAliveParticles,
                (void*) &totalAliveParticles,
                1, MPI_LONG, MPI_SUM, Partition::getComm());

#ifndef USE_VTK_COSMO
#ifdef DEBUG
  cout << "Rank " << setw(3) << this->myProc
       << " #alive = " << this->numberOfAliveParticles << endl;
#endif

  if (this->myProc == MASTER) {
    cout << "TotalAliveParticles " << totalAliveParticles << endl;
  }
#endif

  MPI_Type_free(&CosmoParticleType);

  }
}
#endif // USE_SERIAL_COSMO

/////////////////////////////////////////////////////////////////////////
//
// Each processor reads 0 or more files, a buffer at a time, and shares
// the particles by passing the buffer round robin to every other processor
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readParticlesRoundRobin(int reserveQ)
{
  // Find how many input files there are and deal them between the processors
  // Calculates the max number of files per processor and max number of
  // particles per file so that buffering can be done
  // For round robin sharing determine where to send and receive buffers from
  partitionInputFiles();

  // Compute the total number of particles in the problem
  // Compute the maximum number of particles in any one file to set buffer size
  findFileParticleCount();

  // If there is only one input file we don't have to do MPI messaging
  // because each processor will read that same file and extract only
  // the particles in range
  if (this->numberOfFiles == 1) {
    if (this->inputType == RECORD) {
      readFromRecordFile();
    } else {
      readFromBlockFile();
    }
  } else {

  // MPI buffer size might limit the number of particles read from a file
  // and passed round robin
  // Largest file will have a number of buffer chunks to send if it is too large
  // Every processor must send that number of chunks even if its own file
  // does not have that much information

  if (ENFORCE_MAX_READ == true && this->maxParticles > MAX_READ) {
    this->maxRead = MAX_READ;
    this->maxReadsPerFile = (this->maxParticles / this->maxRead) + 1;
  } else {
    this->maxRead = this->maxParticles;
    this->maxReadsPerFile = 1;
  }

  // Allocate space to hold buffer information for reading of files
  // Mass is constant use that float to store the tag
  // Number of particles is the first integer in the buffer
  int bufferSize = sizeof(int) + (this->maxRead * RECORD_SIZE);
  Message* message1 = new Message(bufferSize);
  Message* message2 = new Message(bufferSize);

  // Allocate space for the data read from the file
  POSVEL_T *fBlock = 0;
  POSVEL_T *lBlock = 0;
  POSVEL_T *vBlock = 0;
  ID_T* iBlock = 0;

  // RECORD format reads one particle at a time
  if (this->inputType == RECORD) {
    fBlock = new POSVEL_T[COSMO_FLOAT];
    iBlock = new ID_T[COSMO_INT];
  }

  // BLOCK format reads all particles at one time for triples
  else if (this->inputType == BLOCK) {
    lBlock = new POSVEL_T[this->maxRead * DIMENSION];
    vBlock = new POSVEL_T[this->maxRead * DIMENSION];
    iBlock = new ID_T[this->maxRead];
  }

  // Reserve particle storage to minimize reallocation
  int reserveSize = (int) (this->maxFiles * this->maxParticles * DEAD_FACTOR);

  // If multiple processors are reading the same file we can reduce size
  reserveSize /= this->processorsPerFile;

  if(reserveQ) {
#ifndef USE_VTK_COSMO
    cout << "readParticlesRoundRobin reserving vectors" << endl;
#endif
    this->xx->reserve(reserveSize);
    this->yy->reserve(reserveSize);
    this->zz->reserve(reserveSize);
    this->vx->reserve(reserveSize);
    this->vy->reserve(reserveSize);
    this->vz->reserve(reserveSize);
    this->ms->reserve(reserveSize);
    this->tag->reserve(reserveSize);
  }

  // Running total and index into particle data on this processor
  this->particleCount = 0;

  // Using the input files assigned to this processor, read the input
  // and push round robin to every other processor
  // this->maxFiles is the maximum number to read on any processor
  // Some processors may have no files to read but must still participate
  // in the round robin distribution

  for (int file = 0; file < this->maxFiles; file++) {

    // Open file to read the data if any for this processor
    ifstream* inStream = 0;
    int firstParticle = 0;
    int numberOfParticles = 0;
    int remainingParticles = 0;

    if ((int)this->inFiles.size() > file) {
      inStream = new ifstream(this->inFiles[file].c_str(), ios::in|ios::binary);

#ifndef USE_VTK_COSMO
      cout << "Rank " << this->myProc << " open file " << inFiles[file]
           << " with " << this->fileParticles[file] << " particles" << endl;
#endif

      // Number of particles read at one time depends on MPI buffer size
      numberOfParticles = this->fileParticles[file];
      if (numberOfParticles > this->maxRead)
        numberOfParticles = this->maxRead;

      // If a file is too large to be passed as an MPI message divide it up
      remainingParticles = this->fileParticles[file];

    } else {
#ifndef USE_VTK_COSMO
      cout << "Rank " << this->myProc << " no file to open " << endl;
#endif
    }

    for (int piece = 0; piece < this->maxReadsPerFile; piece++) {

      // Reset each MPI message for each file read
      message1->reset();
      message2->reset();

      // Processor has a file to read and share via round robin with others
      if (file < (int)this->inFiles.size()) {
        if (this->inputType == RECORD) {
          readFromRecordFile(inStream, firstParticle, numberOfParticles,
                             fBlock, iBlock, message1);
        } else {
          readFromBlockFile(inStream, firstParticle, numberOfParticles,
                           this->fileParticles[file],
                           lBlock, vBlock, iBlock, message1);
        }
        firstParticle += numberOfParticles;
        remainingParticles -= numberOfParticles;
        if (remainingParticles <= 0)
          numberOfParticles = 0;
        else if (remainingParticles < numberOfParticles)
          numberOfParticles = remainingParticles;
      }

      // Processor does not have a file to open but must participate in the
      // round robin with an empty buffer
      else {
        // Store number of particles used in first position
        int zero = 0;
        message1->putValue(&zero);
      }

      // Particles belonging to this processor are put in vectors
      distributeParticles(message1, message2);
    }

    // Can delete the read buffers as soon as last file is read because
    // information has been transferred into the double buffer1
    if (file == (this->maxFiles - 1)) {
      if (this->inputType == RECORD) {
        delete [] fBlock;
        delete [] iBlock;
      } else if (this->inputType == BLOCK) {
        delete [] lBlock;
        delete [] vBlock;
        delete [] iBlock;
      }
    }

    if ((int)this->inFiles.size() > file)
      inStream->close();
  }

  // After all particles have been distributed to vectors the double
  // buffers can be deleted
  delete message1;
  delete message2;

  // Count the particles across processors
  long totalAliveParticles = 0;
#ifdef USE_SERIAL_COSMO
  totalAliveParticles = this->numberOfAliveParticles;
#else
  MPI_Allreduce((void*) &this->numberOfAliveParticles,
                (void*) &totalAliveParticles,
                1, MPI_LONG, MPI_SUM, Partition::getComm());
#endif

#ifndef USE_VTK_COSMO
#ifdef DEBUG
  cout << "Rank " << setw(3) << this->myProc
       << " #alive = " << this->numberOfAliveParticles << endl;
#endif

  if (this->myProc == MASTER) {
    cout << "TotalAliveParticles " << totalAliveParticles << endl;
  }
#endif
  }
}

/////////////////////////////////////////////////////////////////////////
//
// Using the base name of the data, go to the subdirectory and determine
// how many input files there are.  Parcel those files between all the
// processors which will be responsible for actually reading 0 or more.
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::partitionInputFiles(bool force1PPF)
{
  // Find number of input files for this problem given the base input name
  // Get the subdirectory name containing the input files
  string::size_type dirPos = this->baseFile.rfind("/");
  string subdirectory;
  string baseName;

  // If the directory is not given use the current directory
  if (dirPos == string::npos) {
    subdirectory = "./";
    baseName = this->baseFile;
  } else {
    subdirectory = this->baseFile.substr(0, dirPos + 1);
    baseName = this->baseFile.substr(dirPos + 1);
  }

  // strip everything back to the first non-number
  string::size_type pos = baseName.size() - 1;
  int numbersOK = 1;

  while(numbersOK)
    {
    if(baseName[pos] >= '0' && baseName[pos] <= '9')
      {
      if(pos > 0)
        {
        pos = pos - 1;
        }
      else
        {
        break;
        }
      }
    else
      {
      numbersOK = 0;
      }
    }

  // base name is everything up to the numbers
  baseName = baseName.substr(0, pos + 1);

  // Open the subdirectory and make a list of input files
  DIR* directory = opendir(subdirectory.c_str());
  struct dirent* directoryEntry;
  vector<string> files;

  if (directory != NULL) {
  while ((directoryEntry = readdir(directory)))
    {
    // get the name
    string fileName = directoryEntry->d_name;
    pos = fileName.find(baseName.c_str());

    // if it starts with the base name
    if(pos == 0)
      {
      // check to see if it is all numbers on the end
      pos = baseName.size() + 1;
      numbersOK = 1;

      while(pos < fileName.size())
        {
        if(fileName[pos] < '0' || fileName[pos] > '9')
          {
          numbersOK = 0;
          break;
          }

        pos = pos + 1;
        }

      if(numbersOK)
        {
        fileName = subdirectory + fileName;
        files.push_back(fileName);
        }
      }
    }

  closedir(directory);
  }

  this->numberOfFiles = (int)files.size();

  if (this->numberOfFiles == 0) {
#ifdef USE_VTK_COSMO
    vtkStdString temp = "Processor ";
    temp += this->myProc;
    temp += " found no input files.\n";
    vtkOutputWindowDisplayErrorText(temp.c_str());

    return;
#else
    cout << "Rank " << this->myProc << " found no input files" << endl;
    exit(1);
#endif
  }

#ifndef USE_VTK_COSMO
#ifdef DEBUG
  if (this->myProc == MASTER) {
    for (int i = 0; i < this->numberOfFiles; i++)
      cout << "   File " << i << ": " << files[i] << endl;
  }
#endif
#endif

  // Divide the files between all the processors
  // If there are 1 or more files per processor set the
  // buffering up with a full round robin between all processors
  if (this->numberOfFiles >= this->numProc || force1PPF) {

    // Number of round robin sends to share all the files
    this->processorsPerFile = 1;
    this->numberOfFileSends = this->numProc - 1;
    this->maxFileSends = this->numberOfFileSends;

    // Which files does this processor read
    for (int i = 0; i < this->numberOfFiles; i++)
      if ((i % this->numProc) == this->myProc)
        this->inFiles.push_back(files[i]);

    // Where is the file sent, and where is it received
    if (this->myProc == this->numProc - 1)
      this->nextProc = 0;
    else
      this->nextProc = this->myProc + 1;
    if (this->myProc == 0)
      this->prevProc = this->numProc - 1;
    else
      this->prevProc = this->myProc - 1;
  }

  // If there are more processors than file set up as many round robin loops
  // as possible so that multiple processors read the same file. If the number
  // of files does not divide evenly into the number of processors the last
  // round robin loop will be bigger and some processors will contribute
  // buffers of 0 size to send

  else {

    // Assign the round robin circle (last circle is bigger than others)
    this->processorsPerFile = this->numProc / this->numberOfFiles;
    int numberOfRoundRobinCircles = this->processorsPerFile;
    int myCircle = this->myProc / this->numberOfFiles;
    int extraProcessors = this->numProc -
            (numberOfRoundRobinCircles * this->numberOfFiles);
    if (myCircle == numberOfRoundRobinCircles)
      myCircle--;

    int firstInCircle = myCircle * this->numberOfFiles;
    int lastInCircle = firstInCircle + this->numberOfFiles - 1;
    if (myCircle == (numberOfRoundRobinCircles - 1))
      lastInCircle += extraProcessors;

    // How big is the round robin circle this processor is in
    // What is the biggest round robin circle (needed because of MPI_Barrier)
    this->numberOfFileSends = lastInCircle - firstInCircle;
    this->maxFileSends = this->numberOfFiles + extraProcessors;

    // Which file does this processor read
    int index = this->myProc % this->numberOfFiles;
    if (myCircle == (this->myProc / this->numberOfFiles))
      this->inFiles.push_back(files[index]);

    // Where is the file sent, and where is it received
    if (this->myProc == lastInCircle)
      this->nextProc = firstInCircle;
    else
      this->nextProc = this->myProc + 1;
    if (this->myProc == firstInCircle)
      this->prevProc = lastInCircle;
    else
      this->prevProc = this->myProc - 1;
  }
}

/////////////////////////////////////////////////////////////////////////
//
// Open each input file belonging to this processor and find the number
// of particles for setting buffer sizes
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::findFileParticleCount()
{
  // Compute the total number of particles in the problem
  // Compute the maximum number of particles in any one file to set buffer size
  long numberOfParticles = 0;
  long maxNumberOfParticles = 0;
  int numberOfMyFiles = (int)this->inFiles.size();

  // Each processor counts the particles in its own files
  for (int i = 0; i < numberOfMyFiles; i++) {

    // Open my file
    ifstream *inStream = new ifstream(this->inFiles[i].c_str(), ios::in);
    if (inStream->fail()) {
      delete inStream;
#ifdef USE_VTK_COSMO
      vtkStdString message = "File ";
      message += this->inFiles[i];
      message += " cannot be opened.\n";
      vtkOutputWindowDisplayErrorText(message.c_str());

      this->totalParticles = 0;
      this->maxParticles = 0;
      return;
#else
      cout << "File: " << this->inFiles[i] << " cannot be opened" << endl;
      exit (-1);
#endif
    }

    if (this->inputType == RECORD) {

      // Compute the number of particles from file size
      inStream->seekg(0L, ios::end);
      int numberOfRecords = inStream->tellg() / RECORD_SIZE;
      this->fileParticles.push_back(numberOfRecords);

      numberOfParticles += numberOfRecords;
      if (maxNumberOfParticles < numberOfRecords)
        maxNumberOfParticles = numberOfRecords;
    }

    else if (this->inputType == BLOCK) {

      // Find the number of particles in the header
      readGadgetHeader(inStream);

      int numberOfRecords = this->gadgetParticleCount;
      this->fileParticles.push_back(numberOfRecords);

      numberOfParticles += numberOfRecords;
      if (maxNumberOfParticles < numberOfRecords)
        maxNumberOfParticles = numberOfRecords;
    }

    inStream->close();
    delete inStream;
  }

  // If multiple processors read the same file, just do the reduce on one set
  if (this->processorsPerFile > 1) {
    if (this->myProc >= this->numberOfFiles) {
      numberOfParticles = 0;
      maxNumberOfParticles = 0;
    }
  }

  // Share the information about total particles
#ifdef USE_SERIAL_COSMO
  this->totalParticles = numberOfParticles;
#else
  MPI_Allreduce((void*) &numberOfParticles,
                (void*) &this->totalParticles,
                1, MPI_LONG, MPI_SUM, Partition::getComm());
#endif

  // Share the information about max particles in a file for setting buffer size
#ifdef USE_SERIAL_COSMO
  this->maxParticles = maxNumberOfParticles;
#else
  MPI_Allreduce((void*) &maxNumberOfParticles,
                (void*) &this->maxParticles,
                1, MPI_LONG, MPI_MAX, Partition::getComm());
#endif

  // Share the maximum number of files on a processor for setting the loop
#ifdef USE_SERIAL_COSMO
  this->maxFiles = numberOfMyFiles;
#else
  MPI_Allreduce((void*) &numberOfMyFiles,
                (void*) &this->maxFiles,
                1, MPI_INT, MPI_MAX, Partition::getComm());
#endif

#ifndef USE_VTK_COSMO
#ifdef DEBUG
  if (this->myProc == MASTER) {
    cout << "Total particle count: " << this->totalParticles << endl;
    cout << "Max particle count:   " << this->maxParticles << endl;
  }
#endif
#endif
}

/////////////////////////////////////////////////////////////////////////
//
// Each processor reads 0 or more files, a buffer at a time.
// The particles are processed by seeing if they are in the subextent of
// this processor and are tagged either ALIVE or if dead, by the index of
// the neighbor zone which contains that particle.  That buffer is sent
// round robin to (myProc + 1) % numProc where it is processed and sent on.
// After each processor reads one buffer and sends and receives numProc - 1
// times the next buffer from the file is read.  Must use a double buffering
// scheme so that on each send/recv we switch buffers.
//
// Input files may be BLOCK or RECORD structured
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::distributeParticles(
                Message* message1,      // Send/receive buffers
                Message* message2)      // Send/receive buffers
{
  // Each processor has filled a buffer with particles read from a file
  // or had no particles to read but set the count in the buffer to 0
  // Process the buffer to keep only those within range
  Message* recvMessage = message1;
  Message* sendMessage = message2;

  // Process the original send buffer of particles from the file
  collectLocalParticles(recvMessage, sendMessage);

  // Distribute buffer round robin so that all processors see it
  for (int step = 0; step < this->maxFileSends; step++) {

    if (step < this->numberOfFileSends)
      {
      // Send buffer to the next processor if round robin loop is still active
      sendMessage->send(this->nextProc);

      // Receive buffer from the previous processor
      recvMessage->receive(this->prevProc);
      }

#ifndef USE_SERIAL_COSMO
    MPI_Barrier(Partition::getComm());
#endif

    // Process the send buffer for alive and dead before sending on
    // the particles that were not claimed by this processor
    if (step < this->numberOfFileSends)
      collectLocalParticles(recvMessage, sendMessage);

#ifndef USE_SERIAL_COSMO
    MPI_Barrier(Partition::getComm());
#endif
  }
}

/////////////////////////////////////////////////////////////////////////////
//
// Input file is RECORD structured so read each particle record and populate
// the double buffer in particle order for the rest of the processing
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromRecordFile(
                        ifstream* inStream,     // Stream to read from
                        int firstParticle,      // First particle index
                        int numberOfParticles,  // Number to read this time
                        POSVEL_T* fBlock,       // Buffer for read in data
                        ID_T* iBlock,           // Buffer for read in data
                        Message* message)       // Reordered data
{
  // Store number of particles used in first position
  message->putValue(&numberOfParticles);
  if (numberOfParticles == 0)
    return;

  // Seek to the first particle locations and read
  int skip = RECORD_SIZE * firstParticle;
  inStream->seekg(skip, ios::beg);

  // Store each particle location, velocity, mass and tag (as float) in buffer
  int changeCount = 0;
  for (int p = 0; p < numberOfParticles; p++) {

    // Set file pointer to the requested particle
    inStream->read(reinterpret_cast<char*>(fBlock),
                   COSMO_FLOAT * sizeof(POSVEL_T));

    if (inStream->gcount() != COSMO_FLOAT * sizeof(POSVEL_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // Convert units if requested
    fBlock[0] *= this->distConvertFactor;
    fBlock[2] *= this->distConvertFactor;
    fBlock[4] *= this->distConvertFactor;
    fBlock[6] *= this->massConvertFactor;

    inStream->read(reinterpret_cast<char*>(iBlock),
                   COSMO_INT * sizeof(ID_T));

    if (inStream->gcount() != COSMO_INT * sizeof(ID_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // If the location is not within the bounding box wrap around
    for (int i = 0; i <= 4; i = i + 2) {
      if (fBlock[i] >= this->boxSize) {
#ifndef USE_VTK_COSMO
#ifdef DEBUG
        cout << "Location at " << i << " changed from " << fBlock[i] << endl;
#endif
#endif
        fBlock[i] -= this->boxSize;
        changeCount++;
      }
    }

    // Store location and velocity and mass in message buffer
    // Reorder so that location vector is followed by velocity vector
    message->putValue(&fBlock[0]);
    message->putValue(&fBlock[2]);
    message->putValue(&fBlock[4]);
    message->putValue(&fBlock[1]);
    message->putValue(&fBlock[3]);
    message->putValue(&fBlock[5]);
    message->putValue(&fBlock[6]);

    // Store the integer tag
    message->putValue(&iBlock[0]);
  }
}

/////////////////////////////////////////////////////////////////////////////
//
// Input file is BLOCK structured so read head and each block of data.
// Gadget format:
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_H 4 bytes (size of header)
//    Header (6 types of particles with counts and masses)
//    SKIP_H 4 bytes (size of header)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_L 4 bytes (size of location block in bytes)
//    Block of location data where each particle's x,y,z is stored together
//    SKIP_L 4 bytes (size of location block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_V 4 bytes (size of velocity block in bytes)
//    Block of velocity data where each particle's xv,yv,zv is stored together
//    SKIP_V 4 bytes (size of velocity block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_T 4 bytes (size of tag block in bytes)
//    Block of tag data
//    SKIP_T 4 bytes (size of tag block in bytes)
//
// Reorder the data after it is read into the same structure as the
// RECORD data so that the rest of the code does not have to be changed
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromBlockFile(
                        ifstream* inStream,     // Stream to read from
                        int firstParticle,      // First particle index
                        int numberOfParticles,  // Number to read this time
                        int totParticles,       // Total particles in file
                        POSVEL_T* lBlock,       // Buffer for read of location
                        POSVEL_T* vBlock,       // Buffer for read of velocity
                        ID_T* iBlock,           // Buffer for read in data
                        Message* message)       // Reordered data
{
  // Store number of particles used in first position
  message->putValue(&numberOfParticles);
  if (numberOfParticles == 0)
    return;

  // Calculate skips to first location, velocity and tag
  int skipToLocation = 0;
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;		// Size of header
  skipToLocation += GADGET_HEADER_SIZE;		// Header
  skipToLocation += GADGET_SKIP;		// Size of header
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;		// Size of location block

  int skipToVelocity = skipToLocation;
  skipToVelocity += DIMENSION * sizeof(POSVEL_T) * totParticles;
  skipToVelocity += GADGET_SKIP;		// Size of location block
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToVelocity += GADGET_SKIP;		// Size of velocity block

  int skipToTag = skipToVelocity;
  skipToTag += DIMENSION * sizeof(POSVEL_T) * totParticles;
  skipToTag += GADGET_SKIP;			// Size of velocity block
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToTag += GADGET_SKIP;			// Size of tag block

  // Seek to the first requested particle location and read triples
  inStream->seekg(skipToLocation, ios::beg);
  int skip = (DIMENSION * sizeof(POSVEL_T) * firstParticle);
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) lBlock, sizeof(POSVEL_T),
                 DIMENSION * numberOfParticles, inStream);

  // Convert units of distance
  for (int i = 0; i < DIMENSION*numberOfParticles; i++)
    lBlock[i] *= this->distConvertFactor;

  // If the location is not within the bounding box wrap around
  for (int i = 0; i < DIMENSION*numberOfParticles; i++) {
    if (lBlock[i] >= this->boxSize)
      lBlock[i] -= this->boxSize;
  }

  // Seek to first requested particle velocity and read triples
  inStream->seekg(skipToVelocity, ios::beg);
  skip = (DIMENSION * sizeof(POSVEL_T) * firstParticle); // skip to velocity
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) vBlock, sizeof(POSVEL_T),
                 DIMENSION * numberOfParticles, inStream);

  // Seek to first requested particle tag and read
  inStream->seekg(skipToTag, ios::beg);
  skip = sizeof(ID_T) * firstParticle;             // skip to tag
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) iBlock, sizeof(ID_T),
                 numberOfParticles, inStream);

  // Store the locations in the message buffer in record order
  // so that the same distribution method for RECORD will work
  int particlesRemaining = numberOfParticles;
  int curParticle = firstParticle;
  int type = 0;
  int indx = 0;
  int tagindx = 0;

  // When more than one gadget particle type is in the file we must
  // know what mass to assign to the current piece read
  while (particlesRemaining > 0) {

    // Set particle type and mass based on current particle
    while (type < NUM_GADGET_TYPES && curParticle >= this->gadgetStart[type])
      type++;
    type--;

    POSVEL_T particleMass =
      (POSVEL_T) this->gadgetHeader.mass[type] * this->massConvertFactor;

    // Place particles of this type and mass in the buffer
    int numLeftInType = this->gadgetHeader.npart[type] - 
                        (curParticle - this->gadgetStart[type]);
    int count = min(particlesRemaining, numLeftInType);

    for (int p = 0; p < count; p++) {

      // Locations
      message->putValue(&lBlock[indx]);           // X location
      message->putValue(&lBlock[indx+1]);         // Y location
      message->putValue(&lBlock[indx+2]);         // Z location

      // Velocities
      message->putValue(&vBlock[indx]);           // X velocity
      message->putValue(&vBlock[indx+1]);         // Y velocity
      message->putValue(&vBlock[indx+2]);         // Z velocity

      // Mass
      message->putValue(&particleMass);

      // Id tag
      message->putValue(&iBlock[tagindx]);
      indx += DIMENSION;
      tagindx++;
    }

    // Do we have more particles of a different type
    particlesRemaining -= count;
    curParticle += count;
  }
}

#ifndef USE_SERIAL_COSMO
/////////////////////////////////////////////////////////////////////////////
//
// Input file is RECORD structured so read each particle record and populate
// the double buffer in particle order for the rest of the processing
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromRecordFile(
                        ifstream* inStream,     // Stream to read from
                        int firstParticle,      // First particle index
                        int numberOfParticles,  // Number to read this time
                        POSVEL_T* fBlock,       // Buffer for read in data
                        ID_T* iBlock,           // Buffer for read in data
                        std::vector< std::vector<CosmoParticle> > &pByProc)
{
  if (numberOfParticles == 0)
    return;

  // Seek to the first particle locations and read
  int skip = RECORD_SIZE * firstParticle;
  inStream->seekg(skip, ios::beg);

  // Store each particle location, velocity, mass and tag (as float) in buffer
  int changeCount = 0;
  for (int p = 0; p < numberOfParticles; p++) {

    // Set file pointer to the requested particle
    inStream->read(reinterpret_cast<char*>(fBlock),
                   COSMO_FLOAT * sizeof(POSVEL_T));

    if (inStream->gcount() != COSMO_FLOAT * sizeof(POSVEL_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // Convert units if requested
    fBlock[0] *= this->distConvertFactor;
    fBlock[2] *= this->distConvertFactor;
    fBlock[4] *= this->distConvertFactor;
    fBlock[6] *= this->massConvertFactor;

    inStream->read(reinterpret_cast<char*>(iBlock),
                   COSMO_INT * sizeof(ID_T));

    if (inStream->gcount() != COSMO_INT * sizeof(ID_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // If the location is not within the bounding box wrap around
    for (int i = 0; i <= 4; i = i + 2) {
      if (fBlock[i] >= this->boxSize) {
#ifndef USE_VTK_COSMO
#ifdef DEBUG
        cout << "Location at " << i << " changed from " << fBlock[i] << endl;
#endif
#endif
        fBlock[i] -= this->boxSize;
        changeCount++;
      }
    }

    // Figure out to which rank this particle belongs
    float sizeX = this->boxSize / this->layoutSize[0];
    float sizeY = this->boxSize / this->layoutSize[1];
    float sizeZ = this->boxSize / this->layoutSize[2];
    int coords[3] = { (int) (fBlock[0]/sizeX),
                      (int) (fBlock[2]/sizeY),
                      (int) (fBlock[4]/sizeZ)
                    };
    int home;
    MPI_Cart_rank(Partition::getComm(), coords, &home);

    // Store location and velocity and mass
    // Reorder so that location vector is followed by velocity vector
    CosmoParticle particle = { { fBlock[0],
                                 fBlock[2],
                                 fBlock[4],
                                 fBlock[1],
                                 fBlock[3],
                                 fBlock[5],
                                 fBlock[6]
                               }, {
                                 static_cast<int>(iBlock[0])
                               }
                             };
    pByProc[home].push_back(particle);
  }
}

/////////////////////////////////////////////////////////////////////////////
//
// Input file is BLOCK structured so read head and each block of data.
// Gadget format:
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_H 4 bytes (size of header)
//    Header (6 types of particles with counts and masses)
//    SKIP_H 4 bytes (size of header)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_L 4 bytes (size of location block in bytes)
//    Block of location data where each particle's x,y,z is stored together
//    SKIP_L 4 bytes (size of location block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_V 4 bytes (size of velocity block in bytes)
//    Block of velocity data where each particle's xv,yv,zv is stored together
//    SKIP_V 4 bytes (size of velocity block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_T 4 bytes (size of tag block in bytes)
//    Block of tag data
//    SKIP_T 4 bytes (size of tag block in bytes)
//
// Reorder the data after it is read into the same structure as the
// RECORD data so that the rest of the code does not have to be changed
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromBlockFile(
                        ifstream* inStream,     // Stream to read from
                        int firstParticle,      // First particle index
                        int numberOfParticles,  // Number to read this time
                        int totParticles,       // Total particles in file
                        POSVEL_T* lBlock,       // Buffer for read of location
                        POSVEL_T* vBlock,       // Buffer for read of velocity
                        ID_T* iBlock,           // Buffer for read in data
                        std::vector< std::vector<CosmoParticle> > &pByProc)
{
  if (numberOfParticles == 0)
    return;

  // Calculate skips to first location, velocity and tag
  int skipToLocation = 0;
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;		// Size of header
  skipToLocation += GADGET_HEADER_SIZE;		// Header
  skipToLocation += GADGET_SKIP;		// Size of header
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;		// Size of location block

  int skipToVelocity = skipToLocation;
  skipToVelocity += DIMENSION * sizeof(POSVEL_T) * totParticles;
  skipToVelocity += GADGET_SKIP;		// Size of location block
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToVelocity += GADGET_SKIP;		// Size of velocity block

  int skipToTag = skipToVelocity;
  skipToTag += DIMENSION * sizeof(POSVEL_T) * totParticles;
  skipToTag += GADGET_SKIP;			// Size of velocity block
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToTag += GADGET_SKIP;			// Size of tag block

  // Seek to the first requested particle location and read triples
  inStream->seekg(skipToLocation, ios::beg);
  int skip = (DIMENSION * sizeof(POSVEL_T) * firstParticle);
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) lBlock, sizeof(POSVEL_T),
                 DIMENSION * numberOfParticles, inStream);

  // Convert units of distance
  for (int i = 0; i < DIMENSION*numberOfParticles; i++)
    lBlock[i] *= this->distConvertFactor;

  // If the location is not within the bounding box wrap around
  for (int i = 0; i < DIMENSION*numberOfParticles; i++) {
    if (lBlock[i] >= this->boxSize)
      lBlock[i] -= this->boxSize;
  }

  // Seek to first requested particle velocity and read triples
  inStream->seekg(skipToVelocity, ios::beg);
  skip = (DIMENSION * sizeof(POSVEL_T) * firstParticle); // skip to velocity
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) vBlock, sizeof(POSVEL_T),
                 DIMENSION * numberOfParticles, inStream);

  // Seek to first requested particle tag and read
  inStream->seekg(skipToTag, ios::beg);
  skip = sizeof(ID_T) * firstParticle;             // skip to tag
  inStream->seekg(skip, ios::cur);

  readData(this->gadgetSwap, (void*) iBlock, sizeof(ID_T),
                 numberOfParticles, inStream);

  // Store the locations in the message buffer in record order
  // so that the same distribution method for RECORD will work
  int particlesRemaining = numberOfParticles;
  int curParticle = firstParticle;
  int type = 0;
  int indx = 0;
  int tagindx = 0;

  // When more than one gadget particle type is in the file we must
  // know what mass to assign to the current piece read
  while (particlesRemaining > 0) {

    // Set particle type and mass based on current particle
    while (type < NUM_GADGET_TYPES && curParticle >= this->gadgetStart[type])
      type++;
    type--;

    POSVEL_T particleMass =
      (POSVEL_T) this->gadgetHeader.mass[type] * this->massConvertFactor;

    // Place particles of this type and mass in the buffer
    int numLeftInType = this->gadgetHeader.npart[type] - 
                        (curParticle - this->gadgetStart[type]);
    int count = min(particlesRemaining, numLeftInType);

    float sizeX = this->boxSize / this->layoutSize[0];
    float sizeY = this->boxSize / this->layoutSize[1];
    float sizeZ = this->boxSize / this->layoutSize[2];

    for (int p = 0; p < count; p++) {
      // Figure out to which rank this particle belongs
      int coords[3] = { (int) (lBlock[indx]/sizeX),
                        (int) (lBlock[indx+1]/sizeY),
                        (int) (lBlock[indx+2]/sizeZ)
                      };
      int home;
      MPI_Cart_rank(Partition::getComm(), coords, &home);

      // Store location and velocity and mass
      // Reorder so that location vector is followed by velocity vector
      CosmoParticle particle = { { lBlock[indx],   // X location
                                   lBlock[indx+1], // Y location
                                   lBlock[indx+2], // Z location
                                   vBlock[indx],   // X velocity
                                   vBlock[indx+1], // Y velocity
                                   vBlock[indx+2], // Z velocity
                                   particleMass
                                 }, {
                                   static_cast<int>(iBlock[tagindx])
                                 }
                               };
      pByProc[home].push_back(particle);

      indx += DIMENSION;
      tagindx++;
    }

    // Do we have more particles of a different type
    particlesRemaining -= count;
    curParticle += count;
  }
}
#endif // USE_SERIAL_COSMO

/////////////////////////////////////////////////////////////////////////////
//
// Process the data buffer of particles to choose those which are ALIVE
// or DEAD on this processor.  Do wraparound tests to populate as for a
// 3D torus.  Dead particle status is the zone id of the neighbor processor
// which contains it as an ALIVE particle.
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::collectLocalParticles(
                Message* recvMessage,      // Read particles and extract
                Message* sendMessage)      // Other particles copied here
{
  // In order to read a buffer, reset position to the beginning
  recvMessage->reset();
  sendMessage->reset();

  int recvParticles;
  int sendParticles = 0;
  recvMessage->getValue(&recvParticles);
  sendMessage->putValue(&sendParticles);

  POSVEL_T loc[DIMENSION], vel[DIMENSION], mass;
  ID_T id;

  // Test each particle in the buffer to see if it is ALIVE or DEAD
  // If it is DEAD assign it to the neighbor zone that it is in
  // Check all combinations of wraparound

  for (int i = 0; i < recvParticles; i++) {
    for (int dim = 0; dim < DIMENSION; dim++)
      recvMessage->getValue(&loc[dim]);
    for (int dim = 0; dim < DIMENSION; dim++)
      recvMessage->getValue(&vel[dim]);
    recvMessage->getValue(&mass);
    recvMessage->getValue(&id);

    // Is the particle ALIVE on this processor
    if ((loc[0] >= minAlive[0] && loc[0] < maxAlive[0]) &&
        (loc[1] >= minAlive[1] && loc[1] < maxAlive[1]) &&
        (loc[2] >= minAlive[2] && loc[2] < maxAlive[2])) {

          this->xx->push_back(loc[0]);
          this->yy->push_back(loc[1]);
          this->zz->push_back(loc[2]);
          this->vx->push_back(vel[0]);
          this->vy->push_back(vel[1]);
          this->vz->push_back(vel[2]);
          this->ms->push_back(mass);
          this->tag->push_back(id);

          this->numberOfAliveParticles++;
          this->particleCount++;
    } else {

      // Pass the particle along to the next processor in send buffer
      sendParticles++;
      for (int dim = 0; dim < DIMENSION; dim++)
        sendMessage->putValue(&loc[dim]);
      for (int dim = 0; dim < DIMENSION; dim++)
        sendMessage->putValue(&vel[dim]);
      sendMessage->putValue(&mass);
      sendMessage->putValue(&id);
    }
  }
  // Overwrite the send buffer first word with the known number of particles
  sendMessage->putValueAtPosition(&sendParticles, 0);
}

/////////////////////////////////////////////////////////////////////////
//
// Each processor reads 1 file or gets a pointer to data eventually
// As the particle is read it will be stored as an alive particle on this
// processor and will be checked about neighbor ranges to see if it must
// be exchanged
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readParticlesOneToOne(int reserveQ)
{
  // File name is the base file name with processor id appended
  // Because MPI Cartesian topology is used the arrangement of files in
  // physical space must follow the rule of last dimension varies fastest
  ostringstream fileName;
  fileName << this->baseFile << this->myProc;
  this->inFiles.push_back(fileName.str());

  // Compute the total number of particles in the problem
  // Compute the maximum number of particles in any one file to set buffer size
  findFileParticleCount();

  // Reserve particle storage to minimize reallocation
  int reserveSize = (int) (this->maxParticles * DEAD_FACTOR);

  if(reserveQ) {
#ifndef USE_VTK_COSMO
    cout << "readParticlesOneToOne reserving vectors" << endl;
#endif
    this->xx->reserve(reserveSize);
    this->yy->reserve(reserveSize);
    this->zz->reserve(reserveSize);
    this->vx->reserve(reserveSize);
    this->vy->reserve(reserveSize);
    this->vz->reserve(reserveSize);
    this->ms->reserve(reserveSize);
    this->tag->reserve(reserveSize);
  }

  // Running total and index into particle data on this processor
  this->particleCount = 0;

  // Read the input file storing particles immediately because all are alive
  if (this->inputType == RECORD) {
    readFromRecordFile();
  } else {
    readFromBlockFile();
  }
}

/////////////////////////////////////////////////////////////////////////////
//
// Input file is RECORD structured so read each particle record and populate
// the vectors of particles marking all as ALIVE
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromRecordFile()
{
  // Only one file per processor named in index 0
  ifstream inStream(this->inFiles[0].c_str(), ios::in);
  int numberOfParticles = this->fileParticles[0];

#ifndef USE_VTK_COSMO
  cout << "Rank " << this->myProc << " open file " << this->inFiles[0]
       << " with " << numberOfParticles << " particles" << endl;
#endif

  POSVEL_T* fBlock = new POSVEL_T[COSMO_FLOAT];
  ID_T* iBlock = new ID_T[COSMO_INT];

  // Store each particle location, velocity and tag
  for (int i = 0; i < numberOfParticles; i++) {

    // Set file pointer to the requested particle
    inStream.read(reinterpret_cast<char*>(fBlock),
                   COSMO_FLOAT * sizeof(POSVEL_T));

    if (inStream.gcount() != COSMO_FLOAT * sizeof(POSVEL_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      inStream.close();
      delete [] fBlock;
      delete [] iBlock;

      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // Convert units if requested
    fBlock[0] *= this->distConvertFactor;
    fBlock[2] *= this->distConvertFactor;
    fBlock[4] *= this->distConvertFactor;
    fBlock[6] *= this->massConvertFactor;

    inStream.read(reinterpret_cast<char*>(iBlock),
                   COSMO_INT * sizeof(ID_T));

    if (inStream.gcount() != COSMO_INT * sizeof(ID_T)) {
#ifdef USE_VTK_COSMO
      vtkOutputWindowDisplayErrorText("Premature end-of-file.\n");
      inStream.close();
      delete [] fBlock;
      delete [] iBlock;

      return;
#else
      cout << "Premature end-of-file" << endl;
      exit (-1);
#endif
    }

    // Store information in buffer if within range on this processor
    if ((fBlock[0] >= minAlive[0] && fBlock[0] <= maxAlive[0]) &&
        (fBlock[2] >= minAlive[1] && fBlock[2] <= maxAlive[1]) &&
        (fBlock[4] >= minAlive[2] && fBlock[4] <= maxAlive[2])) {

      this->xx->push_back(fBlock[0]);
      this->vx->push_back(fBlock[1]);
      this->yy->push_back(fBlock[2]);
      this->vy->push_back(fBlock[3]);
      this->zz->push_back(fBlock[4]);
      this->vz->push_back(fBlock[5]);
      this->ms->push_back(fBlock[6]);
      this->tag->push_back(iBlock[0]);

      this->numberOfAliveParticles++;
      this->particleCount++;
    }
  }

  inStream.close();
  delete [] fBlock;
  delete [] iBlock;
}

/////////////////////////////////////////////////////////////////////////////
//
// Input file is BLOCK structured so read head and each block of data.
// Gadget format:
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_H 4 bytes (size of header)
//    Header (6 types of particles with counts and masses)
//    SKIP_H 4 bytes (size of header)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_L 4 bytes (size of location block in bytes)
//    Block of location data where each particle's x,y,z is stored together
//    SKIP_L 4 bytes (size of location block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_V 4 bytes (size of velocity block in bytes)
//    Block of velocity data where each particle's xv,yv,zv is stored together
//    SKIP_V 4 bytes (size of velocity block in bytes)
//
//    SKIP_GADGET_2 has extra 16 bytes
//    SKIP_T 4 bytes (size of tag block in bytes)
//    Block of tag data
//    SKIP_T 4 bytes (size of tag block in bytes)
//
// Reorder the data after it is read into the same structure as the
// RECORD data so that the rest of the code does not have to be changed
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readFromBlockFile()
{
  // Only one file per processor named in index 0
  ifstream inStream(this->inFiles[0].c_str(), ios::in);
  int numberOfParticles = this->fileParticles[0];

#ifndef USE_VTK_COSMO
  cout << "Rank " << this->myProc << " open file " << this->inFiles[0]
       << " with " << numberOfParticles << " particles" << endl;
#endif

  // Calculate skips to first location, velocity and tag
  int skipToLocation = 0;
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;                // Size of header
  skipToLocation += GADGET_HEADER_SIZE;         // Header
  skipToLocation += GADGET_SKIP;                // Size of header
  if (this->gadgetFormat == GADGET_2)
    skipToLocation += GADGET_2_SKIP;
  skipToLocation += GADGET_SKIP;                // Size of location block

  // Allocate blocks to read into
  POSVEL_T* lBlock = new POSVEL_T[numberOfParticles * DIMENSION];
  POSVEL_T* vBlock = new POSVEL_T[numberOfParticles * DIMENSION];
  ID_T* iBlock = new ID_T[numberOfParticles];

  // Seek to particle locations and read triples
  inStream.seekg(skipToLocation, ios::beg);
  readData(this->gadgetSwap, (void*) lBlock, sizeof(POSVEL_T),
           DIMENSION * numberOfParticles, &inStream);

  // Convert locations
  for (int p = 0; p < DIMENSION * numberOfParticles; p++)
    lBlock[p] *= this->distConvertFactor;

  // Seek to particle velocities and read triples
  inStream.seekg((2 * GADGET_SKIP), ios::cur);
  readData(this->gadgetSwap, (void*) vBlock, sizeof(POSVEL_T),
           DIMENSION * numberOfParticles, &inStream);


  // Seek to particle tags and read
  inStream.seekg((2 * GADGET_SKIP), ios::cur);
  readData(this->gadgetSwap, (void*) iBlock, sizeof(ID_T),
           numberOfParticles, &inStream);

  // Store mass, locations, velocities and tags into arrays if in range
  // Range test is needed because this code is used for ONE_TO_ONE where all
  // particles must be added, and by one single input file over many
  // processors where messaging is not needed, but some particles don't belong
  int indx = 0;
  int tagindx = 0;
  for (int type = 0; type < NUM_GADGET_TYPES; type++) {

    POSVEL_T particleMass =
      (POSVEL_T) this->gadgetHeader.mass[type] * this->massConvertFactor;

    for (int p = 0; p < this->gadgetHeader.npart[type]; p++) {

      if ((lBlock[indx] >= minAlive[0] && lBlock[indx] < maxAlive[0]) &&
          (lBlock[indx+1] >= minAlive[1] && lBlock[indx+1] < maxAlive[1]) &&
          (lBlock[indx+2] >= minAlive[2] && lBlock[indx+2] < maxAlive[2])) {

        this->xx->push_back(lBlock[indx]);
        this->yy->push_back(lBlock[indx+1]);
        this->zz->push_back(lBlock[indx+2]);
        this->vx->push_back(vBlock[indx]);
        this->vy->push_back(vBlock[indx+1]);
        this->vz->push_back(vBlock[indx+2]);
        this->ms->push_back(particleMass);
        this->tag->push_back(iBlock[tagindx]);

        this->numberOfAliveParticles++;
        this->particleCount++;
      }
      indx += DIMENSION;
      tagindx++;
    }
  }

  delete [] lBlock;
  delete [] vBlock;
  delete [] iBlock;
  inStream.close();
}

/////////////////////////////////////////////////////////////////////////////
//
// Read the Gadget header from the stream
// Gadget file may be Gadget-1 format with no block indicators or
// Gadget-2 format with size of block 4 byte integers surrounding each block
// Data may be big or little endian which we can tell by checking that
// the header size is 256 in the first 4 bytes
//
/////////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readGadgetHeader(ifstream* gStr)
{
  this->gadgetSwap = false;
  this->gadgetFormat = 1;
  int blockSize, blockSize2;
  string gadget2;

  // Set the gadget format type by reading the first 4 byte integer
  // If it is not "256" or "65536" then gadget-2 format with 16 bytes in front
  readData(this->gadgetSwap, (void*) &blockSize, GADGET_SKIP, 1, gStr);
  if (blockSize != GADGET_HEADER_SIZE && blockSize != GADGET_HEADER_SIZE_SWP) {
    this->gadgetFormat = GADGET_2;
    gadget2 = readString(gStr, GADGET_2_SKIP - GADGET_SKIP);
    readData(this->gadgetSwap, (void*) &blockSize, GADGET_SKIP, 1, gStr);
  }

  // Set the swap type
  if (blockSize != GADGET_HEADER_SIZE) {
    this->gadgetSwap = true;
    blockSize = GADGET_HEADER_SIZE;
  }

  // Read the Gadget header
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.npart[0],
                         sizeof(int), NUM_GADGET_TYPES, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.mass[0],
                         sizeof(double), NUM_GADGET_TYPES, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.time,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.redshift,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_sfr,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_feedback,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.npartTotal[0],
                         sizeof(int), NUM_GADGET_TYPES, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_cooling,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.num_files,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.BoxSize,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.Omega0,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.OmegaLambda,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.HubbleParam,
                         sizeof(double), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_stellarage,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_metals,
                         sizeof(int), 1, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.HighWord[0],
                         sizeof(int), NUM_GADGET_TYPES, gStr);
  readData(this->gadgetSwap, (void*) &this->gadgetHeader.flag_entropy,
                         sizeof(int), 1, gStr);
  string fill = readString(gStr, GADGET_FILL);
  strcpy(&this->gadgetHeader.fill[0], fill.c_str());

  // Read the Gadget header size to verify block
  readData(this->gadgetSwap, (void*) &blockSize2, GADGET_SKIP, 1, gStr);
  if (blockSize != blockSize2)
#ifdef USE_VTK_COSMO
    vtkOutputWindowDisplayErrorText("Mismatch of header size and header structure.\n");
#else
    cout << "Mismatch of header size and header structure" << endl;
#endif

  // Every type particle will have location, velocity and tag so sum up
  this->gadgetParticleCount = 0;
  for (int i = 0; i < NUM_GADGET_TYPES; i++) {
    this->gadgetStart[i] = this->gadgetParticleCount;
    this->gadgetParticleCount += this->gadgetHeader.npart[i];
  }
}

/////////////////////////////////////////////////////////////////////////
//
// Read in the requested number of characters
//
/////////////////////////////////////////////////////////////////////////

string ParticleDistribute::readString(ifstream* inStr, int size)
{
   char* buffer = new char[size + 1];
   inStr->read(buffer, size);
   buffer[size] = '\0';

   // Make sure string has legal values
   if (isalnum(buffer[0]) == 0)
      buffer[0] = '\0';
   for (int i = 1; i < size; i++)
      if (isprint(buffer[i]) == 0)
         buffer[i] = '\0';

   string retString = buffer;
   delete [] buffer;
   return retString;
}

/////////////////////////////////////////////////////////////////////////
//
// Read in the number of items from the file pointer and
// byte swap if necessary
//
/////////////////////////////////////////////////////////////////////////

void ParticleDistribute::readData(
        bool swap,
        void* data,
        unsigned long dataSize,
        unsigned long dataCount,
        ifstream* inStr)
{
   // Read all the data from the file
   inStr->read(reinterpret_cast<char*>(data), dataSize*dataCount);

   if (swap == true) {

      // Byte swap each integer
      char* dataPtr = (char*) data;
      char temp;
      for (unsigned long item = 0; item < dataCount; item++) {

         // Do a byte-by-byte swap, reversing the order.
         for (unsigned int i = 0; i < dataSize / 2; i++) {
            temp = dataPtr[i];
            dataPtr[i] = dataPtr[dataSize - 1 - i];
            dataPtr[dataSize - 1 - i] = temp;
         }
         dataPtr += dataSize;
      }
   }
}