Interpolation.h

Go to the documentation of this file.

#pragma once
 
 
#include "objects/finders/NeighbourFinder.h"
#include "objects/wrappers/SharedPtr.h"
#include "quantities/Storage.h"
#include "sph/kernel/Kernel.h"
#include "system/Factory.h"
#include "system/Settings.h"
#include "thread/Scheduler.h"
 
NAMESPACE_SPH_BEGIN
 
template <typename Type>
class IInterpolant : public Noncopyable {
public:
    virtual Type interpolate(const Vector& pos) const = 0;
};
 
template <typename Type>
class SphInterpolant : public IInterpolant<Type> {
protected:
    ArrayView<const Vector> r;
    ArrayView<const Float> m, rho;
    ArrayView<const Type> quantity;
 
    AutoPtr<IBasicFinder> finder;
    LutKernel<3> kernel;
    Float searchRadius = 0._f;
 
public:
    SphInterpolant(const Storage& storage,
        const QuantityId id,
        const OrderEnum deriv,
        const RunSettings& settings = RunSettings::getDefaults())
        : SphInterpolant(storage, id, deriv, Factory::getFinder(settings), Factory::getKernel<3>(settings)) {}
 
    SphInterpolant(const Storage& storage,
        const QuantityId id,
        const OrderEnum deriv,
        AutoPtr<IBasicFinder>&& finder,
        LutKernel<3>&& kernel)
        : finder(std::move(finder))
        , kernel(std::move(kernel)) {
        SPH_ASSERT(this->finder);
        r = storage.getValue<Vector>(QuantityId::POSITION);
        tie(m, rho) = storage.getValues<Float>(QuantityId::MASS, QuantityId::DENSITY);
        quantity = storage.getAll<Type>(id)[int(deriv)];
        this->build(storage);
    }
 
    virtual Type interpolate(const Vector& pos) const override {
        Type interpol(0._f);
        Array<NeighbourRecord> neighs;
        finder->findAll(pos, searchRadius, neighs);
        for (Size k = 0; k < neighs.size(); ++k) {
            const Size j = neighs[k].index;
            if (neighs[k].distanceSqr < sqr(kernel.radius() * r[j][H])) {
                interpol += m[j] / rho[j] * quantity[j] * kernel.value(pos - r[j], r[j][H]);
            }
        }
        return interpol;
    }
 
private:
    void build(const Storage& storage) {
        ArrayView<const Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
        finder->build(SEQUENTIAL, r);
        searchRadius = 0._f;
        for (Size i = 0; i < r.size(); ++i) {
            searchRadius = max(searchRadius, kernel.radius() * r[i][H]);
        }
    }
};
 
template <typename Type>
class CorrectedSphInterpolant : public SphInterpolant<Type> {
public:
    using SphInterpolant<Type>::SphInterpolant;
 
    virtual Type interpolate(const Vector& pos) const override {
        Type interpol(0._f);
        Float weight = 0._f;
        Array<NeighbourRecord> neighs;
        this->finder->findAll(pos, this->searchRadius, neighs);
        for (Size k = 0; k < neighs.size(); ++k) {
            const Size j = neighs[k].index;
            if (neighs[k].distanceSqr < sqr(this->kernel.radius() * this->r[j][H])) {
                const Float w =
                    this->m[j] / this->rho[j] * this->kernel.value(pos - this->r[j], this->r[j][H]);
                interpol += this->quantity[j] * w;
                weight += w;
            }
        }
        if (weight > 0._f) {
            return interpol / weight;
        } else {
            return Type(0._f);
        }
    }
};
 
NAMESPACE_SPH_END