rbf.hpp

/*---------------------------------------------------------------------------*\
*
*  bitpit
*
*  Copyright (C) 2015-2021 OPTIMAD engineering Srl
*
*  -------------------------------------------------------------------------
*  License
*  This file is part of bitpit.
*
*  bitpit is free software: you can redistribute it and/or modify it
*  under the terms of the GNU Lesser General Public License v3 (LGPL)
*  as published by the Free Software Foundation.
*
*  bitpit 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 Lesser General Public
*  License for more details.
*
*  You should have received a copy of the GNU Lesser General Public License
*  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
*
\*---------------------------------------------------------------------------*/
 
#ifndef __BITPIT_RBF_HPP__
#define __BITPIT_RBF_HPP__
 
#include "bitpit_discretization.hpp"
#include <array>
#include <set>
#include <vector>
 
namespace bitpit{
 
enum class RBFBasisFunction {
    CUSTOM     = 0,  
    WENDLANDC2 = 1,  
    LINEAR     = 2,  
    GAUSS90    = 3,  
    GAUSS95    = 4,  
    GAUSS99    = 5,  
    C1C0       = 6,  
    C2C0       = 7,  
    C0C1       = 8,  
    C1C1       = 9,  
    C2C1       = 10, 
    C0C2       = 11, 
    C1C2       = 12, 
    C2C2       = 13, 
    COSINUS    = 14, 
    THINPLATE  = 15, 
};
 
enum class RBFMode {
    INTERP = 1, 
    PARAM  = 2  
};
 
class RBFKernel{
    class LinearPolynomial : public ReconstructionPolynomial
    {
private:
        int m_dim;                          
        int m_fields;                       
    public:
        LinearPolynomial();
        void clear();
        void setDimension(int dim);
        void setDataCount(int fields);
        void evalBasis(const double *x, double *basis);
        void initialize();
    };
 
private:
    int     m_fields;                               
    RBFMode m_mode;                                 
    std::vector<double>  m_supportRadii;            
    RBFBasisFunction m_typef;                       
    double  (*m_fPtr)(double);
 
    std::vector<double>                 m_error;    
protected:
    std::vector<std::vector<double>>    m_value;    
    std::vector<std::vector<double>>    m_weight;   
    std::vector<bool>                   m_activeNodes;   
    int m_maxFields;                                
    int m_nodes;                                    
    bool m_polyEnabled;                             
    std::vector<int> m_polyActiveBasis;             
    LinearPolynomial m_polynomial;                  
public:
    RBFKernel();
    RBFKernel(const RBFKernel & other);
 
    virtual ~RBFKernel() = default;
 
    void                    setFunction(RBFBasisFunction);
    void                    setFunction(double (&funct)(double ));
 
    RBFBasisFunction        getFunctionType();
    int                     getDataCount();
    int                     getActiveCount();
    std::vector<int>        getActiveSet();
 
    void                    enablePolynomial(bool enable = true);
    int                     getPolynomialDimension();
    int                     getPolynomialWeightsCount();
 
    bool                    isActive(int );
 
    bool                    activateNode(int );
    bool                    activateNode(const std::vector<int> &);
    void                    activateAllNodes();
    bool                    deactivateNode(int );
    bool                    deactivateNode(const std::vector<int> &);
    void                    deactivateAllNodes();
 
    void                    setSupportRadius(double);
    void                    setSupportRadius(const std::vector<double> &);
    double                  getSupportRadius();
    double                  getSupportRadius(int);
 
    void                    setMode(RBFMode mode);
    RBFMode                 getMode();
 
    void                    setDataToNode (int , const std::vector<double> &);
    void                    setDataToAllNodes(int , const std::vector<double> &);
 
    int                     addData();
    int                     addData(const std::vector<double> &);
    bool                    removeData(int);
    bool                    removeData(std::vector<int> &);
    void                    removeAllData();
 
    void                    fitDataToNodes();
    void                    fitDataToNodes(int);
 
    std::vector<double>     evalRBF(const std::array<double,3> &);
    std::vector<double>     evalRBF(int jnode);
    double                  evalBasis(double);
    double                  evalBasisPair(int i, int j);
 
    int                     solve();
    int                     greedy(double);
 
    const std::vector<std::vector<double>> & getWeights() const;
 
protected:
    double                  evalError();
    int                     addGreedyPoint();
    int                     solveLSQ();
    void                    swap(RBFKernel & x) noexcept;
 
private:
 
    virtual double calcDist(int i, int j) = 0;
    virtual double calcDist(const std::array<double,3> & point, int j)                  = 0;
    virtual std::vector<double> evalPolynomialBasis(int i)                              = 0;
    virtual std::vector<double> evalPolynomialBasis(const std::array<double,3> &point)  = 0;
    virtual void initializePolynomialActiveBasis()                                      = 0;
    virtual void initializePolynomial()                                                 = 0;
};
 
class RBF : public RBFKernel {
 
protected:
    std::vector<std::array<double,3>>   m_node;     
public:
    ~RBF();
    RBF(RBFBasisFunction = RBFBasisFunction::WENDLANDC2);
    RBF(const RBF & other);
    RBF & operator=(RBF other);
 
    int                     getTotalNodesCount();
 
    int                     addNode(const std::array<double,3> &);
    std::vector<int>        addNode(const std::vector<std::array<double,3>> &);
    bool                    removeNode(int);
    bool                    removeNode(std::vector<int> &);
    void                    removeAllNodes();
 
protected:
    void     swap(RBF & x) noexcept;
 
private:
    double calcDist(int i, int j) override;
    double calcDist(const std::array<double,3> & point, int j) override;
    void initializePolynomialActiveBasis() override;
    void initializePolynomial() override;
    std::vector<double> evalPolynomialBasis(int i) override;
    std::vector<double> evalPolynomialBasis(const std::array<double, 3> &point) override;
};
 
namespace rbf
{
    double                  wendlandc2(double);
    double                  linear(double);
    double                  gauss90(double);
    double                  gauss95(double);
    double                  gauss99(double);
    double                  c1c0(double);
    double                  c2c0(double);
    double                  c0c1(double);
    double                  c1c1(double);
    double                  c2c1(double);
    double                  c0c2(double);
    double                  c1c2(double);
    double                  c2c2(double);
    double                  cosinus(double);
    double                  thinplate(double);
} // namespace rbf
}
 
#endif