Commit 602d9b65 authored by BRAMAS Berenger's avatar BRAMAS Berenger

Add a new way to find neighbor/interaction morton indexes in nearly constant time

parent 31637f5b
#ifndef FNEIGHBORINDEXES_HPP
#define FNEIGHBORINDEXES_HPP
#include "../Utils/FGlobal.hpp"
/**
* @brief The FAbstractNeighborIndex class is an abstraction of
* what should propose classes that can find the neighbor morton indexes
* from a given index.
* These kind of classes are useful in the P2P and the M2L
* (in the find neighbors functions)
*/
class FAbstractNeighborIndex{
public:
/* will return -1 or 0 if the cell is at the border */
virtual int minX() = 0;
/* will return -1 or 0 if the cell is at the border */
virtual int minY() = 0;
/* will return -1 or 0 if the cell is at the border */
virtual int minZ() = 0;
/* will return 1 or 0 if the cell is at the border */
virtual int maxX() = 0;
/* will return 1 or 0 if the cell is at the border */
virtual int maxY() = 0;
/* will return 1 or 0 if the cell is at the border */
virtual int maxZ() = 0;
/* given the neighbor relative position returns the morton index */
virtual MortonIndex getIndex(const int inX, const int inY, const int inZ) = 0;
};
#include "FTreeCoordinate.hpp"
/**
* This class compute the neigh position from the tree coordinates
*/
class FCoordinateNeighborIndex : public FAbstractNeighborIndex{
int currentMinX;
int currentMinY;
int currentMinZ;
int currentMaxX;
int currentMaxY;
int currentMaxZ;
const MortonIndex mindex;
const int level;
const FTreeCoordinate coord;
public:
FCoordinateNeighborIndex(const MortonIndex inMindex, const int inLevel)
: mindex(inMindex), level(inLevel), coord(mindex, level) {
currentMinX = (coord.getX()==0? 0 : -1);
currentMinY = (coord.getY()==0? 0 : -1);
currentMinZ = (coord.getZ()==0? 0 : -1);
const int limite = FMath::pow2(level);
currentMaxX = (coord.getX()== (limite-1)? 0 : 1);
currentMaxY = (coord.getY()== (limite-1)? 0 : 1);
currentMaxZ = (coord.getZ()== (limite-1)? 0 : 1);
}
int minX() override{
return currentMinX;
}
int minY() override{
return currentMinY;
}
int minZ() override{
return currentMinZ;
}
int maxX() override{
return currentMaxX;
}
int maxY() override{
return currentMaxY;
}
int maxZ() override{
return currentMaxZ;
}
virtual MortonIndex getIndex(const int inX, const int inY, const int inZ) override{
return FTreeCoordinate(inX+coord.getX(), inY+coord.getY(), inZ+coord.getZ()).getMortonIndex(level);
}
};
/**
* This class returns the neigh indexes from bitwise operations.
*/
class FBitsNeighborIndex : public FAbstractNeighborIndex{
int currentMinX;
int currentMinY;
int currentMinZ;
int currentMaxX;
int currentMaxY;
int currentMaxZ;
const MortonIndex mindex;
const int level;
const MortonIndex flagZ;
const MortonIndex flagY;
const MortonIndex flagX;
MortonIndex mindexes[3];
static const int NbBitsFlag = MaxTreeHeight;
static MortonIndex BuildFlag(){
MortonIndex flag = 0;
int counter = NbBitsFlag;
while(counter--){
flag = (flag<<3) | 1;
}
return flag;
}
public:
FBitsNeighborIndex(const MortonIndex inMindex, const int inLevel)
: mindex(inMindex), level(inLevel),
flagZ(BuildFlag()), flagY(flagZ<<1), flagX(flagZ<<2){
currentMinX = (((mindex&flagX)^flagX) == flagX ? 0 : -1);
currentMinY = (((mindex&flagY)^flagY) == flagY ? 0 : -1);
currentMinZ = (((mindex&flagZ)^flagZ) == flagZ ? 0 : -1);
const int nbShiftLimit = (NbBitsFlag-level)*3;
const MortonIndex limiteX = (flagX>>nbShiftLimit);
const MortonIndex limiteY = (flagY>>nbShiftLimit);
const MortonIndex limiteZ = (flagZ>>nbShiftLimit);
currentMaxX = ((mindex&limiteX) == limiteX? 0 : 1);
currentMaxY = ((mindex&limiteY) == limiteY? 0 : 1);
currentMaxZ = ((mindex&limiteZ) == limiteZ? 0 : 1);
{
MortonIndex mindex_minus = mindex;
MortonIndex flag_minus = (((mindex&flagX)^flagX) == flagX ? 0 : 4)
| (((mindex&flagY)^flagY) == flagY ? 0 : 2)
| (((mindex&flagZ)^flagZ) == flagZ ? 0 : 1);
while(flag_minus){
const MortonIndex prevflag_minus = flag_minus;
flag_minus = (flag_minus& ~mindex_minus)<<3;
mindex_minus = (mindex_minus^prevflag_minus);
}
mindexes[0] = mindex_minus;
}
{
MortonIndex mindex_plus = mindex;
MortonIndex flag_plus = ((mindex&limiteX) == limiteX? 0 : 4)
| ((mindex&limiteY) == limiteY? 0 : 2)
| ((mindex&limiteZ) == limiteZ? 0 : 1);
while(flag_plus){
const MortonIndex prevflag_plus = flag_plus;
flag_plus = (flag_plus& mindex_plus)<<3;
mindex_plus = (mindex_plus^prevflag_plus);
}
mindexes[2] = mindex_plus;
}
mindexes[1] = mindex;
}
int minX() override{
return currentMinX;
}
int minY() override{
return currentMinY;
}
int minZ() override{
return currentMinZ;
}
int maxX() override{
return currentMaxX;
}
int maxY() override{
return currentMaxY;
}
int maxZ() override{
return currentMaxZ;
}
virtual MortonIndex getIndex(const int inX, const int inY, const int inZ) override{
return (mindexes[inX+1]&flagX) | (mindexes[inY+1]&flagY) | (mindexes[inZ+1]&flagZ);
}
};
#endif // FNEIGHBORINDEXES_HPP
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#include "FUTester.hpp"
#include "Containers/FTreeCoordinate.hpp"
#include "Containers/FNeighborIndexes.hpp"
/**
* This file is a unit test for the FNeigborIndexes classes
*/
/** this class test the list container */
class TestIndexes : public FUTester<TestIndexes> {
void MortonLimite(){
for(int idxLevel = 1 ; idxLevel < MaxTreeHeight ; ++idxLevel){
const int limit = FMath::pow2(idxLevel);
for(int x = 0 ; x < 3 ; ++x){
const int xbox = (x*(limit-1))/2;
for(int y = 0 ; y < 3 ; ++y){
const int ybox = (y*(limit-1))/2;
for(int z = 0 ; z < 3 ; ++z){
const int zbox = (z*(limit-1))/2;
const MortonIndex mindex = FTreeCoordinate(xbox, ybox, zbox).getMortonIndex(idxLevel);
FCoordinateNeighborIndex coordindexes(mindex, idxLevel);
FBitsNeighborIndex bitsindex(mindex, idxLevel);
uassert(coordindexes.minX() == bitsindex.minX());
uassert(coordindexes.minY() == bitsindex.minY());
uassert(coordindexes.minZ() == bitsindex.minZ());
uassert(coordindexes.maxX() == bitsindex.maxX());
uassert(coordindexes.maxY() == bitsindex.maxY());
uassert(coordindexes.maxZ() == bitsindex.maxZ());
for(int idxX = coordindexes.minX() ; idxX <= coordindexes.maxX() ; ++idxX){
for(int idxY = coordindexes.minY() ; idxY <= coordindexes.maxY() ; ++idxY){
for(int idxZ = coordindexes.minZ() ; idxZ <= coordindexes.maxZ() ; ++idxZ){
if(idxX || idxY || idxZ){
uassert(coordindexes.getIndex(idxX, idxY, idxZ)
== bitsindex.getIndex(idxX, idxY, idxZ));
const MortonIndex neigh = bitsindex.getIndex(idxX, idxY, idxZ);
const FTreeCoordinate neighCoord(neigh, idxLevel);
uassert(xbox + idxX == neighCoord.getX());
uassert(ybox + idxY == neighCoord.getY());
uassert(zbox + idxZ == neighCoord.getZ());
}
}
}
}
}
}
}
}
}
void Morton(){
srand(0);
for(int idxLevel = 1 ; idxLevel < MaxTreeHeight ; ++idxLevel){
const int limit = FMath::pow2(idxLevel);
for(int idxTest = 0 ; idxTest < 100 ; ++idxTest){
const int xbox = int(drand48()*double(limit));
const int ybox = int(drand48()*double(limit));
const int zbox = int(drand48()*double(limit));
const MortonIndex mindex = FTreeCoordinate(xbox, ybox, zbox).getMortonIndex(idxLevel);
FCoordinateNeighborIndex coordindexes(mindex, idxLevel);
FBitsNeighborIndex bitsindex(mindex, idxLevel);
uassert(coordindexes.minX() == bitsindex.minX());
uassert(coordindexes.minY() == bitsindex.minY());
uassert(coordindexes.minZ() == bitsindex.minZ());
uassert(coordindexes.maxX() == bitsindex.maxX());
uassert(coordindexes.maxY() == bitsindex.maxY());
uassert(coordindexes.maxZ() == bitsindex.maxZ());
for(int idxX = coordindexes.minX() ; idxX <= coordindexes.maxX() ; ++idxX){
for(int idxY = coordindexes.minY() ; idxY <= coordindexes.maxY() ; ++idxY){
for(int idxZ = coordindexes.minZ() ; idxZ <= coordindexes.maxZ() ; ++idxZ){
if(idxX || idxY || idxZ){
uassert(coordindexes.getIndex(idxX, idxY, idxZ)
== bitsindex.getIndex(idxX, idxY, idxZ));
const MortonIndex neigh = bitsindex.getIndex(idxX, idxY, idxZ);
const FTreeCoordinate neighCoord(neigh, idxLevel);
uassert(xbox + idxX == neighCoord.getX());
uassert(ybox + idxY == neighCoord.getY());
uassert(zbox + idxZ == neighCoord.getZ());
}
}
}
}
}
}
}
void MortonAll(){
{
const int idxLevel = 5;
const int limit = FMath::pow2(idxLevel);
for(int xbox = 0 ; xbox < limit ; ++xbox){
for(int ybox = 0 ; ybox < limit ; ++ybox){
for(int zbox = 0 ; zbox < limit ; ++zbox){
const MortonIndex mindex = FTreeCoordinate(xbox, ybox, zbox).getMortonIndex(idxLevel);
FCoordinateNeighborIndex coordindexes(mindex, idxLevel);
FBitsNeighborIndex bitsindex(mindex, idxLevel);
uassert(coordindexes.minX() == bitsindex.minX());
uassert(coordindexes.minY() == bitsindex.minY());
uassert(coordindexes.minZ() == bitsindex.minZ());
uassert(coordindexes.maxX() == bitsindex.maxX());
uassert(coordindexes.maxY() == bitsindex.maxY());
uassert(coordindexes.maxZ() == bitsindex.maxZ());
for(int idxX = coordindexes.minX() ; idxX <= coordindexes.maxX() ; ++idxX){
for(int idxY = coordindexes.minY() ; idxY <= coordindexes.maxY() ; ++idxY){
for(int idxZ = coordindexes.minZ() ; idxZ <= coordindexes.maxZ() ; ++idxZ){
if(idxX || idxY || idxZ){
uassert(coordindexes.getIndex(idxX, idxY, idxZ)
== bitsindex.getIndex(idxX, idxY, idxZ));
const MortonIndex neigh = bitsindex.getIndex(idxX, idxY, idxZ);
const FTreeCoordinate neighCoord(neigh, idxLevel);
uassert(xbox + idxX == neighCoord.getX());
uassert(ybox + idxY == neighCoord.getY());
uassert(zbox + idxZ == neighCoord.getZ());
}
}
}
}
}
}
}
}
}
// set test
void SetTests(){
AddTest(&TestIndexes::MortonLimite,"Test NeighborIndexes at limits");
AddTest(&TestIndexes::Morton,"Test NeighborIndexes");
AddTest(&TestIndexes::MortonAll,"Test All NeighborIndexes");
}
};
// You must do this
TestClass(TestIndexes)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment