diff options
Diffstat (limited to 'src/AMRHierLib/AMRGrid.cc')
| -rw-r--r-- | src/AMRHierLib/AMRGrid.cc | 612 |
1 files changed, 612 insertions, 0 deletions
diff --git a/src/AMRHierLib/AMRGrid.cc b/src/AMRHierLib/AMRGrid.cc new file mode 100644 index 0000000..8e429b9 --- /dev/null +++ b/src/AMRHierLib/AMRGrid.cc @@ -0,0 +1,612 @@ +//------------------------------------------------------------------------------ +// +// 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(); +} |