LookupMap.h

Go to the documentation of this file.

#pragma once
 
 
#include "objects/containers/Array.h"
#include "objects/geometry/Box.h"
#include "objects/geometry/Vector.h"
 
NAMESPACE_SPH_BEGIN
 
class LookupMap : public Noncopyable {
private:
    Array<Array<Size>> storage;
 
    Box tightBox;
    Box boundingBox;
 
    Size dimensionSize;
 
public:
    LookupMap() = default;
 
    LookupMap(const Size n)
        : storage(pow<3>(n))
        , dimensionSize(n) {}
 
    void update(ArrayView<const Vector> points) {
        tightBox = Box();
        for (const Vector& v : points) {
            tightBox.extend(v);
        }
        boundingBox = this->extendBox(tightBox, 1.e-6_f);
        for (auto& a : storage) {
            a.clear();
        }
        // put particles into voxels
        for (Size i = 0; i < points.size(); ++i) {
            Indices idxs = this->map(points[i]);
            (*this)(idxs).push(i);
        }
    }
 
    INLINE bool empty() const {
        return storage.empty();
    }
 
    INLINE Vector clamp(const Vector& pos) const {
        return tightBox.clamp(pos);
    }
 
    LookupMap& operator=(LookupMap&& other) {
        storage = std::move(other.storage);
        dimensionSize = other.dimensionSize;
        boundingBox = other.boundingBox;
        return *this;
    }
 
    INLINE const Array<Size>& operator()(const Indices& idxs) const {
        const Size idx = idxs[X] * sqr(dimensionSize) + idxs[Y] * dimensionSize + idxs[Z];
        SPH_ASSERT(unsigned(idx) < unsigned(storage.size()));
        return storage[idx];
    }
 
    INLINE Array<Size>& operator()(const Indices& idxs) {
        const Size idx = idxs[X] * sqr(dimensionSize) + idxs[Y] * dimensionSize + idxs[Z];
        SPH_ASSERT(unsigned(idx) < unsigned(storage.size()));
        return storage[idx];
    }
 
    INLINE Box voxel(const Indices idxs) const {
        Vector lower(boundingBox.lower() + boundingBox.size() * Vector(idxs) / dimensionSize);
        Vector upper(boundingBox.lower() + boundingBox.size() * Vector(idxs + Indices(1)) / dimensionSize);
        return Box(lower, upper);
    }
 
    INLINE Vector getVoxelSize() const {
        return boundingBox.size() / dimensionSize;
    }
 
    INLINE Size getDimensionSize() const {
        return dimensionSize;
    }
 
    INLINE Indices map(const Vector& v) const {
        SPH_ASSERT(boundingBox.size()[X] > 0._f && boundingBox.size()[Y] > 0._f && boundingBox.size()[Z] > 0._f);
        SPH_ASSERT(dimensionSize >= 1);
        Vector idxs = (v - boundingBox.lower()) / (boundingBox.size()) * dimensionSize;
        // Ordinarily, idxs are never <0 or >=dimensionSize, BUT in case the points are slightly move, this
        // can happen. In this case, we want to return a valid result rather than crashing, even though some
        // neighbours might be missed if the tree is not rebuilt
        return Indices(Sph::clamp(idxs, Vector(0._f), Vector(dimensionSize - 1._f)));
    }
 
private:
    Box extendBox(const Box& box, const Float eps) const {
        const Vector extension =
            max(eps * box.size(), eps * abs(box.lower()), eps * abs(box.upper()), Vector(eps));
        Box extendedBox = box;
        extendedBox.extend(box.upper() + extension);
        extendedBox.extend(box.lower() - extension);
        return extendedBox;
    }
};
 
NAMESPACE_SPH_END