//------------------------------------------------------------------------------
//
// Project:	AMR Visualization
// Module:	$RCSfile$
// Language:	C++
// Date:	$Date$
// Author:	$Author$
// Version:	$Revision$
//
//------------------------------------------------------------------------------

#include "AMRGrid.hh"

// Standard C/C++ Includes
#include <list>
#include <algorithm>

// Own basic includes
#include <assume.hh>
#include <interval.h>

// AMR related includes
#include <AmrGridReader.hh>
#include <AMRLevel.hh>

namespace AMRGridAuxFct {
  inline char *dataTypeToStr(IObase::DataType dt) {
    switch (dt) {
      case 0:
	return "Byte/Int8"; break;
      case 1:
	return "Int16"; break;
      case 2:
	return "Int32"; break;
      case 3:
	return "Int64"; break;
      case 4:
	return "Float32"; break;
      case 5:
	return "Float64"; break;
      case 6:
	return "uInt8/uChar"; break;
      case 7:
	return "uInt16"; break;
      case 8:
	return "uInt32"; break;
      case 9:
	return "uInt64"; break;
      case 10:
	return "Char/Char8/String"; break;
      case 11: 
	return "Unicode/Char16"; break;
      default:
	return "UNKNOWN/ERROR"; break;
    }
  }

  inline int ceilDiv(int divident, int divisor)
  {
    return (divident / divisor) + ((divident % divisor) == 0 ? 0 : 1);
  }

  std::ostream& operator<<(std::ostream &os, const int v[3])
  {
    os << "(" << v[0] << ", " << v[1] << ", " << v[2] << ")";
    return os;
  }
  std::ostream& operator<<(std::ostream &os, const double v[3])
  {
    os << "(" << v[0] << ", " << v[1] << ", " << v[2] << ")";
    return os;
  }
}

AMRGrid::GridProperty::~GridProperty()
{
}

AMRGrid::IntersectionMap::IntersectionMap(int iRes, int jRes, int kRes, int subCellIRes, int subCellJRes, int subCellKRes, MapType mapType) : mMapType(mapType), mIRes(iRes), mJRes(jRes), mKRes(kRes), mIntersectingGridIndices(new int[mIRes*mJRes*mKRes]), mSubCellIRes(subCellIRes), mSubCellJRes(subCellJRes), mSubCellKRes(subCellKRes), mRefinementMaps(new IntersectionMap*[mIRes*mJRes*mKRes])
{
  init();
}

AMRGrid::IntersectionMap::~IntersectionMap()
{
  delete[] mIntersectingGridIndices;
  for (int i=0; i<mIRes*mJRes*mKRes; ++i) 
    delete mRefinementMaps[i];
  delete[] mRefinementMaps;
}

void AMRGrid::IntersectionMap::init()
{
  for(int i=0; i<mIRes*mJRes*mKRes; ++i) {
    mIntersectingGridIndices[i] = NoGrid;
    mRefinementMaps[i] = NULL; 
  }
}

AMRGrid::IntersectionMap& AMRGrid::IntersectionMap::refineCell(int i, int j, int k)
{
  assert(0<=i && i<mIRes && 0<=j && j<mJRes && 0<=k && k<mKRes);
  int idx = idxForTriple(i, j, k);
  if (mRefinementMaps[idx]) return *mRefinementMaps[idx];
  else return *(mRefinementMaps[idx] = new IntersectionMap(mSubCellIRes, mSubCellJRes, mSubCellKRes, 1, 1, 1, RefinementMap));
}

void AMRGrid::IntersectionMap::createBlockRefined(int i, int j, int k, int subIMin, int subIMax, int subJMin, int subJMax, int subKMin, int subKMax, AMRHierarchy::GridId gridIdx)
{
  IntersectionMap &refinedCell = refineCell(i, j, k);
  for (int currI=subIMin; currI<=subIMax; ++currI) {
    for (int currJ=subJMin; currJ<=subJMax; ++currJ) {
      for (int currK=subKMin; currK<=subKMax; ++currK) {
	refinedCell(currI, currJ, currK) = gridIdx;
      }
    }
  }      
}

AMRGrid::AMRGrid(AMRHierarchy::GridId index, AMRHierarchy *aMRHierarchy, AmrGridReader *amrGridReader) : mGridIndex(index), mAMRHierarchy(aMRHierarchy), mRefiningGrids(), mRefinedGrids(), mNextGridInLevel(0), mPrevGridInLevel(0), mAmrGridReader(amrGridReader), mValueRange(0.0,0.0)
{
  level = -1;
  data = NULL;
}

AMRGrid::~AMRGrid()
{
  unloadCells();
  for (std_ext::hash_map<const char*, GridProperty*>::iterator it = mGridProperties.begin(); it != mGridProperties.end(); ++it)
    delete it->second;
}

bool AMRGrid::loadCells()
{
  if (datatype!=IObase::Float32) {
    std::cerr << "Error: Grid data is in " << AMRGridAuxFct::dataTypeToStr(IObase::DataType(datatype)) << " format. Currently only Float32 datasets are supported" << std::endl;
    return false;
  }
  if (data==NULL)
    mAmrGridReader->getGridData(*this, mGridIndex);
  return true;
}

void AMRGrid::loadCellsForAllRefiningGrids()
{
  for (AMRRefiningGridIterator refGridIt=refiningGridsBegin(); refGridIt!=refiningGridsEnd(); ++refGridIt)
    (*refGridIt)->loadCells();
}

bool AMRGrid::overlaps(const AMRGrid &otherGrid) const
{
  for (int i=0; i<rank; ++i) {
    int sr1 = spatialrefinement[i];
    int pr1 = gridplacementrefinement[i];
    int sr2 = otherGrid.spatialrefinement[i];
    int pr2 = otherGrid.gridplacementrefinement[i];
    // Increase refinement to sr1*pr1*sr2*pr2
    int min1 = iorigin[i]*sr1*pr2*sr2; // Already refined by pr1
    int min2 = otherGrid.iorigin[i]*pr1*sr1*sr2; // Already refined by pr2
    int ext1 = dims[i]*pr1*pr2*sr2; // Already refined by sr1
    int ext2 = otherGrid.dims[i]*sr1*pr1*pr2; // Already refined by sr2
    if (!intervalsOverlap(min1, min1+ext1-1, min2, min2+ext2-1)) return false;
  }
  return true;
}

void AMRGrid::printInfo(std::ostream &os) const
{
  using namespace AMRGridAuxFct;
  os << "Grid at level " << level << "; origin " << origin << "; dims: " << dims << std::endl;
  os << "Extend: " << minext << " - " << maxext << std::endl;
  os << "Gridspacing: " <<  delta << std::endl;
  os << "Refinement: " << timerefinement << "; " << spatialrefinement << std::endl;
  os << "IOrigin: " << iorigin << " with a placement refinement of " << gridplacementrefinement << std::endl;
  os << "NBytes: " << nbytes << " Datatype: " << dataTypeToStr(IObase::Int2DataType(datatype)) << "; Itemsize: " << IObase::sizeOf(IObase::Int2DataType(datatype)) << std::endl;
}

void AMRGrid::printInfo() const
{
  printInfo(std::cout);
}

AMRGrid::IntersectionMap* AMRGrid::computeIntersectionMap(IntersectionMap::MapType mapType, const PrimitiveBitSet& gridsOfInterest, bool considerAllGrids) const
{
  if (!refinementAvailable() || (!considerAllGrids && gridsOfInterest.empty())) return new IntersectionMap(dims[0], dims[1], dims[2], 1, 1, 1, mapType);
  if (mapType == IntersectionMap::RefinementMap) {
    std::cerr << "Error: Illegal mapType (RefinementMap) in AMRGrid::computeIntersectionMap()" << std::endl;
    return NULL;
  }
  std::list<AMRGrid*> refiningGridsOfInterest;
  int calculationRefinement[3] = { 1, 1, 1 };
  for (AMRRefiningGridConstIterator refiningCandIter=refiningGridsBegin();
      refiningCandIter!=refiningGridsEnd();
      ++refiningCandIter) {
    if (considerAllGrids || gridsOfInterest.getBitValue((*refiningCandIter)->index())) {
      refiningGridsOfInterest.push_back(*refiningCandIter);
      // <ATTENTION> The following choice of calculationRefinement assumes,
      // that spatial refinement and grid placement refinement of all
      // refining grids are integer multiples of each other and of spatial
      // and grid placment refinement of the current grid.
      calculationRefinement[0] = std::max(calculationRefinement[0],
	  std::max((*refiningCandIter)->gridplacementrefinement[0], (*refiningCandIter)->spatialrefinement[0]));
      calculationRefinement[1] = std::max(calculationRefinement[1],
	  std::max((*refiningCandIter)->gridplacementrefinement[1], (*refiningCandIter)->spatialrefinement[1]));
      calculationRefinement[2] = std::max(calculationRefinement[2],
	  std::max((*refiningCandIter)->gridplacementrefinement[2], (*refiningCandIter)->spatialrefinement[2]));
      // </ATTENTION>
    }
  }

  int thisSpatialFactor [3];
  thisSpatialFactor[0] = calculationRefinement[0]/spatialrefinement[0];
  thisSpatialFactor[1] = calculationRefinement[1]/spatialrefinement[1];
  thisSpatialFactor[2] = calculationRefinement[2]/spatialrefinement[2];
  assert(calculationRefinement[0]%spatialrefinement[0]==0);
  assert(calculationRefinement[1]%spatialrefinement[1]==0);
  assert(calculationRefinement[2]%spatialrefinement[2]==0);

  int thisPlacementFactor[3];
  thisPlacementFactor[0] = calculationRefinement[0]/gridplacementrefinement[0];
  thisPlacementFactor[1] = calculationRefinement[1]/gridplacementrefinement[1];
  thisPlacementFactor[2] = calculationRefinement[2]/gridplacementrefinement[2];
  assert(calculationRefinement[0]%gridplacementrefinement[0]==0);
  assert(calculationRefinement[1]%gridplacementrefinement[1]==0);
  assert(calculationRefinement[2]%gridplacementrefinement[2]==0);

  IntersectionMap* theIntersectionMap = new IntersectionMap(dims[0], dims[1], dims[2],
      thisSpatialFactor[0], thisSpatialFactor[1], thisSpatialFactor[2], mapType);

  for (std::list<AMRGrid*>::iterator refiningIter=refiningGridsOfInterest.begin();
      refiningIter!=refiningGridsOfInterest.end();
      ++refiningIter) {
    int otherSpatialFactor [3];
    otherSpatialFactor[0] = calculationRefinement[0]/(*refiningIter)->spatialrefinement[0];
    otherSpatialFactor[1] = calculationRefinement[1]/(*refiningIter)->spatialrefinement[1];
    otherSpatialFactor[2] = calculationRefinement[2]/(*refiningIter)->spatialrefinement[2];
    assert(calculationRefinement[0]%(*refiningIter)->spatialrefinement[0]==0);
    assert(calculationRefinement[1]%(*refiningIter)->spatialrefinement[1]==0);
    assert(calculationRefinement[2]%(*refiningIter)->spatialrefinement[2]==0);
    int otherPlacementFactor[3];
    otherPlacementFactor[0] = calculationRefinement[0]/(*refiningIter)->gridplacementrefinement[0];
    otherPlacementFactor[1] = calculationRefinement[1]/(*refiningIter)->gridplacementrefinement[1];
    otherPlacementFactor[2] = calculationRefinement[2]/(*refiningIter)->gridplacementrefinement[2];
    assert(calculationRefinement[0]%(*refiningIter)->gridplacementrefinement[0]==0);
    assert(calculationRefinement[1]%(*refiningIter)->gridplacementrefinement[1]==0);
    assert(calculationRefinement[2]%(*refiningIter)->gridplacementrefinement[2]==0);

    int minExt[3];
    minExt[0] = (*refiningIter)->iorigin[0]*otherPlacementFactor[0] - iorigin[0]*thisPlacementFactor[0];
    minExt[1] = (*refiningIter)->iorigin[1]*otherPlacementFactor[1] - iorigin[1]*thisPlacementFactor[1];
    minExt[2] = (*refiningIter)->iorigin[2]*otherPlacementFactor[2] - iorigin[2]*thisPlacementFactor[2];
    // Clip
    if (minExt[0] < 0) minExt[0]=0;
    if (minExt[1] < 0) minExt[1]=0;
    if (minExt[2] < 0) minExt[2]=0;

    int maxExt[3];
    maxExt[0] = minExt[0] + (*refiningIter)->dims[0]*otherSpatialFactor[0] - 1;
    maxExt[1] = minExt[1] + (*refiningIter)->dims[1]*otherSpatialFactor[1] - 1;
    maxExt[2] = minExt[2] + (*refiningIter)->dims[2]*otherSpatialFactor[2] - 1;
    // Clip
    if (maxExt[0] > dims[0] * thisSpatialFactor[0] - 1) maxExt[0] = dims[0] * thisSpatialFactor[0] - 1;
    if (maxExt[1] > dims[1] * thisSpatialFactor[1] - 1) maxExt[1] = dims[1] * thisSpatialFactor[1] - 1;
    if (maxExt[2] > dims[2] * thisSpatialFactor[2] - 1) maxExt[2] = dims[2] * thisSpatialFactor[2] - 1;

    if (mapType == IntersectionMap::RefinedInnerCells) {
      minExt[0]++;
      minExt[1]++;
      minExt[2]++;
      maxExt[0]--;
      maxExt[1]--;
      maxExt[2]--;
      if (!(maxExt[0]>minExt[0] && maxExt[1]>minExt[1] && maxExt[2]>minExt[2]))
	// Skip this grid, no inner cells!
	continue;
    }

    int minExtCell[3];
    minExtCell[0] = minExt[0] / thisSpatialFactor[0]; // "Rounds" to floor(x), because >=0 !!!
    minExtCell[1] = minExt[1] / thisSpatialFactor[1]; // "Rounds" to floor(x), because >=0 !!!
    minExtCell[2] = minExt[2] / thisSpatialFactor[2]; // "Rounds" to floor(x), because >=0 !!!
    int minExtSubCell[3];
    minExtSubCell[0] = minExt[0] % thisSpatialFactor[0];
    minExtSubCell[1] = minExt[1] % thisSpatialFactor[1];
    minExtSubCell[2] = minExt[2] % thisSpatialFactor[2];
    int maxExtCell[3];
    maxExtCell[0] = maxExt[0] / thisSpatialFactor[0];
    maxExtCell[1] = maxExt[1] / thisSpatialFactor[1];
    maxExtCell[2] = maxExt[2] / thisSpatialFactor[2];
    int maxExtSubCell[3];
    maxExtSubCell[0] = maxExt[0] % thisSpatialFactor[0];
    maxExtSubCell[1] = maxExt[1] % thisSpatialFactor[1];
    maxExtSubCell[2] = maxExt[2] % thisSpatialFactor[2];
    // Create the low-res intersection information
    if (mapType == IntersectionMap::BorderCells) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	  (*theIntersectionMap)(i,j,minExtCell[2]) = (*refiningIter)->mGridIndex;
	  (*theIntersectionMap)(i,j,maxExtCell[2]) = (*refiningIter)->mGridIndex;
	}
      }
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  (*theIntersectionMap)(i,minExtCell[1],k) = (*refiningIter)->mGridIndex;
	  (*theIntersectionMap)(i,maxExtCell[1],k) = (*refiningIter)->mGridIndex;
	}
      }
      for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  (*theIntersectionMap)(minExtCell[0],j,k) = (*refiningIter)->mGridIndex;
	  (*theIntersectionMap)(maxExtCell[0],j,k) = (*refiningIter)->mGridIndex;
	}
      }
    }
    else if (mapType == IntersectionMap::RefinedCells || mapType == IntersectionMap::RefinedInnerCells) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	  for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	    (*theIntersectionMap)(i,j,k) = (*refiningIter)->mGridIndex;
	  }
	}
      }
    }
    /*
    // Assumption: Each grid covers at least two cells in each extend
    assume(minExtCell[0]!=maxExtCell[0]);
    assume(minExtCell[1]!=maxExtCell[1]);
    assume(minExtCell[2]!=maxExtCell[2]);
     */

    // Create subcells, where necessary
    int subCellIRes = thisSpatialFactor[0];
    int subCellJRes = thisSpatialFactor[1];
    int subCellKRes = thisSpatialFactor[2];

    if (minExtSubCell[2]) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subIMin = (i==minExtCell[0] ? minExtSubCell[0] : 0);
	  int subIMax = (i==maxExtCell[0] ? maxExtSubCell[0] : subCellIRes-1);
	  int subJMin = (j==minExtCell[1] ? minExtSubCell[1] : 0);
	  int subJMax = (j==maxExtCell[1] ? maxExtSubCell[1] : subCellJRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(i, j, minExtCell[2], subIMin, subIMax, subJMin, subJMax, minExtSubCell[2], subCellKRes-1, (*refiningIter)->mGridIndex);
	}
      }
    }
    if (maxExtSubCell[2]!=subCellKRes-1) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subIMin = (i==minExtCell[0] ? minExtSubCell[0] : 0);
	  int subIMax = (i==maxExtCell[0] ? maxExtSubCell[0] : subCellIRes-1);
	  int subJMin = (j==minExtCell[1] ? minExtSubCell[1] : 0);
	  int subJMax = (j==maxExtCell[1] ? maxExtSubCell[1] : subCellJRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(i, j, maxExtCell[2], subIMin, subIMax, subJMin, subJMax, 0, maxExtSubCell[2], (*refiningIter)->mGridIndex);
	}
      }
    }
    if (minExtSubCell[1]) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subIMin = (i==minExtCell[0] ? minExtSubCell[0] : 0);
	  int subIMax = (i==maxExtCell[0] ? maxExtSubCell[0] : subCellIRes-1);
	  int subKMin = (k==minExtCell[2] ? minExtSubCell[2] : 0);
	  int subKMax = (k==maxExtCell[2] ? maxExtSubCell[2] : subCellKRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(i, minExtCell[1], k, subIMin, subIMax, minExtSubCell[1], subCellJRes-1, subKMin, subKMax, (*refiningIter)->mGridIndex);
	}
      }
    }
    if (maxExtSubCell[1]!=subCellJRes-1) {
      for (int i=minExtCell[0]; i<=maxExtCell[0]; ++i) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subIMin = (i==minExtCell[0] ? minExtSubCell[0] : 0);
	  int subIMax = (i==maxExtCell[0] ? maxExtSubCell[0] : subCellIRes-1);
	  int subKMin = (k==minExtCell[2] ? minExtSubCell[2] : 0);
	  int subKMax = (k==maxExtCell[2] ? maxExtSubCell[2] : subCellKRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(i, maxExtCell[1], k, subIMin, subIMax, 0, maxExtSubCell[1], subKMin, subKMax, (*refiningIter)->mGridIndex);
	}
      }
    }
    if (minExtSubCell[0]) {
      for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subJMin = (j==minExtCell[1] ? minExtSubCell[1] : 0);
	  int subJMax = (j==maxExtCell[1] ? maxExtSubCell[1] : subCellJRes-1);
	  int subKMin = (k==minExtCell[2] ? minExtSubCell[2] : 0);
	  int subKMax = (k==maxExtCell[2] ? maxExtSubCell[2] : subCellKRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(minExtCell[0], j, k, minExtSubCell[0], subCellIRes-1, subJMin, subJMax, subKMin, subKMax, (*refiningIter)->mGridIndex);
	}
      }
    }
    if (maxExtSubCell[0]!=subCellIRes-1) {
      for (int j=minExtCell[1]; j<=maxExtCell[1]; ++j) {
	for (int k=minExtCell[2]; k<=maxExtCell[2]; ++k) {
	  // <HACK>
	  // This isn't nice but avoids checking for serveral special cases
	  int subJMin = (j==minExtCell[1] ? minExtSubCell[1] : 0);
	  int subJMax = (j==maxExtCell[1] ? maxExtSubCell[1] : subCellJRes-1);
	  int subKMin = (k==minExtCell[2] ? minExtSubCell[2] : 0);
	  int subKMax = (k==maxExtCell[2] ? maxExtSubCell[2] : subCellKRes-1);
	  // </HACK>
	  theIntersectionMap->createBlockRefined(maxExtCell[0], j, k, 0, maxExtSubCell[0], subJMin, subJMax, subKMin, subKMax, (*refiningIter)->mGridIndex);
	}
      }
    }
  }
  return theIntersectionMap;
}

AMRGrid::AMRRefinementInfo *AMRGrid::getRefinementInfo(const AMRGrid *refiningGrid) const
{
  if(level != refiningGrid->level-1) return 0;

  int calculationRefinement[3];
  calculationRefinement[0] = std::min(refiningGrid->spatialrefinement[0], refiningGrid->gridplacementrefinement[0]);
  calculationRefinement[1] = std::min(refiningGrid->spatialrefinement[1], refiningGrid->gridplacementrefinement[1]);
  calculationRefinement[2] = std::min(refiningGrid->spatialrefinement[2], refiningGrid->gridplacementrefinement[2]);

  // Assume that gridplacementrefinement and spatialrefinement can be expressed
  // as integer ratios with respect to each other and the factors of the parent level
  assume(calculationRefinement[0] % refiningGrid->spatialrefinement[0] == 0); 
  assume(calculationRefinement[1] % refiningGrid->spatialrefinement[1] == 0); 
  assume(calculationRefinement[2] % refiningGrid->spatialrefinement[2] == 0); 
  assume(calculationRefinement[0] % refiningGrid->gridplacementrefinement[0] == 0);
  assume(calculationRefinement[1] % refiningGrid->gridplacementrefinement[1] == 0);
  assume(calculationRefinement[2] % refiningGrid->gridplacementrefinement[2] == 0);
  assume(calculationRefinement[0] % spatialrefinement[0] == 0); 
  assume(calculationRefinement[1] % spatialrefinement[1] == 0); 
  assume(calculationRefinement[2] % spatialrefinement[2] == 0); 
  assume(calculationRefinement[0] % gridplacementrefinement[0] == 0);
  assume(calculationRefinement[1] % gridplacementrefinement[1] == 0);
  assume(calculationRefinement[2] % gridplacementrefinement[2] == 0);

  int refiningPlacementFactor[3];
  refiningPlacementFactor[0] = calculationRefinement[0] / refiningGrid->gridplacementrefinement[0];
  refiningPlacementFactor[1] = calculationRefinement[1] / refiningGrid->gridplacementrefinement[1];
  refiningPlacementFactor[2] = calculationRefinement[2] / refiningGrid->gridplacementrefinement[2];
  int refiningSpatialFactor[3];
  refiningSpatialFactor[0] =  calculationRefinement[0] / refiningGrid->spatialrefinement[0];
  refiningSpatialFactor[1] =  calculationRefinement[1] / refiningGrid->spatialrefinement[1];
  refiningSpatialFactor[2] =  calculationRefinement[2] / refiningGrid->spatialrefinement[2];
  int refinedPlacementFactor[3];
  refinedPlacementFactor[0] = calculationRefinement[0] / gridplacementrefinement[0];
  refinedPlacementFactor[1] = calculationRefinement[1] / gridplacementrefinement[1];
  refinedPlacementFactor[2] = calculationRefinement[2] / gridplacementrefinement[2];
  int refinedSpatialFactor[3];
  refinedSpatialFactor[0] =  calculationRefinement[0] / spatialrefinement[0];
  refinedSpatialFactor[1] =  calculationRefinement[1] / spatialrefinement[1];
  refinedSpatialFactor[2] =  calculationRefinement[2] / spatialrefinement[2];
  int minExt[3];
  minExt[0] = refiningGrid->iorigin[0]*refiningPlacementFactor[0] - iorigin[0]*refinedPlacementFactor[0];
  minExt[1] = refiningGrid->iorigin[1]*refiningPlacementFactor[1] - iorigin[1]*refinedPlacementFactor[1];
  minExt[2] = refiningGrid->iorigin[2]*refiningPlacementFactor[2] - iorigin[2]*refinedPlacementFactor[2];
  int maxExt[3];
  maxExt[0] = minExt[0]+(refiningGrid->dims[0]-1)*refiningSpatialFactor[0];
  maxExt[1] = minExt[1]+(refiningGrid->dims[1]-1)*refiningSpatialFactor[1];
  maxExt[2] = minExt[2]+(refiningGrid->dims[2]-1)*refiningSpatialFactor[2];

  int refinedCellsMin[3];
  refinedCellsMin[0] = minExt[0] / refinedSpatialFactor[0];
  refinedCellsMin[1] = minExt[1] / refinedSpatialFactor[1];
  refinedCellsMin[2] = minExt[2] / refinedSpatialFactor[2];
  int refinedCellsSubMin[3];
  refinedCellsSubMin[0] = minExt[0] % refinedSpatialFactor[0];
  refinedCellsSubMin[1] = minExt[1] % refinedSpatialFactor[1];
  refinedCellsSubMin[2] = minExt[2] % refinedSpatialFactor[2];
  int refinedCellsMax[3];
  refinedCellsMax[0] = maxExt[0] / refinedSpatialFactor[0];
  refinedCellsMax[1] = maxExt[1] / refinedSpatialFactor[1];
  refinedCellsMax[2] = maxExt[2] / refinedSpatialFactor[2];
  int refinedCellsSubMax[3];
  refinedCellsSubMax[0] = maxExt[0] % refinedSpatialFactor[0];
  refinedCellsSubMax[1] = maxExt[1] % refinedSpatialFactor[1];
  refinedCellsSubMax[2] = maxExt[2] % refinedSpatialFactor[2];

  int completelyRefinedCellsMin[3];
  completelyRefinedCellsMin[0] = refinedCellsMin[0] + (refinedCellsSubMin[0] == 0 ? 0 : 1);
  completelyRefinedCellsMin[1] = refinedCellsMin[1] + (refinedCellsSubMin[1] == 0 ? 0 : 1);
  completelyRefinedCellsMin[2] = refinedCellsMin[2] + (refinedCellsSubMin[2] == 0 ? 0 : 1);
  int completelyRefinedCellsMax[3];
  completelyRefinedCellsMax[0] = refinedCellsMax[0] - (refinedCellsSubMax[0] != refinedSpatialFactor[0]-1 ? 1 : 0);
  completelyRefinedCellsMax[1] = refinedCellsMax[1] - (refinedCellsSubMax[1] != refinedSpatialFactor[1]-1 ? 1 : 0);
  completelyRefinedCellsMax[2] = refinedCellsMax[2] - (refinedCellsSubMax[2] != refinedSpatialFactor[2]-1 ? 1 : 0);

  int refiningCellsMin[3];
  refiningCellsMin[0] = (completelyRefinedCellsMin[0] * refinedSpatialFactor[0] - minExt[0]) / refiningSpatialFactor[0];
  refiningCellsMin[1] = (completelyRefinedCellsMin[1] * refinedSpatialFactor[1] - minExt[1]) / refiningSpatialFactor[1];
  refiningCellsMin[2] = (completelyRefinedCellsMin[2] * refinedSpatialFactor[2] - minExt[2]) / refiningSpatialFactor[2];
  assume((completelyRefinedCellsMin[0] * refinedSpatialFactor[0] - minExt[0]) % refiningSpatialFactor[0] == 0);
  assume((completelyRefinedCellsMin[1] * refinedSpatialFactor[1] - minExt[1]) % refiningSpatialFactor[1] == 0);
  assume((completelyRefinedCellsMin[2] * refinedSpatialFactor[2] - minExt[2]) % refiningSpatialFactor[2] == 0);
  int refiningCellsMax[3];
  refiningCellsMax[0] = ((completelyRefinedCellsMax[0]+1)*refinedSpatialFactor[0]-minExt[0]) / refiningSpatialFactor[0] - 1;
  refiningCellsMax[1] = ((completelyRefinedCellsMax[1]+1)*refinedSpatialFactor[1]-minExt[1]) / refiningSpatialFactor[1] - 1;
  refiningCellsMax[2] = ((completelyRefinedCellsMax[2]+1)*refinedSpatialFactor[2]-minExt[2]) / refiningSpatialFactor[2] - 1;
  assume((completelyRefinedCellsMax[0]*refinedSpatialFactor[0] - minExt[0] + 1) % refiningSpatialFactor[0] == 0);
  assume((completelyRefinedCellsMax[1]*refinedSpatialFactor[1] - minExt[1] + 1) % refiningSpatialFactor[1] == 0);
  assume((completelyRefinedCellsMax[2]*refinedSpatialFactor[2] - minExt[2] + 1) % refiningSpatialFactor[2] == 0);

  int mSpatialRefinementWRTParent[3];
  mSpatialRefinementWRTParent[0] = refiningGrid->spatialrefinement[0]/spatialrefinement[0];
  mSpatialRefinementWRTParent[1] = refiningGrid->spatialrefinement[1]/spatialrefinement[1];
  mSpatialRefinementWRTParent[2] = refiningGrid->spatialrefinement[2]/spatialrefinement[2];
  assume(refiningGrid->spatialrefinement[0]%spatialrefinement[0]==0);
  assume(refiningGrid->spatialrefinement[1]%spatialrefinement[1]==0);
  assume(refiningGrid->spatialrefinement[2]%spatialrefinement[2]==0);

  return new AMRRefinementInfo(
      AxisAlignedBox<int>(completelyRefinedCellsMin, completelyRefinedCellsMax),
      AxisAlignedBox<int>(refiningCellsMin, refiningCellsMax),
      AxisAlignedBox<int>(refinedCellsMin, refinedCellsMax),
      AxisAlignedBox<int>(refinedCellsSubMin, refinedCellsSubMax),
      mSpatialRefinementWRTParent
			      );
}

AxisAlignedBox<int> AMRGrid::getRefiningCells(const AMRGrid *otherGrid, int i, int j, int k) const
{
  if (otherGrid->level > level) return otherGrid->getRefiningCells(this, i, j, k);
  assume(spatialrefinement[0] == gridplacementrefinement[0]);
  assume(spatialrefinement[1] == gridplacementrefinement[1]);
  assume(spatialrefinement[2] == gridplacementrefinement[2]);
  int relativeRefinement[3];
  relativeRefinement[0] = spatialrefinement[0] / otherGrid->spatialrefinement[0];
  relativeRefinement[1] = spatialrefinement[1] / otherGrid->spatialrefinement[1];
  relativeRefinement[2] = spatialrefinement[2] / otherGrid->spatialrefinement[2];
  assume(spatialrefinement[0] % otherGrid->spatialrefinement[0] == 0);
  assume(spatialrefinement[1] % otherGrid->spatialrefinement[1] == 0);
  assume(spatialrefinement[2] % otherGrid->spatialrefinement[2] == 0);
  int start[3];
  start[0] = (otherGrid->iorigin[0]+i) * relativeRefinement[0] - iorigin[0];
  start[1] = (otherGrid->iorigin[1]+j) * relativeRefinement[1] - iorigin[1];
  start[2] = (otherGrid->iorigin[2]+k) * relativeRefinement[2] - iorigin[2];
  int end[3];
  end[0] = start[0] + relativeRefinement[0] - 1;
  end[1] = start[1] + relativeRefinement[1] - 1;
  end[2] = start[2] + relativeRefinement[2] - 1;
  return AxisAlignedBox<int>(start, end);
}

AMRHierarchy::GridId AMRGrid::getCellRefiningGridId(int i, int j, int k) const
{
  // Get handle for next level
  AMRLevel *nextLevel = mAMRHierarchy->getLevel(level+1);

  if (nextLevel) {
    // Convert into absoulte level coordinates of next level
    // <ATTENTION> This assumes spatialrefinment == gridplacementrefinment </ATTENTION>
    int absIdxInNextLevel[3];
    absIdxInNextLevel[xAxis] = (i + iorigin[xAxis])*(nextLevel->getSpatialRefinement(xAxis)/spatialrefinement[xAxis]);
    absIdxInNextLevel[yAxis] = (j + iorigin[yAxis])*(nextLevel->getSpatialRefinement(yAxis)/spatialrefinement[yAxis]);
    absIdxInNextLevel[zAxis] = (k + iorigin[zAxis])*(nextLevel->getSpatialRefinement(zAxis)/spatialrefinement[zAxis]);

    return nextLevel->getGridId(absIdxInNextLevel);
  }
  else 
    return AMRGrid::NoGrid;
}

void AMRGrid::addProperty(const char* name, AMRGrid::GridProperty* property)
{
  mGridProperties[name] = property;
}

AMRGrid::GridProperty* AMRGrid::getProperty(const char *name) const
{
  std_ext::hash_map<const char*, GridProperty*>::const_iterator pos = mGridProperties.find(name);
  if (pos == mGridProperties.end())
    return 0;
  else
    return pos->second;
}

void AMRGrid::removeProperty(const char* name)
{
  mGridProperties.erase(name);
}

void AMRGrid::computeValueRange()
{
  bool cellsAlreadyLoaded = (data != 0);
  if (loadCells()) {
    assert(data);
    mValueRange.min = getValue(0);
    mValueRange.max = mValueRange.min;
    for (int i = 0; i < dims[0] * dims[1] * dims[2]; ++i) {
      float val = getValue(i);
      if ( val < mValueRange.min) 
	mValueRange.min = val;
      else if (val > mValueRange.max) 
	mValueRange.max = val;
    }
  }
  else {
    std::cerr << "Error: Couldn't calculate value range. No data!" << std::endl;
  }
  if (!cellsAlreadyLoaded) unloadCells();
}