ReconstructScalar.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 "ReconstructFields.hpp"
#include <unordered_set>
namespace mimmo{
 
ReconstructScalar::ReconstructScalar(MPVLocation loc){
    m_name = "mimmo.ReconstructScalar";
    m_overlapCriterium = OverlapMethod::MAX;
    m_loc = loc;
    if(m_loc == MPVLocation::UNDEFINED || m_loc == MPVLocation::INTERFACE) m_loc= MPVLocation::POINT;
}
 
ReconstructScalar::ReconstructScalar(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.ReconstructScalar";
    m_overlapCriterium = OverlapMethod::MAX;
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    std::string fallback_loc  = "-1";
    std::string input_loc = rootXML.get("DataLocation", fallback_loc);
 
    input = bitpit::utils::string::trim(input);
    input_loc = bitpit::utils::string::trim(input_loc);
 
    int loc = std::stoi(input_loc);
    if(loc > 0 && loc < 3){
        m_loc  =static_cast<MPVLocation>(loc);
    }else{
        m_loc = MPVLocation::POINT;
    }
 
    if(input == "mimmo.ReconstructScalar"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
ReconstructScalar::~ReconstructScalar(){}
 
ReconstructScalar::ReconstructScalar(const ReconstructScalar & other):BaseManipulation(other){
    m_loc = other.m_loc;
    m_overlapCriterium = other.m_overlapCriterium;
    m_subpatch = other.m_subpatch;
    m_result = other.m_result;
    m_subresults = other.m_subresults;
}
 
void ReconstructScalar::swap(ReconstructScalar & x ) noexcept
{
    std::swap(m_loc, x.m_loc);
    std::swap(m_overlapCriterium,x.m_overlapCriterium);
    std::swap(m_subpatch, x.m_subpatch);
    std::swap(m_subresults, x.m_subresults);
    //std::swap(m_result, x.m_result);
    m_result.swap(x.m_result);
    BaseManipulation::swap(x);
};
 
 
OverlapMethod
ReconstructScalar::getOverlapCriteriumENUM(){
    return m_overlapCriterium;
}
 
int    ReconstructScalar::getOverlapCriterium(){
    return static_cast<int>(m_overlapCriterium);
}
 
int
ReconstructScalar::getNData(){
    return m_subpatch.size();
}
 
dmpvector1D *
ReconstructScalar::getResultField(){
    return  &m_result;
};
 
std::vector<dmpvector1D *>
ReconstructScalar::getResultFields(){
 
    std::vector<dmpvector1D *> res;
    res.reserve(m_subresults.size());
    for(dmpvector1D & ffref : m_subresults){
        res.push_back(&ffref);
    }
    return(res);
};
 
void
ReconstructScalar::setOverlapCriteriumENUM( OverlapMethod funct){
    setOverlapCriterium(static_cast<int>(funct));
};
 
void
ReconstructScalar::setOverlapCriterium( int funct){
    if(funct <1 ||funct > 4)    return;
    m_overlapCriterium = static_cast<OverlapMethod>(funct);
};
 
 
void
ReconstructScalar::addData( dmpvector1D  * field){
    if(!field)  return;
    if(field->getGeometry()== nullptr && field->getDataLocation() != m_loc) return;
    if(field->getGeometry()->getType()==3 && m_loc==MPVLocation::CELL){
        (*m_log)<<"warning in "<<m_name<<" : trying to add field referred to a Point Cloud, while Class has Data Location referred to CELLS. Do Nothing."<<std::endl;
        return;
    }
    m_subpatch.push_back(*field);
 
};
 
void
ReconstructScalar::removeData(mimmo::MimmoSharedPointer<MimmoObject> patch){
    //in m_subpatch remove progressively all pierced vector elements in last position
    //which links towards a geometry of type patch.
    while(m_subpatch.back().getGeometry() == patch){
        m_subpatch.pop_back();
    }
 
    //start searching from the begin on element linking patch
    std::vector<dmpvector1D>::iterator it = m_subpatch.begin();
 
    while(it != m_subpatch.end()){
        if (it->getGeometry() == patch){
            *it = m_subpatch.back();
            m_subpatch.pop_back();
            //be always sure the last element is clean and does not link
            // the patch
            while(m_subpatch.back().getGeometry() == patch){
                m_subpatch.pop_back();
            }
        }
        ++it; //forward the iterator.
    }
};
 
void
ReconstructScalar::removeAllData(){
    m_subpatch.clear();
    m_result.clear();
    m_subresults.clear();
};
 
void
ReconstructScalar::clear(){
    BaseManipulation::clear();
    removeAllData();
    m_overlapCriterium = OverlapMethod::MAX;
}
 
void
ReconstructScalar::plotData(){
 
    m_result.setName("field");
    write(m_result.getGeometry(), m_result);
 
};
 
void
ReconstructScalar::plotSubData(int i){
 
    std::string originalname = m_name;
    m_name = originalname + "SubPatch"+std::to_string(i);
    m_subresults[i].setName("field");
    write(m_subresults[i].getGeometry(), m_subresults[i]);
    m_name = originalname;
 
};
 
void
ReconstructScalar::execute(){
 
    if(getGeometry() == nullptr){
        (*m_log)<<m_name + " : nullptr pointer to linked geometry found"<<std::endl;
        throw std::runtime_error("Error in "+m_name + " : nullptr pointer to linked geometry found");
    }
 
    //Overlap fields
    m_result.clear();
    m_result.setGeometry(getGeometry());
    m_result.setDataLocation(m_loc);
 
    m_subresults.clear();
 
    std::unordered_set<long> idsTarget;
    {
        livector1D ids = idsGeoDataLocation(getGeometry());
        idsTarget.insert(ids.begin(), ids.end());
    }
 
    bitpit::PiercedVector<int> counter;
    for (int i=0; i<getNData(); i++){
        dmpvector1D* pv = &m_subpatch[i];
        livector1D ids = idsGeoDataLocation(pv->getGeometry());
        for (auto ID: ids){
            if(idsTarget.count(ID)==0) continue;
 
            if (!m_result.exists(ID)){
                m_result.insert(ID, (*pv)[ID]);
                counter.insert(ID, 1);
            }
            else{
                overlapFields(ID, (*pv)[ID]);
                counter[ID]++;
            }
        }
    }
 
    if (m_overlapCriterium == OverlapMethod::AVERAGE){
        long int ID;
        auto itend = m_result.end();
        for (auto it=m_result.begin(); it!=itend; ++it){
            ID = it.getId();
            (*it) = (*it) / double(counter[ID]);
        }
    }
 
    //Update field on whole geometry
    if(!m_result.completeMissingData(0.0) ){
        (*m_log)<<"Warning in "<<m_name<<". It seems a reconstruct field with values uncoherent with target geometry is generated."<<std::endl;
    }
 
    //Create subresults
    m_subresults.resize(getNData());
    for (int i=0; i<getNData(); i++){
        dmpvector1D* pv = &m_subpatch[i];
        m_subresults[i].setGeometry(pv->getGeometry());
        m_subresults[i].setDataLocation(m_loc);
        livector1D ids = idsGeoDataLocation(pv->getGeometry());
        for (auto ID : ids){
            if(idsTarget.count(ID)>0){
                m_subresults[i].insert(ID, m_result[ID]);
            }
        }
    }
}
 
void     ReconstructScalar::plotOptionalResults(){
    plotData();
    for (int i=0; i<getNData(); i++){
        plotSubData(i);
    }
}
 
//DEVELOPERS REMIND if more overlap methods are added refer to this method to implement them
void
ReconstructScalar::overlapFields(long int ID, double & locField){
 
    switch(m_overlapCriterium){
    case OverlapMethod::MAX :
    {
        double actual = m_result[ID];
        double dummy = locField;
 
        if (actual < dummy)
            m_result[ID] = locField;
    }
    break;
    case OverlapMethod::MIN :
    {
        double actual = m_result[ID];
        double dummy = locField;
 
        if (actual > dummy)
            m_result[ID] = locField;
    }
    break;
    case OverlapMethod::AVERAGE :
        m_result[ID] += locField;
        break;
    case OverlapMethod::SUM :
        m_result[ID] += locField;
        break;
    default : //never been reached
        break;
    }
 
};
 
void
ReconstructScalar::buildPorts(){
 
    bool built = true;
 
    //input
    built = (built && createPortIn<mimmo::MimmoSharedPointer<MimmoObject>, ReconstructScalar>(&m_geometry, M_GEOM, true));
    built = (built && createPortIn<dmpvector1D*, ReconstructScalar>(this, &mimmo::ReconstructScalar::addData, M_SCALARFIELD, true));
 
    //output
    built = (built && createPortOut<dmpvector1D*, ReconstructScalar>(this, &ReconstructScalar::getResultField, M_SCALARFIELD));
    built = (built && createPortOut<mimmo::MimmoSharedPointer<MimmoObject>, ReconstructScalar>(&m_geometry, M_GEOM));
    built = (built && createPortOut<std::vector<dmpvector1D*>, ReconstructScalar>(this, &mimmo::ReconstructScalar::getResultFields, M_VECSCALARFIELDS));
    m_arePortsBuilt = built;
};
 
void
ReconstructScalar::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    if(slotXML.hasOption("DataLocation")){
        std::string input = slotXML.get("DataLocation");
        input = bitpit::utils::string::trim(input);
        int temp = -1;
        if(!input.empty()){
            std::stringstream ss(input);
            ss>>temp;
        }
        if(int(m_loc) != temp){
            (*m_log)<<"Error absorbing DataLocation in "<<m_name<<". Class and read locations mismatch"<<std::endl;
            if (temp == 0) (*m_log)<<"XML DataLocation in "<<m_name<<" is set to 0-UNDEFINED"<<std::endl;
            if (temp == 3) (*m_log)<<"XML DataLocation in "<<m_name<<" is set to 3-INTERFACE, not supported for now."<<std::endl;
            throw std::runtime_error (m_name + " : xml absorbing failed.");
        }
    }
    //start absorbing
    BaseManipulation::absorbSectionXML(slotXML, name);
 
    if(slotXML.hasOption("OverlapCriterium")){
        std::string input = slotXML.get("OverlapCriterium");
        input = bitpit::utils::string::trim(input);
        int value = 1;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
            value = std::min(std::max(1, value),4);
        }
        setOverlapCriterium(value);
    }
};
 
void
ReconstructScalar::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
    BaseManipulation::flushSectionXML(slotXML, name);
 
    slotXML.set("DataLocation", std::to_string(int(m_loc)));
    int value = static_cast<int>(m_overlapCriterium);
    slotXML.set("OverlapCriterium", std::to_string(value));
 
};
 
livector1D ReconstructScalar::idsGeoDataLocation(mimmo::MimmoSharedPointer<MimmoObject> geo){
 
    if (m_loc == MPVLocation::POINT) return geo->getVertices().getIds();
    if (m_loc == MPVLocation::CELL)  return geo->getCells().getIds();
    return livector1D(0);
}
 
}