MRBF.cpp

/*---------------------------------------------------------------------------*\
 *
 *  mimmo
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of mimmo.
 *
 *  mimmo 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.
 *
 *  mimmo 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 mimmo. If not, see <http://www.gnu.org/licenses/>.
 *
 \ *---------------------------------------------------------------------------*/
 
#include "MRBF.hpp"
 
namespace mimmo{
 
    //TODO study how to manipulate supportRadius of RBF to define a local/global smoothing of RBF
 
MRBF::MRBF(MRBFSol mode){
    m_name = "mimmo.MRBF";
    setMode(mode);
    m_maxFields=-1;
    m_tol = 1.0E-06;
    m_bfilter = false;
    m_supportRadiusValue = -1.0;
    m_srIsReal = false;
    m_functype = -1;
    m_isCompact = false;
    m_rbfgeometry = nullptr;
    m_rbfdispl = nullptr;
    m_rbfScalarDispl = nullptr;
    m_rbfSupportRadii = nullptr;
    m_diagonalFactor = 1.0;
    m_areScalarResults = false;
};
 
MRBF::MRBF(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.MRBF";
    m_maxFields=-1;
    m_tol = 1.0E-06;
    m_bfilter = false;
    m_supportRadiusValue = -1.0;
    m_srIsReal = false;
    m_functype = -1;
    m_isCompact = false;
    m_rbfgeometry = nullptr;
    m_rbfdispl = nullptr;
    m_rbfScalarDispl = nullptr;
    m_rbfSupportRadii = nullptr;
    m_diagonalFactor = 1.0;
    m_areScalarResults = false;
 
    setMode(MRBFSol::NONE);
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
    if(input == "mimmo.MRBF"){
        std::string fallback_mode = "0";
        std::string input2 = rootXML.get("Mode", fallback_mode);
        input2 = bitpit::utils::string::trim(input2);
        int mode_int = std::stoi(input2);
        mode_int = std::min(2, std::max(0, mode_int));
        setMode(static_cast<MRBFSol>(mode_int));
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
MRBF::~MRBF(){};
 
MRBF::MRBF(const MRBF & other):BaseManipulation(other), bitpit::RBF(other){
    m_tol = other.m_tol;
    m_solver = other.m_solver;
    m_bfilter = other.m_bfilter;
    m_supportRadiusValue  = other.m_supportRadiusValue;
    m_srIsReal = other.m_srIsReal;
    m_supportRadii = other.m_supportRadii;
    m_functype = other.m_functype;
    m_isCompact = other.m_isCompact;
    if(m_bfilter)    m_filter = other.m_filter;
    m_rbfgeometry = other.m_rbfgeometry;
    m_rbfdispl = other.m_rbfdispl;
    m_rbfScalarDispl = other.m_rbfScalarDispl;
    m_rbfSupportRadii = other.m_rbfSupportRadii;
    m_diagonalFactor = other.m_diagonalFactor;
    m_areScalarResults = other.m_areScalarResults;
};
 
MRBF& MRBF::operator=(MRBF other){
    swap(other);
    return *this;
};
 
void MRBF::swap(MRBF & x) noexcept
{
    std::swap(m_tol, x.m_tol);
    std::swap(m_solver, x.m_solver);
    m_filter.swap(x.m_filter);
    std::swap(m_bfilter, x.m_bfilter);
    std::swap(m_supportRadiusValue , x.m_supportRadiusValue);
    std::swap(m_srIsReal, x.m_srIsReal);
    std::swap(m_supportRadii, x.m_supportRadii);
    std::swap(m_functype, x.m_functype);
    std::swap(m_isCompact, x.m_isCompact);
    m_displ.swap(x.m_displ);
    m_scalarDispl.swap(x.m_scalarDispl);
    std::swap(m_effectiveSR, x.m_effectiveSR);
    std::swap(m_rbfgeometry, x.m_rbfgeometry);
    std::swap(m_rbfdispl, x.m_rbfdispl);
    std::swap(m_rbfScalarDispl, x.m_rbfScalarDispl);
    std::swap(m_rbfSupportRadii, x.m_rbfSupportRadii);
    std::swap(m_diagonalFactor, x.m_diagonalFactor);
    std::swap(m_areScalarResults, x.m_areScalarResults);
 
    RBF::swap(x);
 
    BaseManipulation::swap(x);
}
 
void
MRBF::buildPorts(){
    bool built = true;
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, MRBF>(&m_geometry, M_GEOM, true));
    built = (built && createPortIn<dvecarr3E, MRBF>(this, &mimmo::MRBF::setNode, M_COORDS));
    built = (built && createPortIn<dvecarr3E, MRBF>(this, &mimmo::MRBF::setDisplacements, M_DISPLS));
    built = (built && createPortIn<std::vector<double>, MRBF>(this, &mimmo::MRBF::setScalarDisplacements, M_DATAFIELD));
    built = (built && createPortIn<std::vector<double>, MRBF>(this, &mimmo::MRBF::setVariableSupportRadii, M_DATAFIELD2));
    built = (built && createPortIn<dmpvector1D*, MRBF>(this, &mimmo::MRBF::setFilter, M_FILTER));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, MRBF>(this, &mimmo::MRBF::setNode, M_GEOM2));
    built = (built && createPortIn<dmpvecarr3E*, MRBF>(this, &mimmo::MRBF::setDisplacements, M_VECTORFIELD));
    built = (built && createPortIn<dmpvector1D*, MRBF>(this, &mimmo::MRBF::setScalarDisplacements, M_SCALARFIELD));
    built = (built && createPortIn<dmpvector1D*, MRBF>(this, &mimmo::MRBF::setVariableSupportRadii, M_SCALARFIELD2));
 
    built = (built && createPortOut<dmpvecarr3E*, MRBF>(this, &mimmo::MRBF::getDisplacements, M_GDISPLS));
    built = (built && createPortOut<dmpvector1D*, MRBF>(this, &mimmo::MRBF::getScalarDisplacements, M_SCALARFIELD));
    built = (built && createPortOut<MimmoSharedPointer<MimmoObject>, MRBF>(this, &BaseManipulation::getGeometry, M_GEOM));
    m_arePortsBuilt = built;
};
 
dvecarr3E*
MRBF::getNodes(){
    return(&m_node);
}
 
MRBFSol
MRBF::getMode(){
    return m_solver;
};
 
dmpvector1D*
MRBF::getFilter(){
    return &m_filter;
};
 
bool
MRBF::isSupportRadiusReal(){
    return m_srIsReal;
}
 
bool
MRBF::isVariableSupportRadiusSet(){
    return !m_supportRadii.empty();
}
 
dvector1D &
MRBF::getEffectivelyUsedSupportRadii(){
    return m_effectiveSR;
}
 
double
MRBF::getDiagonalFactor(){
    return m_diagonalFactor;
}
 
int
MRBF::getFunctionType(){
    if(m_functype < 0){
        return static_cast<int>(RBF::getFunctionType());
    }else{
        return m_functype;
    }
}
 
 
dmpvecarr3E*
MRBF::getDisplacements(){
    if(m_areScalarResults)  return nullptr;
    else                    return &m_displ;
};
 
dmpvector1D*
MRBF::getScalarDisplacements(){
    if(m_areScalarResults)  return &m_scalarDispl ;
    else                    return nullptr;
};
 
 
bool
MRBF::areResultsInScalarMode(){
    return m_areScalarResults;
}
 
 
bool
MRBF::isCompact(){
    return m_isCompact;
}
 
int
MRBF::addNode(darray3E node){
    return(RBF::addNode(node));
};
 
std::vector<int>
MRBF::addNode(dvecarr3E nodes){
    return(RBF::addNode(nodes));
};
 
void
MRBF::setNode(darray3E node){
    removeAllNodes();
    RBF::addNode(node);
    m_rbfgeometry = nullptr;
};
 
void
MRBF::setNode(dvecarr3E nodes){
    removeAllNodes();
    RBF::addNode(nodes);
    m_rbfgeometry = nullptr;
};
 
void
MRBF::setNode(MimmoSharedPointer<MimmoObject> geometry){
    if(!geometry)    return ;
    m_rbfgeometry = geometry;
};
 
void
MRBF::setFilter(dmpvector1D * filter){
    if(!filter) return;
    m_filter.clear();
    m_bfilter = !(filter->empty());
    m_filter = *filter;
};
 
ivector1D
MRBF::checkDuplicatedNodes(double tol){
    ivector1D marked;
    int sizeEff = getTotalNodesCount();
    if( sizeEff == 0 ) return marked;
 
    bvector1D check(sizeEff, false);
 
    for(int i=0; i<sizeEff; ++i){
        for(int j=i+1; j<sizeEff; ++j){
            double dist = norm2(m_node[j] - m_node[i]);
            if(!check[j] && dist <= tol){
                marked.push_back(i);
                check[j] = true;
            }
        }
    }
    return(marked);
}
 
bool
MRBF::removeDuplicatedNodes(ivector1D * list){
    ivector1D marked;
    if(list==nullptr){
        marked = checkDuplicatedNodes();
        list = &marked;
    }
    return(removeNode(*list));
}
 
 
void
MRBF::setSupportRadiusLocal(double suppR_){
    suppR_ = std::max(-1.0,suppR_);
    m_supportRadiusValue = suppR_;
    m_srIsReal = false;
    m_supportRadii.clear();
}
 
 
void
MRBF::setSupportRadiusReal(double suppR_){
    suppR_ = std::max(-1.0,suppR_);
    m_supportRadiusValue = suppR_;
    m_srIsReal = true;
    m_supportRadii.clear();
}
 
void
MRBF::setSupportRadiusValue(double suppR_){
    setSupportRadiusReal(suppR_);
}
 
void
MRBF::setVariableSupportRadii(dvector1D sradii){
    if(sradii.empty() || m_solver != MRBFSol::NONE) return;
    m_supportRadii = sradii;
    m_supportRadiusValue = -1.0;
    m_srIsReal = false;
}
 
void
MRBF::setVariableSupportRadii(dmpvector1D* sradii){
    if(!sradii || m_solver != MRBFSol::NONE) return;
    m_rbfSupportRadii = sradii;
    m_supportRadiusValue = -1.0;
    m_srIsReal = false;
}
 
void
MRBF::setDiagonalFactor(double diagonalFactor){
    m_diagonalFactor = std::min(std::max(diagonalFactor, 0.), 1.);
}
 
void
MRBF::setTol(double tol){
    m_tol = tol;
}
 
void
MRBF::setDisplacements(dvecarr3E displ){
    int size = displ.size();
 
    removeAllData();
 
    dvector1D temp(size);
    for(int loc=0; loc<3; ++loc){
        for(int i=0; i<size; ++i){
            temp[i] = displ[i][loc];
        }
        addData(temp);
    }
 
    m_rbfdispl = nullptr;
    m_areScalarResults = false;
}
 
void
MRBF::setDisplacements(dmpvecarr3E* displ){
    if (!displ) return;
    m_rbfdispl = displ;
    m_areScalarResults = false;
}
 
 
void
MRBF::setScalarDisplacements(dvector1D displ){
    removeAllData();
    addData(displ);
 
    m_rbfScalarDispl = nullptr;
    m_areScalarResults = true;
 
}
 
void
MRBF::setScalarDisplacements(dmpvector1D* displ){
    if (!displ) return;
    m_rbfScalarDispl = displ;
}
 
void
MRBF::setFunction( const MRBFBasisFunction &bfunc, bool isCompact )
{
    switch(bfunc){
 
    case( MRBFBasisFunction::HEAVISIDE10):
            RBF::setFunction(heaviside10);
            m_functype = static_cast<int>(MRBFBasisFunction::HEAVISIDE10);
            break;
 
    case( MRBFBasisFunction::HEAVISIDE50):
            RBF::setFunction(heaviside50);
            m_functype = static_cast<int>(MRBFBasisFunction::HEAVISIDE50);
            break;
 
    case( MRBFBasisFunction::HEAVISIDE100):
            RBF::setFunction(heaviside100);
            m_functype = static_cast<int>(MRBFBasisFunction::HEAVISIDE100);
            break;
 
    case( MRBFBasisFunction::HEAVISIDE1000):
            RBF::setFunction(heaviside1000);
            m_functype = static_cast<int>(MRBFBasisFunction::HEAVISIDE1000);
            break;
 
    case( MRBFBasisFunction::DSIGMOID):
            RBF::setFunction(dsigmoid);
            m_functype = static_cast<int>(MRBFBasisFunction::DSIGMOID);
            break;
 
    default:
            bitpit::RBF::setFunction(bitpit::RBFBasisFunction::WENDLANDC2);
            m_functype = -1;
            break;
    }
 
    m_isCompact = isCompact;
}
 
void
MRBF::setFunction( const bitpit::RBFBasisFunction &bfunc, bool isCompact)
{
    bitpit::RBF::setFunction(bfunc);
    m_functype = -1;
    m_isCompact = isCompact;
}
 
void
MRBF::setCompactSupport(bool isCompact)
{
    m_isCompact = isCompact;
}
 
void
MRBF::clear(){
    BaseManipulation::clear();
    clearFilter();
    m_tol = 0.00001;
    m_supportRadiusValue = -1.0;
    m_srIsReal = false;
    m_supportRadii.clear();
};
 
void
MRBF::clearFilter(){
    m_filter.clear();
    m_bfilter = false;
};
 
 
 
void
MRBF::execute(){
 
    MimmoSharedPointer<MimmoObject> container = getGeometry();
    if(container == nullptr){
        (*m_log)<<m_name + " : nullptr pointer to linked geometry found"<<std::endl;
        throw std::runtime_error(m_name + "nullptr pointer to linked geometry found");
    }
 
    // If RBF nodes passed by MimmoObject initialize RBFs
    if (m_rbfgeometry){
        if (!initRBFwGeometry())
        {
            // Initialization not valid (log message written in method). Skip execution.
            return;
        }
    }
 
    //prepare m_displ or m_scalarDispl according to m_areScalarResults;
    m_displ.clear();
    m_scalarDispl.clear();
    if(m_areScalarResults){
        m_scalarDispl.setDataLocation(mimmo::MPVLocation::POINT);
        m_scalarDispl.reserve(getGeometry()->getNVertices());
        m_scalarDispl.setGeometry(getGeometry());
    }else{
        m_displ.setDataLocation(mimmo::MPVLocation::POINT);
        m_displ.reserve(getGeometry()->getNVertices());
        m_displ.setGeometry(getGeometry());
    }
 
    //resize displacements.
    int size = 0;
    int sizeF = getDataCount();
 
    for (int i=0; i<sizeF; i++){
 
        if(m_solver == MRBFSol::NONE)    size = m_weight[i].size();
        else                            size = m_value[i].size();
 
        if(size != getTotalNodesCount()){
            (*m_log) << "warning: " << getName() << " has displacements of " << i << " field with size (" << size << ") that does not fit number of RBF nodes ("<< getTotalNodesCount() << ")" << std::endl;
            fitDataToNodes(i);
        }
    }
 
    //compute the support radius in m_effectiveSR
    //and push homogeneous support radius info to the base class,
    //in case of Mode WHOLE/GREEDY
    computeEffectiveSupportRadiusList();
 
    //calculate weights for interpolation modes. This is not required
    // in parameterization mode MRBFSol::NONE.
    if (m_solver == MRBFSol::WHOLE)    solve();
    if (m_solver == MRBFSol::GREEDY)    greedy(m_tol);
 
 
    // Prepare the list of vertices to be used during rbf evaluations
    std::unordered_set<long> activeMeshVertices;
 
    // Use bounding box to decide if use the whole geometry or filter the vertices
    std::array<double,3> boxmin, boxmax;
    container->getBoundingBox(boxmin, boxmax);
    // Use diagonal length vs. maximum support radius
    double dlength = norm2(boxmax-boxmin);
    double maximumRadius = 0.;
    for (double radius : getEffectivelyUsedSupportRadii())
    {
        if (radius > maximumRadius){
            maximumRadius = radius;
        }
    }
    bool useWholeGeometry = false;
    // Use the whole geometry if maximum radius is greater than the diagonal of bounding box multiplied for a custom factor
    if (maximumRadius > (m_diagonalFactor*dlength))
    {
        useWholeGeometry = true;
    }
 
    // Set the active vertices of target geometry
    if (isCompact() && !useWholeGeometry){
        // Fill the list of vertices with them included in rbf radii
        if(container->getKdTreeSyncStatus() != SyncStatus::SYNC){
            container->buildKdTree();
        }
        int nnodes = getTotalNodesCount();
        std::vector<long> ids;
        bitpit::KdTree<3,bitpit::Vertex,long>* tree = container->getKdTree();
        for(int i=0; i<nnodes; ++i){
            bitpit::Vertex vertexNode(bitpit::Vertex::NULL_ID, m_node[i]);
            tree->hNeighbors(&vertexNode, m_effectiveSR[i], &ids, nullptr);
            activeMeshVertices.insert(ids.begin(), ids.end());
        }
    }
    else{
        // Use all the geometry vertices
        std::vector<long> ids = container->getVerticesIds();
        activeMeshVertices.insert(ids.begin(), ids.end());
    }
 
    // get deformation using own class evalRBF.
    std::array<double,3> tempValue;
    std::vector<double>  resultValue;
    for(const long &id: activeMeshVertices){
        bitpit::Vertex & vertex = container->getPatch()->getVertex(id);
        resultValue =evalRBF(vertex.getCoords());
        if(m_areScalarResults) {
            m_scalarDispl.insert(id, resultValue[0]);
        }else{
            std::copy_n(resultValue.begin(), 3,tempValue.begin());
            m_displ.insert(id, tempValue);
        }
    }
 
    //apply m_filter if it's active;
    if(m_bfilter){
        checkFilter();
        if(m_areScalarResults){
            for (auto it=m_scalarDispl.begin(); it!=m_scalarDispl.end(); ++it){
                (*it) *= m_filter.at(it.getId());
            }
        }else{
            for (auto it=m_displ.begin(); it!=m_displ.end(); ++it){
                (*it) *= m_filter.at(it.getId());
            }
        }
    }
 
 
    if(m_areScalarResults){
        m_scalarDispl.completeMissingData(0.0);
    }else{
        m_displ.completeMissingData({{0.0,0.0,0.0}});
    }
};
 
void
MRBF::apply(){
    if(m_geometry == nullptr) return;
 
    if(m_areScalarResults){
 
        //I can directly apply the "scalar" displacement field only in case of surface meshes.
        // check it out if targeted geometry is a surface.
        if(m_geometry->getType() != 1){
            (*m_log)<<"WARNING "<<m_name<<" : while in SCALAR mode, the class can apply a displacement field only in case of surface geometry (MimmoObject type 1). Skip apply(). "<<std::endl;
            return;
        }
        m_geometry->updateAdjacencies();
        //retrieve local normals and recover 3D comp m_displ
        // multiplying them for the scalar field input vertex by vertex.
        bitpit::SurfaceKernel* skernel = static_cast<bitpit::SurfaceKernel*>(m_geometry->getPatch());
        m_displ.clear();
        m_displ.reserve(m_geometry->getNVertices());
        bitpit::ConstProxyVector<long> verts;
        std::size_t size;
        long idN;
        for(const bitpit::Cell & cell: m_geometry->getCells()){
            verts = cell.getVertexIds();
            size = verts.size();
            for(std::size_t i=0; i<size; ++i){
                idN = verts[i];
                if(!m_displ.exists(idN)){
                    m_displ.insert(idN, skernel->evalVertexNormal(cell.getId(), i));
                    m_displ[idN] *= m_scalarDispl[idN];
                }
            }
        }
    }//endif m_areScalarResults;
 
    _apply(m_displ); //base manipulation utility method.
 
    if(m_areScalarResults)  m_displ.clear();
}
 
 
void
MRBF::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    std::string input;
 
    BaseManipulation::absorbSectionXML(slotXML, name);
 
    m_tol = 1.0E-6;
    if(slotXML.hasOption("Tolerance")){
        input = slotXML.get("Tolerance");
        input = bitpit::utils::string::trim(input);
        double value = m_tol;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss >> value;
            if(value > 0.0)    setTol(value);
        }
    };
 
    if(slotXML.hasOption("RBFShape")){
        input = slotXML.get("RBFShape");
        int value =1;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss>>value;
        }
        value = std::max(1, value);
        if(value > 100){
            if (value > 105){
                value = 105;
            }
            MRBFBasisFunction shapetype = static_cast<MRBFBasisFunction>(value);
            setFunction(shapetype);
        }
        else{
            if(value > 13){
                value =1;
            }
            bitpit::RBFBasisFunction shapetype = static_cast<bitpit::RBFBasisFunction>(value);
            setFunction(shapetype);
        }
    };
 
    if(slotXML.hasOption("SupportRadiusLocal")){
        input = slotXML.get("SupportRadiusLocal");
        double value = -1.0;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss >> value;
            setSupportRadiusLocal(value);
        }
    };
    //LEGACY ENTRY old version of SupportRadiusLocal
    if(slotXML.hasOption("SupportRadius")){
        input = slotXML.get("SupportRadius");
        double value = -1.0;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss >> value;
            setSupportRadiusLocal(value);
        }
    };
 
    if(slotXML.hasOption("SupportRadiusReal")){
        input = slotXML.get("SupportRadiusReal");
        double value = -1.0;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss >> value;
            setSupportRadiusReal(value);
        }
    };
 
    if(slotXML.hasOption("Compact")){
        std::string input = slotXML.get("Compact");
        input = bitpit::utils::string::trim(input);
        bool value = false;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
        }
        setCompactSupport(value);
    }
 
    m_diagonalFactor = 1.0;
    if(slotXML.hasOption("DiagonalFactor")){
        input = slotXML.get("DiagonalFactor");
        input = bitpit::utils::string::trim(input);
        double value = m_diagonalFactor;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss >> value;
            if(value >= 0.0 && value <=1.0)    setDiagonalFactor(value);
        }
    };
 
}
 
void
MRBF::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::flushSectionXML(slotXML, name);
 
    slotXML.set("Mode", std::to_string(static_cast<int>(m_solver)));
 
    if(m_tol != 1.0E-6){
        std::stringstream ss;
        ss<<std::scientific<<m_tol;
        slotXML.set("Tolerance", ss.str());
    }
    int type = getFunctionType();
    if(type != static_cast<int>(bitpit::RBFBasisFunction::CUSTOM)){
        slotXML.set("RBFShape", std::to_string(type));
    }
 
    if(m_supportRadii.empty()){
        std::stringstream ss;
        ss<<std::scientific<<m_supportRadiusValue;
        if(m_srIsReal){
            slotXML.set("SupportRadiusReal", ss.str());
        }else{
            slotXML.set("SupportRadiusLocal", ss.str());
        }
    }
 
    if(isCompact()){
        slotXML.set("Compact", std::to_string(1));
    }
 
    if(m_diagonalFactor != 1.0){
        std::stringstream ss;
        ss<<std::scientific<<m_diagonalFactor;
        slotXML.set("DiagonalFactor", ss.str());
    }
 
}
 
void
MRBF::checkFilter(){
    bool check = m_filter.getDataLocation() == mimmo::MPVLocation::POINT;
    check = check && m_filter.completeMissingData(0.0);
    check = check && m_filter.getGeometry() == getGeometry();
 
    if (!check){
        m_log->setPriority(bitpit::log::Verbosity::DEBUG);
        (*m_log)<<"Not valid filter found in "<<m_name<<". Proceeding with default unitary field"<<std::endl;
        m_log->setPriority(bitpit::log::Verbosity::NORMAL);
 
        m_filter.clear();
        m_filter.setGeometry(m_geometry);
        m_filter.setDataLocation(mimmo::MPVLocation::POINT);
        m_filter.reserve(getGeometry()->getNVertices());
        for (const auto & vertex : getGeometry()->getVertices()){
            m_filter.insert(vertex.getId(), 1.0);
        }
    }
}
 
void
MRBF::setWeight(dvector2D value){
    if(m_solver != MRBFSol::NONE)    return;
 
    int size = value.size();
 
    removeAllData();
 
    dvector1D temp(size);
    int sizeLoc = 0;
    if(!(value.empty()))    sizeLoc = value[0].size();
    for(int loc=0; loc<sizeLoc; ++loc){
        for(int i=0; i<size; ++i){
            temp[i] = value[i][loc];
        }
        addData(temp);
    }
}
 
void
MRBF::plotCloud(std::string directory, std::string filename, int counterFile, bool binary, bool deformed){
 
    if(deformed && m_areScalarResults){
        (*m_log)<<"WARNING "<<m_name<<": plotCloud skipping print of deformed RBF cloud, while the class is in Scalar mode"<<std::endl;
        return;
    }
 
    int nnodes = getTotalNodesCount();
    dvecarr3E* nodes_ = getNodes();
    dvecarr3E nodes(nnodes);
    if(deformed){
        dvecarr3E data(nnodes);
        for(int i=0; i<nnodes; ++i){
            for(int j=0; j<3; ++j){
                if(m_solver == MRBFSol::NONE)   data[i][j] = m_weight[j][i];
                else                            data[i][j] = m_value[j][i];
            }
        }
        for(int i=0; i<nnodes; ++i){
            nodes[i] = (*nodes_)[i] + data[i];
        }
    }else{
        for(int i=0; i<nnodes; ++i){
            nodes[i] = (*nodes_)[i];
        }
    }
 
    bitpit::VTKFormat codex = bitpit::VTKFormat::ASCII;
    if(binary){codex=bitpit::VTKFormat::APPENDED;}
 
    ivector1D conn(nnodes);
    {
        int counter = 0;
        for(auto & val: conn){
            val = counter;
            ++counter;
        }
    }
    bitpit::VTKUnstructuredGrid vtk(directory, filename, bitpit::VTKElementType::VERTEX);
 
#if MIMMO_ENABLE_MPI
    // only rank 0 writes (nodes are shared!)
    if (getRank() == 0){
#endif
        vtk.setGeomData( bitpit::VTKUnstructuredField::POINTS, nodes) ;
        vtk.setGeomData( bitpit::VTKUnstructuredField::CONNECTIVITY, conn) ;
        vtk.setDimensions(conn.size(), nnodes);
        //use connectivity as node labels also;
        vtk.addData("labels",bitpit::VTKFieldType::SCALAR, bitpit::VTKLocation::POINT, conn) ;
        vtk.setCodex(codex);
        if(counterFile>=0){vtk.setCounter(counterFile);}
        vtk.write();
#if MIMMO_ENABLE_MPI
    }
#endif
};
 
void
MRBF::plotOptionalResults(){
 
    std::string dir = m_outputPlot;
    std::string nameCloud = m_name;
    std::string nameCloudD = m_name+"_moved";
 
    plotCloud(dir, nameCloud, getId(), true, false );
    plotCloud(dir, nameCloudD, getId(), true, true );
}
 
void
MRBF::computeEffectiveSupportRadiusList(){
 
    m_effectiveSR.clear();
    int nodesCount = getTotalNodesCount();
    //First step check variable support radii.
    //if not empty, pass as it is to m_effectiveSR, resize to the total count of nodes and return.
    // if class is in MRBFSol::WHOLE/GREEDY m_supportRadii should already be empty
    // see control of setVariableSupportRadii
    if(!m_supportRadii.empty()){
        m_effectiveSR = m_supportRadii;
        m_effectiveSR.resize(nodesCount, std::numeric_limits<double>::min());
        return;
    }
 
    // now you are working with singular supportRadius
    double candidateRadius = 0.0;
    //check if it's negative.
    if(m_supportRadiusValue < std::numeric_limits<double>::min()) {
    //get maximum weight/value displ and assign support radius a 3 times this value.
        int dataCount = getDataCount();
        dvector1D data(dataCount);
        for(int i=0; i<nodesCount; ++i){
            for(int j=0; j<dataCount; ++j){
                if(m_solver == MRBFSol::NONE)    data[j] = m_weight[j][i];
                else                             data[j] = m_value[j][i];
            }
            candidateRadius = std::max(candidateRadius, norm2(data));
        }
        candidateRadius *=3.0;
    }else{
 
        candidateRadius = m_supportRadiusValue;
 
        // check if its local(based on AABB)
        if (!m_srIsReal){
            double bboxDiag;
            darray3E pmin, pmax;
            getGeometry()->getBoundingBox(pmin, pmax);
            bboxDiag= norm2(pmax - pmin);
            candidateRadius *= bboxDiag;
        }
    }
 
    //create homogeneous vector of support radii with this candidateRadius;
    m_effectiveSR.resize(nodesCount, candidateRadius);
 
    //If you are in WHOLE/GREEDY mode, push to the class base
    // the info on the candidate radius
    if(m_solver != MRBFSol::NONE){
        setSupportRadius(candidateRadius);
    }
}
 
std::vector<double>
MRBF::evalRBF( const std::array<double,3> &point){
 
    int datasize = getDataCount();
    std::vector<double> values(datasize, 0.0);
    int                 i, j;
    double              dist, basis;
 
    for( i=0; i<m_nodes; ++i ){
        if( m_activeNodes[i] ) {
 
            dist = norm2(point - m_node[i]) / m_effectiveSR[i];
            basis = evalBasis( dist );
 
            for( j=0; j<datasize; ++j) {
                values[j] += basis * m_weight[j][i];
            }
        }
    }
 
    return values;
}
 
void
MRBF::setMode(MRBFSol solver){
    m_solver = solver;
    if (m_solver == MRBFSol::NONE)    RBF::setMode(bitpit::RBFMode::PARAM);
    else                            RBF::setMode(bitpit::RBFMode::INTERP);
};
 
bool
MRBF::initRBFwGeometry(){
 
    // If rbf displ is null and Scalar mode if off skip check on linked geometry
    if(m_rbfdispl != nullptr && !m_areScalarResults){
        // If displacementes have unlinked geometry skip
        if (!m_rbfdispl->getGeometry()){
            (*m_log)<<m_name + " : null RBF geometry linked by displacements DOF vectors. Skip object."<<std::endl;
            return false;
        }
        else if (m_rbfdispl->getGeometry() != m_rbfgeometry){
            (*m_log)<<m_name + " : RBF DOF vector displacements not linked to RBF geometry. Skip object."<<std::endl;
            return false;
        }
    }
 
    // If rbfScalarDispl is null and Scalar mode is on skip check on linked geometry
    if(m_rbfScalarDispl != nullptr && m_areScalarResults){
        // If displacementes have unlinked geometry skip
        if (!m_rbfScalarDispl->getGeometry()){
            (*m_log)<<m_name + " : null RBF geometry linked by displacements DOF scalars. Skip object."<<std::endl;
            return false;
        }
        else if (m_rbfScalarDispl->getGeometry() != m_rbfgeometry){
            (*m_log)<<m_name + " : RBF DOF scalar displacements not linked to RBF geometry. Skip object."<<std::endl;
            return false;
        }
    }
 
    // If rbf support radii is null skip check on linked geometry
    if (m_rbfSupportRadii){
        if (!m_rbfSupportRadii->getGeometry()){
            (*m_log)<<m_name + " : null RBF geometry linked by support radii vector. Skip object."<<std::endl;
            return false;
        }
        else if (m_rbfSupportRadii->getGeometry() != m_rbfgeometry){
            (*m_log)<<m_name + " : RBF support radii not linked to RBF geometry. Skip object."<<std::endl;
            return false;
        }
    }
 
    //clear the previous data into MRBF
    removeAllNodes();
    removeAllData();
 
    //fill data locally from m_rbfgeometry, m_rbfdispl/m_rbfScalarDispl and m_rbfSupportRadii
    std::map<long, std::array<double,3> > nodes, displs;
    std::map<long, double> radii;
    for(const bitpit::Vertex & vert : m_rbfgeometry->getVertices()){
        long id = vert.getId();
        nodes[id] = vert.getCoords();
 
        if(m_areScalarResults){
            if(m_rbfScalarDispl && m_rbfScalarDispl->exists(id)){
                displs[id].fill(0.0);
                displs[id][0] = m_rbfScalarDispl->at(id);
            }
        }else{
            if(m_rbfdispl && m_rbfdispl->exists(id)){
                displs[id] = m_rbfdispl->at(id);
            }
        }
        if (m_rbfSupportRadii && m_rbfSupportRadii->exists(id)){
            radii[id] = m_rbfSupportRadii->at(id);
        }
    }
 
#if MIMMO_ENABLE_MPI
    //MPI stuffs - serialize nodes and displs data so that
    // every rank has exactly the same data chunks.
 
    //Send my own and receive mpvres from others.
 
    //prepare my own send buffer.
    mimmo::OBinaryStream myrankNodeDataBuffer, myrankDisplDataBuffer, myrankRadiiDataBuffer;
    myrankNodeDataBuffer << (std::size_t)nodes.size();
    for (auto& tuple: nodes){
        myrankNodeDataBuffer << tuple.first;
        myrankNodeDataBuffer << tuple.second;
    }
    myrankDisplDataBuffer << (std::size_t)displs.size();
    for (auto& tuple: displs){
        myrankDisplDataBuffer << tuple.first;
        myrankDisplDataBuffer << tuple.second;
    }
    myrankRadiiDataBuffer << (std::size_t)radii.size();
    for (auto& tuple: radii){
        myrankRadiiDataBuffer << tuple.first;
        myrankRadiiDataBuffer << tuple.second;
    }
 
    long myrankNodeDataBufferSize = myrankNodeDataBuffer.getSize();
    long myrankDisplDataBufferSize = myrankDisplDataBuffer.getSize();
    long myrankRadiiDataBufferSize = myrankRadiiDataBuffer.getSize();
 
    for (int sendRank=0; sendRank<m_nprocs; sendRank++){
 
        if (m_rank != sendRank){
 
            std::size_t nP;
            std::array<double,3> val;
            double radius;
            long int id;
 
            // receive node data from other ranks.
            long nodeBufferSize;
            MPI_Recv(&nodeBufferSize, 1, MPI_LONG, sendRank, 900, m_communicator, MPI_STATUS_IGNORE);
            mimmo::IBinaryStream nodeBuffer(nodeBufferSize);
            MPI_Recv(nodeBuffer.data(), nodeBuffer.getSize(), MPI_CHAR, sendRank, 910, m_communicator, MPI_STATUS_IGNORE);
 
            nodeBuffer >> nP;
            for (std::size_t i = 0; i < nP; ++i) {
                nodeBuffer >> id;
                nodeBuffer >> val;
                if(nodes.count(id) == 0)    nodes.insert({{id, val}});
            }
 
            // receive displ data from other ranks.
            long displBufferSize;
            MPI_Recv(&displBufferSize, 1, MPI_LONG, sendRank, 920, m_communicator, MPI_STATUS_IGNORE);
            mimmo::IBinaryStream displBuffer(displBufferSize);
            MPI_Recv(displBuffer.data(), displBuffer.getSize(), MPI_CHAR, sendRank, 930, m_communicator, MPI_STATUS_IGNORE);
 
            displBuffer >> nP;
            for (std::size_t i = 0; i < nP; ++i) {
                displBuffer >> id;
                displBuffer >> val;
                if(displs.count(id) == 0)    displs.insert({{id, val}});
            }
 
            // receive radii data from other ranks.
            long radiiBufferSize;
            MPI_Recv(&radiiBufferSize, 1, MPI_LONG, sendRank, 920, m_communicator, MPI_STATUS_IGNORE);
            mimmo::IBinaryStream radiiBuffer(radiiBufferSize);
            MPI_Recv(radiiBuffer.data(), radiiBuffer.getSize(), MPI_CHAR, sendRank, 930, m_communicator, MPI_STATUS_IGNORE);
 
            radiiBuffer >> nP;
            for (std::size_t i = 0; i < nP; ++i) {
                radiiBuffer >> id;
                radiiBuffer >> radius;
                if(radii.count(id) == 0)     radii.insert({{id, radius}});
            }
 
        }else{
            //send to all other except me the def data.
            for (int recvRank=0; recvRank<m_nprocs; recvRank++){
                if (m_rank != recvRank){
                    MPI_Send(&myrankNodeDataBufferSize, 1, MPI_LONG, recvRank, 900, m_communicator);
                    MPI_Send(myrankNodeDataBuffer.data(), myrankNodeDataBuffer.getSize(), MPI_CHAR, recvRank, 910, m_communicator);
                    MPI_Send(&myrankDisplDataBufferSize, 1, MPI_LONG, recvRank, 920, m_communicator);
                    MPI_Send(myrankDisplDataBuffer.data(), myrankDisplDataBuffer.getSize(), MPI_CHAR, recvRank, 930, m_communicator);
                    MPI_Send(&myrankRadiiDataBufferSize, 1, MPI_LONG, recvRank, 920, m_communicator);
                    MPI_Send(myrankRadiiDataBuffer.data(), myrankRadiiDataBuffer.getSize(), MPI_CHAR, recvRank, 930, m_communicator);
                }
            }
        }
    }// end external sendrank loop
 
#endif
 
    dvecarr3E nodeRawList(nodes.size());
    std::array<std::vector<double>,3> displList;
    displList.fill(dvector1D(nodes.size(), 0.));
 
    // Resize support radii to the effective serialized size
    // Empty if support radii not passed as pierced vector
    m_supportRadii.clear();
    m_supportRadii.resize(radii.size());
 
    int count(0);
    for(auto & tuple : nodes){
        long id = tuple.first;
        nodeRawList[count] = tuple.second;
        if(displs.count(id) > 0){
            for(int k=0; k<3; ++k){
                displList[k][count] = displs[id][k];
            }
        }
        if(radii.count(id) > 0){
            m_supportRadii[count] = radii[id];
        }
        ++count;
    }
 
    RBF::addNode(nodeRawList);
    int numberOfData = 1 + int(!m_areScalarResults)*2;
    for(int loc=0; loc<numberOfData; ++loc){
        addData(displList[loc]);
    }
 
    // Return valid initialization
    return true;
}
 
double
heaviside10( double dist )
{
    return 1./(1.+std::exp(-10.*(1.-dist)));
}
 
double
heaviside50( double dist )
{
    return 1./(1.+std::exp(-50.*(1.-dist)));
}
 
double
heaviside100( double dist )
{
    return 1./(1.+std::exp(-100.*(1.-dist)));
}
 
double
heaviside1000( double dist )
{
    return 1./(1.+std::exp(-1000.*(1.-dist)));
}
 
double
dsigmoid( double dist )
{
    return 4.*(std::exp(-10.*dist)/std::pow((1.+std::exp(-10.*dist)),int(2)));
}
 
}