From 01d02185fdb1c11c099393d9030c4924aad2cf2e Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Tue, 28 Nov 2023 11:25:03 +0100 Subject: [PATCH 1/9] rbf: add linear polynomial term in interpolation and integration test rbf 00004 --- CMakeLists.txt | 2 +- src/RBF/rbf.cpp | 362 ++++++++++++++++-- src/RBF/rbf.hpp | 44 ++- test/integration_tests/RBF/CMakeLists.txt | 1 + test/integration_tests/RBF/test_RBF_00004.cpp | 141 +++++++ 5 files changed, 514 insertions(+), 36 deletions(-) create mode 100644 test/integration_tests/RBF/test_RBF_00004.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 90461b64a4..69a1003f88 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -580,7 +580,7 @@ set(SURFUNSTRUCTURED_DEPS "common;patchkernel;lineunstructured;IO") set(VOLCARTESIAN_DEPS "common;patchkernel") set(VOLOCTREE_DEPS "common;PABLO;patchkernel") set(VOLUNSTRUCTURED_DEPS "common;patchkernel") -set(RBF_DEPS "operators;IO") +set(RBF_DEPS "operators;IO;discretization") set(DISCRETIZATION_DEPS "common;containers;LA;patchkernel") set(LEVELSET_DEPS "common;communications;SA;CG;surfunstructured;voloctree;volcartesian;volunstructured;IO") set(POD_DEPS "IO;common;containers;voloctree") diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index c09dd499c3..a6c4fb7ee2 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -48,6 +48,8 @@ namespace bitpit { * Some internal methods of the class can change their behaviour according * to the class mode selected. Please check documentation of each single method to * appreciate the differences. Default mode of the class is the "INTERP" RBFMode:: + * In case of interpolation the contribution of a linear polynomial is added + * to regularize the interpolation problem. * The actual abstract class does not implement the type of node you want to use; it * can be a 3D point or an unstructured mesh. Anyway what you have to do in deriving your own * RBF class is to specify a type of node you want to use and implement an evaluator @@ -73,6 +75,9 @@ RBFKernel::RBFKernel() m_activeNodes.clear(); setFunction( RBFBasisFunction::WENDLANDC2); + + m_polyActiveBasis.clear(); + m_polynomial.clear(); } /*! @@ -83,7 +88,8 @@ RBFKernel::RBFKernel(const RBFKernel & other) m_supportRadii(other.m_supportRadii), m_typef(other.m_typef), m_fPtr(other.m_fPtr), m_error(other.m_error), m_value(other.m_value), m_weight(other.m_weight), m_activeNodes(other.m_activeNodes), - m_maxFields(other.m_maxFields), m_nodes(other.m_nodes) + m_maxFields(other.m_maxFields), m_nodes(other.m_nodes), + m_polyActiveBasis(other.m_polyActiveBasis), m_polynomial(other.m_polynomial) { } @@ -105,6 +111,8 @@ void RBFKernel::swap(RBFKernel &other) noexcept std::swap(m_activeNodes, other.m_activeNodes); std::swap(m_maxFields, other.m_maxFields); std::swap(m_nodes, other.m_nodes); + std::swap(m_polyActiveBasis, other.m_polyActiveBasis); + std::swap(m_polynomial, other.m_polynomial); } /*! @@ -612,6 +620,19 @@ std::vector RBFKernel::evalRBF( const std::array &point) } } + // If INTERP mode add the polynomial contribution + if (m_mode == RBFMode::INTERP) { + for (j = 0; j < m_fields; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + std::vector basis = evalPolynomialBasis(point); + int z = 0; + for (int iactive : m_polyActiveBasis) { + values[j] += polynomialCoefficients[iactive] * basis[z]; + z++; + } + } + } + return values; } @@ -644,6 +665,19 @@ std::vector RBFKernel::evalRBF(int jnode) } } + // If INTERP mode add the polynomial contribution + if (m_mode == RBFMode::INTERP) { + for (j = 0; j < m_fields; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + std::vector basis = evalPolynomialBasis(jnode); + int z = 0; + for (int iactive : m_polyActiveBasis) { + values[j] += polynomialCoefficients[iactive] * basis[z]; + z++; + } + } + } + return values; } @@ -660,9 +694,18 @@ int RBFKernel::solve() return -1; } + // Initialize polynomial + initializePolynomial(); + + // Check which parameter of the polynomial has to be activated + // in order to avoid undetermined system + initializePolynomialActiveBasis(); + double dist; - int nS = getActiveCount(); + int nActive = getActiveCount(); + int nPoly = m_polyActiveBasis.size(); + int nS = nActive + nPoly; int nrhs = getDataCount(); int lda = nS; @@ -676,24 +719,39 @@ int RBFKernel::solve() std::vector A(lda * nS); std::vector b(ldb * nrhs); - int k = 0; for (int j = 0; j < nrhs; ++j) { - for( const auto & i : activeSet ) { - b[k] = m_value[j][i]; - ++k; + for (const auto &i : activeSet) { + int k = j * ldb + i; + b[k] = m_value[j][i]; } } - k=0; - for( const auto &i : activeSet ) { - for( const auto &j : activeSet ){ - dist = calcDist(j, i) / getSupportRadius(i); //order by column! - A[k] = evalBasis( dist ); + int k = 0; + for (const auto &i : activeSet) { + for (const auto &j : activeSet) { + dist = calcDist(j, i) / getSupportRadius(i); //order by column! + int row = k % nActive; + int col = k / nActive; + A[(col * nS) + row] = evalBasis(dist); k++; } } - info = LAPACKE_dgesv( LAPACK_COL_MAJOR, nS, nrhs, A.data(), lda, ipiv.data(), b.data(), ldb ); + // Filling terms given by the polynomial block. + // Symmetric terms set in the loop. + if (nPoly != 0) { + for (int i = 0; i < nActive; ++i) { + std::vector polynomialTerms = evalPolynomialBasis(activeSet[i]); + int j = 0; + for (const double &val : polynomialTerms) { + A[(j + nActive) * nS + i] = val; + A[i * nS + (j + nActive)] = val; + j++; + } + } + } + + info = LAPACKE_dgesv(LAPACK_COL_MAJOR, nS, nrhs, A.data(), lda, ipiv.data(), b.data(), ldb); if( info > 0 ) { printf( "The diagonal element of the triangular factor of the linear system matrix \n" ); @@ -704,16 +762,23 @@ int RBFKernel::solve() m_weight.resize(nrhs); - k=0; for (int j = 0; j < nrhs; ++j) { - m_weight[j].resize(m_nodes,0); - - for( const auto &i : activeSet ) { - m_weight[j][i] = b[k]; + m_weight[j].resize(m_nodes, 0); + int k = 0; + for (const auto &i : activeSet) { + m_weight[j][i] = b[j * ldb + k]; ++k; } } + for (int j = 0; j < nrhs; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + int i = 0; + for (int iactive : m_polyActiveBasis) { + polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; + i++; + } + } return 0; } @@ -744,7 +809,6 @@ int RBFKernel::greedy( double tolerance) } m_error[i] = norm2(local); - } std::ios::fmtflags streamFlags(log::cout().flags()); @@ -766,7 +830,6 @@ int RBFKernel::greedy( double tolerance) } else { errorFlag = 1; break; - } } @@ -869,7 +932,6 @@ int RBFKernel::addGreedyPoint( ) } ++i; - } return index; @@ -928,9 +990,18 @@ int RBFKernel::solveLSQ() return -1; } + // Initialize polynomial + initializePolynomial(); + + // Check which parameter of the polynomial has to be activated + // in order to avoid undetermined system + initializePolynomialActiveBasis(); + double dist; - int nR = getActiveCount(); + int nActive = getActiveCount(); + int nPoly = getPolynomialWeightsCount(); + int nS = nActive + nPoly; int nP = m_nodes; int nrhs = getDataCount(); @@ -940,34 +1011,62 @@ int RBFKernel::solveLSQ() double rcond = -1.0; m = nP; - n = nR; + n = nS; - lda = m; - ldb = std::max(n,m); + lda = m + nPoly; + ldb = m + nPoly; std::vector A(lda * n); std::vector b(ldb * nrhs, 0.0); std::vector s(m); for (int j = 0; j < nrhs; ++j) { - for (int i = 0; i < nP; ++i) { + for (int i = 0; i < m; ++i) { int k = j * ldb + i; - b[k] = m_value[j][i]; + b[k] = m_value[j][i]; } } int k = 0; - for( const auto &i : activeSet ) { + for (const auto &i : activeSet) { for (int j = 0; j < nP; ++j) { - dist = calcDist(j, i) / getSupportRadius(i); //in column format - A[k] = evalBasis( dist ); + dist = calcDist(j, i) / getSupportRadius(i); //order by column! + int row = k % nP; + int col = k / nP; + A[(col * ldb) + row] = evalBasis(dist); k++; } } - info = LAPACKE_dgelsd( LAPACK_COL_MAJOR, nP, nR, nrhs, A.data(), lda, b.data(), ldb, s.data(), rcond, &rank ); + // Filling column terms given by the polynomial block. + if (nPoly != 0) { + // Loop on columns + for (int i = 0; i < nActive; ++i) { + std::vector polynomialTerms = evalPolynomialBasis(activeSet[i]); + int j = 0; + for (const double &val : polynomialTerms) { + A[lda * i + (nP + j)] = val; + j++; + } + } + } + + // Filling row terms given by the polynomial block. + if (nPoly != 0) { + // Loop on rows + for (int i = 0; i < nP; ++i) { + std::vector polynomialTerms = evalPolynomialBasis(i); + int j = 0; + for (const double &val : polynomialTerms) { + A[(lda * (nActive + j) + i)] = val; + j++; + } + } + } - if( info > 0 ) { + info = LAPACKE_dgelsd(LAPACK_COL_MAJOR, lda, nS, nrhs, A.data(), lda, b.data(), ldb, s.data(), rcond, &rank); + + if (info > 0) { return 1; } @@ -984,9 +1083,110 @@ int RBFKernel::solveLSQ() } } + for (int j = 0; j < nrhs; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + int i = 0; + for (int iactive : m_polyActiveBasis) { + polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; + i++; + } + } + return 0; } +/*! + * \return The dimension of active polynomial terms + */ +int RBFKernel::getPolynomialWeightsCount() +{ + int nTerms = m_polyActiveBasis.size(); + return nTerms; +} + +/*! + * Default constructor of LinearPolynomial. + * Dimension and number of fields set to zero. + */ +RBFKernel::LinearPolynomial::LinearPolynomial() +{ + setDefault(); +} + +/*! + * Set to default values the LinearPolynomial class. + * Dimension and number of fields set to zero. + */ +void RBFKernel::LinearPolynomial::setDefault() +{ + m_dim = 0; + m_fields = 0; +} + +/*! + * Clear LinearPolynomial object. + * Clear structures and set default values restored. + */ +void RBFKernel::LinearPolynomial::clear() +{ + ReconstructionPolynomial::clear(true); + m_dim = 0; + m_fields = 0; +} + + +/*! + * Initialization of polynomial object. + * Dimension and number of fields has to be already set by the user. + * The polynomial is defined using the global coordinates and centered in + * the origin of the global reference system. + * Polynomial coefficients of full polynomial are initialized to zero values. + */ +void RBFKernel::LinearPolynomial::initialize() +{ + std::array m_origin; + m_origin.fill(0.); + ReconstructionPolynomial::initialize(1, m_dim, m_origin, m_fields, true); + double *coeffs = getCoefficients(); + for (auto i = 0; i < getCoefficientCount(); i++) { + coeffs[i] = 0.; + } +} + +/*! + * Set the space dimension (number of variables) of the multivariate polynomial. + * \param[in] dim Space variable dimension + */ +void RBFKernel::LinearPolynomial::setDimension(int dim) +{ + m_dim = (dim > 0) ? (dim < 4 ? dim : 3) : 0; +} + +/*! + * Set the number of fields (equations) retained by the polynomial class. + * \param[in] fields Number of fields + */ +void RBFKernel::LinearPolynomial::setDataCount(int fields) +{ + m_fields = std::max(0, fields); +} + +/*! + * Evaluate the polynomial basis in a point of the space variables. + * \param[in] x Target point where the basis is evaluated. + * \param[out] basis Basis values + */ +void RBFKernel::LinearPolynomial::evalBasis(const double *x, double *basis) +{ + // Set 0-th degree coefficients + basis[0] = 1.; + + // Set 1-st degree coefficients + for (int i = 0; i < m_dim; ++i) { + basis[i + 1] = x[i]; + } +} + // RBF derived class implementation /*! @@ -1181,7 +1381,105 @@ double RBF::calcDist(int i, int j) */ double RBF::calcDist(const std::array& point, int j) { - return norm2(point-m_node[j]); + return norm2(point - m_node[j]); +} + +/*! + * Initialize activation of monomials terms of linear polynomial part + * If all the nodes are aligned on a plane normal to a cartesian coordinate + * the related coordinate if disable in the polynomial function (no dependency + * can be found). + */ +void RBF::initializePolynomialActiveBasis() +{ + std::vector activeSet(RBFKernel::getActiveSet()); + + // Initialize active terms for an only constant linear polynomial + m_polyActiveBasis.clear(); + m_polyActiveBasis.insert(0); + + // Initialize reference coordinates with the first node + std::array coord(m_node[0]); + + // Check if at least one node has one of the coordinates + // different from the reference ones and enable the related + // polynomial term + std::set coordinatesToCheck({0, 1, 2}); + double tol = 1.0e-12; + for (const auto &i : activeSet) { + std::array point = m_node[i]; + for (auto it = coordinatesToCheck.begin(); it != coordinatesToCheck.end();) { + if (!utils::DoubleFloatingEqual()(coord[*it], point[*it], tol, tol)) { + m_polyActiveBasis.insert(*it+1); + it = coordinatesToCheck.erase(it); + } else { + ++it; + } + } + } +} + +/*! + * Initialize polynomial terms. + */ +void RBF::initializePolynomial() +{ + m_polynomial.clear(); + m_polynomial.setDimension(3); + m_polynomial.setDataCount(getDataCount()); + m_polynomial.initialize(); +} + +/*! + * Compute monomials basis values of linear polynomial part on a RBF node + * @param[in] i i-th RBF node in the list + * \return Values of polynomial basis + */ +std::vector RBF::evalPolynomialBasis(int i) +{ + int nPoly = m_polyActiveBasis.size(); + std::vector result(nPoly); + + if (nPoly < 1) + return result; + + const std::array &point = m_node[i]; + + std::vector completeBasis(m_polynomial.getCoefficientCount()); + m_polynomial.evalBasis(point.data(), completeBasis.data()); + + int k = 0; + for (int j : m_polyActiveBasis) { + result[k] = completeBasis[j]; + k++; + } + + return result; +} + +/*! + * Compute monomials basis values of linear polynomial part on a target point + * @param[in] point point on where to evaluate the basis + * \return Values of polynomial basis + */ +std::vector RBF::evalPolynomialBasis(const std::array &point) +{ + int nPoly = m_polyActiveBasis.size(); + std::vector result(nPoly); + + if (nPoly < 1) + return result; + + std::vector completeBasis(m_polynomial.getCoefficientCount()); + m_polynomial.evalBasis(point.data(), completeBasis.data()); + + int k = 0; + for (int j : m_polyActiveBasis) { + result[k] = completeBasis[j]; + k++; + } + + return result; } // RBF NAMESPACE UTILITIES diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 712cef4913..7228089426 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -27,6 +27,8 @@ #include #include +#include +#include "reconstruction.hpp" namespace bitpit{ @@ -64,6 +66,28 @@ enum class RBFMode { }; class RBFKernel{ + /*! + * Linear Polynomial class. A linear multivariate polynomial in function of the + * variables defining the RBF nodes is used to regularize the interpolation + * problem through RBFs. Used only when RBFMode::PARAM is enabled. + * The LinearPolynomial class derives from ReconstructionPolynomial class + * that is intrisecally defined to work in a three dimensional space. + */ + class LinearPolynomial : public ReconstructionPolynomial + { +private: + int m_dim; /**< Polynomial variables space dimension */ + int m_fields; /**< Number of equations of the polynomial object */ + + public: + LinearPolynomial(); + void setDefault(); + void clear(); + void setDimension(int dim); + void setDataCount(int fields); + void evalBasis(const double *x, double *basis); + void initialize(); + }; private: int m_fields; /** m_error; /**> m_value; /**< displ value to be interpolated on RBF nodes */ std::vector> m_weight; /**< weight of your RBF interpolation */ std::vector m_activeNodes; /** the i-th node is used/not used during RBF evaluation).*/ int m_maxFields; /**< fix the maximum number of fields that can be added to your class*/ int m_nodes; /** m_polyActiveBasis; /**< Active terms of linear polynomial, -1 is constant, 0,1,2 the system coordinates */ + + LinearPolynomial m_polynomial; public: RBFKernel(); @@ -95,6 +122,10 @@ class RBFKernel{ int getActiveCount(); std::vector getActiveSet(); + void setPolynomialDimension(int dim); + int getPolynomialDimension(); + int getPolynomialWeightsCount(); + bool isActive(int ); bool activateNode(int ); @@ -143,8 +174,11 @@ class RBFKernel{ private: virtual double calcDist(int i, int j) = 0; - virtual double calcDist(const std::array & point, int j) = 0; - + virtual double calcDist(const std::array & point, int j) = 0; + virtual std::vector evalPolynomialBasis(int i) = 0; + virtual std::vector evalPolynomialBasis(const std::array &point) = 0; + virtual void initializePolynomialActiveBasis() = 0; + virtual void initializePolynomial() = 0; }; class RBF : public RBFKernel { @@ -172,6 +206,10 @@ class RBF : public RBFKernel { private: double calcDist(int i, int j) override; double calcDist(const std::array & point, int j) override; + void initializePolynomialActiveBasis() override; + void initializePolynomial() override; + std::vector evalPolynomialBasis(int i) override; + std::vector evalPolynomialBasis(const std::array &point) override; }; /*! diff --git a/test/integration_tests/RBF/CMakeLists.txt b/test/integration_tests/RBF/CMakeLists.txt index ecf456132c..d7e5911c4e 100644 --- a/test/integration_tests/RBF/CMakeLists.txt +++ b/test/integration_tests/RBF/CMakeLists.txt @@ -30,6 +30,7 @@ set(TESTS "") list(APPEND TESTS "test_RBF_00001") list(APPEND TESTS "test_RBF_00002") list(APPEND TESTS "test_RBF_00003") +list(APPEND TESTS "test_RBF_00004") # Test extra modules set(TEST_EXTRA_MODULES "") diff --git a/test/integration_tests/RBF/test_RBF_00004.cpp b/test/integration_tests/RBF/test_RBF_00004.cpp new file mode 100644 index 0000000000..2e3b8bf4bf --- /dev/null +++ b/test/integration_tests/RBF/test_RBF_00004.cpp @@ -0,0 +1,141 @@ +/*---------------------------------------------------------------------------*\ + * + * 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 . + * +\*---------------------------------------------------------------------------*/ + +#include +#include +#include +#include +#if BITPIT_ENABLE_MPI==1 +#include +#endif + +#include "bitpit_RBF.hpp" + +using namespace bitpit; + +//Description. Directly Interpolate a scalar (z-displacement) field + +/*! +* Subtest 001 +* +* Testing field interpolation. +*/ +int subtest_001() +{ + // Variables + std::vector< std::array > controlNodes; + std::vector< std::array > points; + std::vector< double > zDispl ; + + // fill control points + int nCP = 4; + controlNodes.resize(nCP); + controlNodes[0] = std::array({0., 0., 0.}); + controlNodes[1] = std::array({0., 2., 0.}); + controlNodes[2] = std::array({2., 2., 0.}); + controlNodes[3] = std::array({2., 0., 0.}); + + // fill displacements of control points + zDispl.resize(nCP); + zDispl[0] = 1.; + zDispl[1] = 1.; + zDispl[2] = 2.; + zDispl[3] = 1.; + + + //RBF used to directly interpolate zDispl field + bitpit::RBF paraMorph; + RBFBasisFunction funct = RBFBasisFunction::WENDLANDC2; + paraMorph.setFunction(funct); + paraMorph.setSupportRadius(0.5) ; + + std::vector nIndex = paraMorph.addNode(controlNodes); + int nData = paraMorph.addData(zDispl); + int err = paraMorph.solve() ; + if(err > 0 ) return 1; + + // create moving points + points.resize(nCP+1); + points[0] = std::array({0., 0., 0.}); + points[1] = std::array({0., 2., 0.}); + points[2] = std::array({2., 2., 0.}); + points[3] = std::array({2., 0., 0.}); + points[4] = std::array({1., 1., 0.}); + + std::vector disp; + for( auto & point : points){ + std::cout << point << " -----> "; + disp = paraMorph.evalRBF(point); + point[2] += disp[0] ; + std::cout << point << std::endl; + }; + + double val = points[4][2]; + return (!bitpit::utils::DoubleFloatingEqual()(val, 1.25)); +} + +// ========================================================================== // +// MAIN // +// ========================================================================== // +int main(int argc, char *argv[]) +{ + // ====================================================================== // + // INITIALIZE MPI // + // ====================================================================== // +#if BITPIT_ENABLE_MPI==1 + MPI_Init(&argc,&argv); +#else + BITPIT_UNUSED(argc); + BITPIT_UNUSED(argv); +#endif + + // ====================================================================== // + // VARIABLES DECLARATION // + // ====================================================================== // + + // Local variabels + int status = 0; + + // ====================================================================== // + // RUN SUB-TESTS // + // ====================================================================== // + try { + status = subtest_001(); + if (status != 0) { + return (10 + status); + } + } catch (const std::exception &exception) { + bitpit::log::cout() << exception.what(); + exit(1); + } + + // ====================================================================== // + // FINALIZE MPI // + // ====================================================================== // +#if BITPIT_ENABLE_MPI==1 + MPI_Finalize(); +#endif + + return status; +} From 3fb90ae3e0b6bc12bd1f325f404bd1f35c6d817f Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Mon, 4 Dec 2023 16:21:40 +0100 Subject: [PATCH 2/9] rbf: add enable/disable polynomial and minor fixes --- src/RBF/rbf.cpp | 102 ++++++++++++++++++++++++++++-------------------- src/RBF/rbf.hpp | 6 +-- 2 files changed, 63 insertions(+), 45 deletions(-) diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index a6c4fb7ee2..319d601a0a 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -48,8 +48,8 @@ namespace bitpit { * Some internal methods of the class can change their behaviour according * to the class mode selected. Please check documentation of each single method to * appreciate the differences. Default mode of the class is the "INTERP" RBFMode:: - * In case of interpolation the contribution of a linear polynomial is added - * to regularize the interpolation problem. + * In case of interpolation the contribution of a linear polynomial can be added + * to regularize the interpolation problem by enabling the polynomial usage. * The actual abstract class does not implement the type of node you want to use; it * can be a 3D point or an unstructured mesh. Anyway what you have to do in deriving your own * RBF class is to specify a type of node you want to use and implement an evaluator @@ -76,6 +76,7 @@ RBFKernel::RBFKernel() setFunction( RBFBasisFunction::WENDLANDC2); + m_polyEnabled = false; m_polyActiveBasis.clear(); m_polynomial.clear(); } @@ -89,7 +90,8 @@ RBFKernel::RBFKernel(const RBFKernel & other) m_fPtr(other.m_fPtr), m_error(other.m_error), m_value(other.m_value), m_weight(other.m_weight), m_activeNodes(other.m_activeNodes), m_maxFields(other.m_maxFields), m_nodes(other.m_nodes), - m_polyActiveBasis(other.m_polyActiveBasis), m_polynomial(other.m_polynomial) + m_polyEnabled(other.m_polyEnabled), m_polyActiveBasis(other.m_polyActiveBasis), + m_polynomial(other.m_polynomial) { } @@ -111,8 +113,9 @@ void RBFKernel::swap(RBFKernel &other) noexcept std::swap(m_activeNodes, other.m_activeNodes); std::swap(m_maxFields, other.m_maxFields); std::swap(m_nodes, other.m_nodes); - std::swap(m_polyActiveBasis, other.m_polyActiveBasis); - std::swap(m_polynomial, other.m_polynomial); + std::swap(m_polyEnabled, other.m_polyEnabled); + std::swap(m_polyActiveBasis, other.m_polyActiveBasis); + std::swap(m_polynomial, other.m_polynomial); } /*! @@ -621,7 +624,7 @@ std::vector RBFKernel::evalRBF( const std::array &point) } // If INTERP mode add the polynomial contribution - if (m_mode == RBFMode::INTERP) { + if (m_mode == RBFMode::INTERP && m_polyEnabled) { for (j = 0; j < m_fields; ++j) { double *polynomialCoefficients = m_polynomial.getCoefficients(j); std::vector basis = evalPolynomialBasis(point); @@ -666,7 +669,7 @@ std::vector RBFKernel::evalRBF(int jnode) } // If INTERP mode add the polynomial contribution - if (m_mode == RBFMode::INTERP) { + if (m_mode == RBFMode::INTERP && m_polyEnabled) { for (j = 0; j < m_fields; ++j) { double *polynomialCoefficients = m_polynomial.getCoefficients(j); std::vector basis = evalPolynomialBasis(jnode); @@ -694,17 +697,19 @@ int RBFKernel::solve() return -1; } - // Initialize polynomial - initializePolynomial(); + if (m_polyEnabled) { + // Initialize polynomial + initializePolynomial(); - // Check which parameter of the polynomial has to be activated - // in order to avoid undetermined system - initializePolynomialActiveBasis(); + // Check which parameter of the polynomial has to be activated + // in order to avoid undetermined system + initializePolynomialActiveBasis(); + } double dist; int nActive = getActiveCount(); - int nPoly = m_polyActiveBasis.size(); + int nPoly = m_polyEnabled ? m_polyActiveBasis.size() : 0; int nS = nActive + nPoly; int nrhs = getDataCount(); @@ -739,7 +744,7 @@ int RBFKernel::solve() // Filling terms given by the polynomial block. // Symmetric terms set in the loop. - if (nPoly != 0) { + if (m_polyEnabled) { for (int i = 0; i < nActive; ++i) { std::vector polynomialTerms = evalPolynomialBasis(activeSet[i]); int j = 0; @@ -771,12 +776,14 @@ int RBFKernel::solve() } } - for (int j = 0; j < nrhs; ++j) { - double *polynomialCoefficients = m_polynomial.getCoefficients(j); - int i = 0; - for (int iactive : m_polyActiveBasis) { - polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; - i++; + if (m_polyEnabled) { + for (int j = 0; j < nrhs; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + int i = 0; + for (int iactive : m_polyActiveBasis) { + polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; + i++; + } } } return 0; @@ -990,17 +997,19 @@ int RBFKernel::solveLSQ() return -1; } - // Initialize polynomial - initializePolynomial(); + if (m_polyEnabled) { + // Initialize polynomial + initializePolynomial(); - // Check which parameter of the polynomial has to be activated - // in order to avoid undetermined system - initializePolynomialActiveBasis(); + // Check which parameter of the polynomial has to be activated + // in order to avoid undetermined system + initializePolynomialActiveBasis(); + } double dist; int nActive = getActiveCount(); - int nPoly = getPolynomialWeightsCount(); + int nPoly = m_polyEnabled ? getPolynomialWeightsCount() : 0; int nS = nActive + nPoly; int nP = m_nodes; int nrhs = getDataCount(); @@ -1030,20 +1039,20 @@ int RBFKernel::solveLSQ() int k = 0; for (const auto &i : activeSet) { for (int j = 0; j < nP; ++j) { - dist = calcDist(j, i) / getSupportRadius(i); //order by column! - int row = k % nP; - int col = k / nP; + dist = calcDist(j, i) / getSupportRadius(i); //order by column! + int row = k % nP; + int col = k / nP; A[(col * ldb) + row] = evalBasis(dist); k++; } } // Filling column terms given by the polynomial block. - if (nPoly != 0) { + if (m_polyEnabled) { // Loop on columns for (int i = 0; i < nActive; ++i) { std::vector polynomialTerms = evalPolynomialBasis(activeSet[i]); - int j = 0; + int j = 0; for (const double &val : polynomialTerms) { A[lda * i + (nP + j)] = val; j++; @@ -1052,11 +1061,11 @@ int RBFKernel::solveLSQ() } // Filling row terms given by the polynomial block. - if (nPoly != 0) { + if (m_polyEnabled) { // Loop on rows for (int i = 0; i < nP; ++i) { std::vector polynomialTerms = evalPolynomialBasis(i); - int j = 0; + int j = 0; for (const double &val : polynomialTerms) { A[(lda * (nActive + j) + i)] = val; j++; @@ -1083,16 +1092,26 @@ int RBFKernel::solveLSQ() } } - for (int j = 0; j < nrhs; ++j) { - double *polynomialCoefficients = m_polynomial.getCoefficients(j); - int i = 0; - for (int iactive : m_polyActiveBasis) { - polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; - i++; + if (m_polyEnabled) { + for (int j = 0; j < nrhs; ++j) { + double *polynomialCoefficients = m_polynomial.getCoefficients(j); + int i = 0; + for (int iactive : m_polyActiveBasis) { + polynomialCoefficients[iactive] = b[j * ldb + nActive + i]; + i++; + } } } + return 0; + } - return 0; +/*! + * Enable/disable the use of polynomial term during interpolation. + * \param[in] enable true/false to enable/disable polynomial usage term + */ +void RBFKernel::enablePolynomial(bool enable) +{ + m_polyEnabled = enable; } /*! @@ -1405,11 +1424,10 @@ void RBF::initializePolynomialActiveBasis() // different from the reference ones and enable the related // polynomial term std::set coordinatesToCheck({0, 1, 2}); - double tol = 1.0e-12; for (const auto &i : activeSet) { std::array point = m_node[i]; for (auto it = coordinatesToCheck.begin(); it != coordinatesToCheck.end();) { - if (!utils::DoubleFloatingEqual()(coord[*it], point[*it], tol, tol)) { + if (!utils::DoubleFloatingEqual()(coord[*it], point[*it])) { m_polyActiveBasis.insert(*it+1); it = coordinatesToCheck.erase(it); } else { diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 7228089426..594f69c205 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -104,10 +104,10 @@ class RBFKernel{ std::vector m_activeNodes; /** the i-th node is used/not used during RBF evaluation).*/ int m_maxFields; /**< fix the maximum number of fields that can be added to your class*/ int m_nodes; /** m_polyActiveBasis; /**< Active terms of linear polynomial, -1 is constant, 0,1,2 the system coordinates */ - LinearPolynomial m_polynomial; - public: RBFKernel(); RBFKernel(const RBFKernel & other); @@ -122,7 +122,7 @@ class RBFKernel{ int getActiveCount(); std::vector getActiveSet(); - void setPolynomialDimension(int dim); + void enablePolynomial(bool enable = true); int getPolynomialDimension(); int getPolynomialWeightsCount(); From 15106c4d5306c0787db9c5682d28d8f2a498beba Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Mon, 4 Dec 2023 16:26:00 +0100 Subject: [PATCH 3/9] rbf: add thin plate function and modify test 00004 --- src/RBF/rbf.cpp | 19 +++++++ src/RBF/rbf.hpp | 5 +- test/integration_tests/RBF/test_RBF_00004.cpp | 54 ++++++++++++++++--- 3 files changed, 69 insertions(+), 9 deletions(-) diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index 319d601a0a..0f96f6e040 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -196,6 +196,11 @@ void RBFKernel::setFunction( RBFBasisFunction bfunc ) m_typef = RBFBasisFunction::COSINUS; break; + case (RBFBasisFunction::THINPLATE): + setFunction(rbf::thinplate); + m_typef = RBFBasisFunction::THINPLATE; + break; + default: setFunction( rbf::wendlandc2); m_typef = RBFBasisFunction::WENDLANDC2; @@ -1701,4 +1706,18 @@ double rbf::c2c2( double dist ) } } +/*! + * Non compact thin plate spline function. + * @param[in] dist distance normalized with respect to support radius + * @return rbf value + */ +double rbf::thinplate(double dist) +{ + if (utils::DoubleFloatingGreater()(dist, 0.)) { + return dist * dist * std::log(dist); + } else { + return 0.; + } } + +} // namespace bitpit diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 594f69c205..980dd5f12c 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -53,6 +53,7 @@ enum class RBFBasisFunction { C1C2 = 12, /**< Compact biquadratic funct, C1 on r=0, C2 on r=1, 0 outside */ C2C2 = 13, /**< Compact poly (degree 5) funct, C2 on r=0, C2 on r=1, 0 outside */ COSINUS = 14, /**< Compact cosinusoidal funct, value of 1 on r=0, 0 outside */ + THINPLATE = 15, /**< Non compact thin plate funct */ }; /*! @@ -232,8 +233,8 @@ namespace rbf double c1c2(double); double c2c2(double); double cosinus(double); -} - + double thinplate(double); +} // namespace rbf } #endif diff --git a/test/integration_tests/RBF/test_RBF_00004.cpp b/test/integration_tests/RBF/test_RBF_00004.cpp index 2e3b8bf4bf..d8dce798eb 100644 --- a/test/integration_tests/RBF/test_RBF_00004.cpp +++ b/test/integration_tests/RBF/test_RBF_00004.cpp @@ -64,14 +64,22 @@ int subtest_001() zDispl[3] = 1.; - //RBF used to directly interpolate zDispl field - bitpit::RBF paraMorph; - RBFBasisFunction funct = RBFBasisFunction::WENDLANDC2; + // RBF used to directly interpolate zDispl field + std::vector supportRadii; + supportRadii.resize(nCP); + supportRadii[0] = 1.; + supportRadii[1] = 0.75; + supportRadii[2] = 1.25; + supportRadii[3] = 0.75; + + bitpit::RBF paraMorph; + RBFBasisFunction funct = RBFBasisFunction::WENDLANDC2; paraMorph.setFunction(funct); - paraMorph.setSupportRadius(0.5) ; - + paraMorph.setSupportRadius(supportRadii); + paraMorph.enablePolynomial(); + std::vector nIndex = paraMorph.addNode(controlNodes); - int nData = paraMorph.addData(zDispl); + paraMorph.addData(zDispl); int err = paraMorph.solve() ; if(err > 0 ) return 1; @@ -92,7 +100,39 @@ int subtest_001() }; double val = points[4][2]; - return (!bitpit::utils::DoubleFloatingEqual()(val, 1.25)); + bool check = !bitpit::utils::DoubleFloatingEqual()(val, 1.25); + + + // Use a different function + paraMorph.removeAllData(); + funct = RBFBasisFunction::THINPLATE; + paraMorph.setFunction(funct); + paraMorph.setSupportRadius(1.); + paraMorph.enablePolynomial(); + + paraMorph.addData(zDispl); + err = paraMorph.solve(); + if (err > 0) + return 1; + + // reset moving points + points[0] = std::array({0., 0., 0.}); + points[1] = std::array({0., 2., 0.}); + points[2] = std::array({2., 2., 0.}); + points[3] = std::array({2., 0., 0.}); + points[4] = std::array({1., 1., 0.}); + + // move points + for (auto &point : points) { + std::cout << point << " -----> "; + disp = paraMorph.evalRBF(point); + point[2] += disp[0]; + std::cout << point << std::endl; + }; + + val = points[4][2]; + check != !bitpit::utils::DoubleFloatingEqual()(val, 1.25); + return check; } // ========================================================================== // From 79a8e9d4066627ff29d45ac6d78c1bca1692ddc5 Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Mon, 4 Dec 2023 16:31:16 +0100 Subject: [PATCH 4/9] rbf: fix documentation of LinearPolynomial --- src/RBF/rbf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 980dd5f12c..79af3c07ff 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -70,7 +70,7 @@ class RBFKernel{ /*! * Linear Polynomial class. A linear multivariate polynomial in function of the * variables defining the RBF nodes is used to regularize the interpolation - * problem through RBFs. Used only when RBFMode::PARAM is enabled. + * problem through RBFs. It can be used only when RBFMode::INTERP is enabled. * The LinearPolynomial class derives from ReconstructionPolynomial class * that is intrisecally defined to work in a three dimensional space. */ From 8e29df70cdc12cf55aa3015cb14dad28792f2898 Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Mon, 4 Dec 2023 19:13:35 +0100 Subject: [PATCH 5/9] rbf: minor fixing --- src/RBF/rbf.cpp | 4 ++-- src/RBF/rbf.hpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index 0f96f6e040..f0734de8a2 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -1134,14 +1134,14 @@ int RBFKernel::getPolynomialWeightsCount() */ RBFKernel::LinearPolynomial::LinearPolynomial() { - setDefault(); + reset(); } /*! * Set to default values the LinearPolynomial class. * Dimension and number of fields set to zero. */ -void RBFKernel::LinearPolynomial::setDefault() +void RBFKernel::LinearPolynomial::reset() { m_dim = 0; m_fields = 0; diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 79af3c07ff..eb251cfe1c 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -25,10 +25,10 @@ #ifndef __BITPIT_RBF_HPP__ #define __BITPIT_RBF_HPP__ -#include +#include "bitpit_discretization.hpp" #include #include -#include "reconstruction.hpp" +#include namespace bitpit{ @@ -82,7 +82,7 @@ class RBFKernel{ public: LinearPolynomial(); - void setDefault(); + void reset(); void clear(); void setDimension(int dim); void setDataCount(int fields); From 181a308841c16213f39e4b349a86bc34e691ccbf Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Wed, 6 Dec 2023 16:52:44 +0100 Subject: [PATCH 6/9] rbf: remove reset method from linear polynomial class --- src/RBF/rbf.cpp | 9 --------- src/RBF/rbf.hpp | 1 - 2 files changed, 10 deletions(-) diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index f0734de8a2..191ff5faa1 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -1133,15 +1133,6 @@ int RBFKernel::getPolynomialWeightsCount() * Dimension and number of fields set to zero. */ RBFKernel::LinearPolynomial::LinearPolynomial() -{ - reset(); -} - -/*! - * Set to default values the LinearPolynomial class. - * Dimension and number of fields set to zero. - */ -void RBFKernel::LinearPolynomial::reset() { m_dim = 0; m_fields = 0; diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index eb251cfe1c..ccf94d3b2b 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -82,7 +82,6 @@ class RBFKernel{ public: LinearPolynomial(); - void reset(); void clear(); void setDimension(int dim); void setDataCount(int fields); From a6a235514914408c1866740b420d95891e57683e Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Wed, 6 Dec 2023 17:16:05 +0100 Subject: [PATCH 7/9] rbf: fixed documentation --- src/RBF/rbf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index ccf94d3b2b..8e04319ded 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -106,7 +106,7 @@ class RBFKernel{ int m_nodes; /** m_polyActiveBasis; /**< Active terms of linear polynomial, -1 is constant, 0,1,2 the system coordinates */ + std::set m_polyActiveBasis; /**< Active terms of linear polynomial, 0 is constant, i+1 the i-th system coordinate */ public: RBFKernel(); From 1ec113e567f093990e34bab4b8b929cf2db232e8 Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Wed, 6 Dec 2023 17:41:40 +0100 Subject: [PATCH 8/9] integration test: fixing integration test rbf 0004 --- test/integration_tests/RBF/test_RBF_00004.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/integration_tests/RBF/test_RBF_00004.cpp b/test/integration_tests/RBF/test_RBF_00004.cpp index d8dce798eb..2878460192 100644 --- a/test/integration_tests/RBF/test_RBF_00004.cpp +++ b/test/integration_tests/RBF/test_RBF_00004.cpp @@ -131,7 +131,7 @@ int subtest_001() }; val = points[4][2]; - check != !bitpit::utils::DoubleFloatingEqual()(val, 1.25); + check = check || !bitpit::utils::DoubleFloatingEqual()(val, 1.25); return check; } From aea7b4ab332df2e3e7cdb226f7ef86e59f8e2739 Mon Sep 17 00:00:00 2001 From: edoardolombardi Date: Thu, 7 Dec 2023 16:59:46 +0100 Subject: [PATCH 9/9] rbf: use vector for m_polyActiveBasis --- src/RBF/rbf.cpp | 17 +++++++---------- src/RBF/rbf.hpp | 2 +- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/RBF/rbf.cpp b/src/RBF/rbf.cpp index 191ff5faa1..076ce26dfb 100644 --- a/src/RBF/rbf.cpp +++ b/src/RBF/rbf.cpp @@ -1411,7 +1411,7 @@ void RBF::initializePolynomialActiveBasis() // Initialize active terms for an only constant linear polynomial m_polyActiveBasis.clear(); - m_polyActiveBasis.insert(0); + m_polyActiveBasis.push_back(0); // Initialize reference coordinates with the first node std::array coord(m_node[0]); @@ -1419,15 +1419,12 @@ void RBF::initializePolynomialActiveBasis() // Check if at least one node has one of the coordinates // different from the reference ones and enable the related // polynomial term - std::set coordinatesToCheck({0, 1, 2}); - for (const auto &i : activeSet) { - std::array point = m_node[i]; - for (auto it = coordinatesToCheck.begin(); it != coordinatesToCheck.end();) { - if (!utils::DoubleFloatingEqual()(coord[*it], point[*it])) { - m_polyActiveBasis.insert(*it+1); - it = coordinatesToCheck.erase(it); - } else { - ++it; + for (int d = 0; d < 2; ++d) { + for (const auto &i : activeSet) { + const std::array &point = m_node[i]; + if (!utils::DoubleFloatingEqual()(coord[d], point[d])) { + utils::addToOrderedVector(d + 1, m_polyActiveBasis); + break; } } } diff --git a/src/RBF/rbf.hpp b/src/RBF/rbf.hpp index 8e04319ded..cd856d48fa 100644 --- a/src/RBF/rbf.hpp +++ b/src/RBF/rbf.hpp @@ -106,7 +106,7 @@ class RBFKernel{ int m_nodes; /** m_polyActiveBasis; /**< Active terms of linear polynomial, 0 is constant, i+1 the i-th system coordinate */ + std::vector m_polyActiveBasis; /**< Active terms of linear polynomial, 0 is constant, i+1 the i-th system coordinate */ public: RBFKernel();