ReconstructVector.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{
 
ReconstructVector::ReconstructVector(MPVLocation loc){
    m_name = "mimmo.ReconstructVector";
    m_overlapCriterium = OverlapMethod::SUM;
    m_loc = loc;
    if(m_loc == MPVLocation::UNDEFINED || m_loc == MPVLocation::INTERFACE) m_loc= MPVLocation::POINT;
}
 
ReconstructVector::ReconstructVector(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.ReconstructVector";
    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.ReconstructVector"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
ReconstructVector::~ReconstructVector(){}
 
ReconstructVector::ReconstructVector(const ReconstructVector & 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 ReconstructVector::swap(ReconstructVector & 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
ReconstructVector::getOverlapCriteriumENUM(){
    return m_overlapCriterium;
}
 
int
ReconstructVector::getOverlapCriterium(){
    return static_cast<int>(m_overlapCriterium);
}
 
int
ReconstructVector::getNData(){
    return m_subpatch.size();
}
 
dmpvecarr3E *
ReconstructVector::getResultField(){
    return  &m_result;
};
 
std::vector<dmpvecarr3E*>
ReconstructVector::getResultFields(){
 
    std::vector<dmpvecarr3E *> res;
    res.reserve(m_subresults.size());
    for(dmpvecarr3E & ffref : m_subresults){
        res.push_back(&ffref);
    }
    return(res);
};
 
void
ReconstructVector::setOverlapCriteriumENUM( OverlapMethod funct){
    setOverlapCriterium(static_cast<int>(funct));
};
 
void
ReconstructVector::setOverlapCriterium( int funct){
    if(funct<1 ||funct>4)    return;
    m_overlapCriterium = static_cast<OverlapMethod>(funct);
};
 
void
ReconstructVector::addData(dmpvecarr3E *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
ReconstructVector::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<dmpvecarr3E>::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
ReconstructVector::removeAllData(){
    m_subpatch.clear();
    m_result.clear();
    m_subresults.clear();
};
 
void
ReconstructVector::clear(){
    BaseManipulation::clear();
    removeAllData();
    m_overlapCriterium = OverlapMethod::SUM;
}
 
void
ReconstructVector::plotData(){
 
    m_result.setName("field");
    write(m_result.getGeometry(), m_result);
 
};
 
void
ReconstructVector::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
ReconstructVector::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++){
        dmpvecarr3E* 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();
            for (int j=0; j<3; j++)
                (*it)[j] = (*it)[j] / double(counter[ID]);
        }
    }
 
    //Update field on whole geometry
    darray3E zero = {{0.0,0.0,0.0}};
 
    if(!m_result.completeMissingData(zero)){
        (*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++){
        dmpvecarr3E* 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
ReconstructVector::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
ReconstructVector::overlapFields(long int ID, darray3E & locField){
 
    switch(m_overlapCriterium){
    case OverlapMethod::MAX :
 
    {
        double actual = norm2(m_result[ID]);
        double dummy = norm2(locField);
 
        if (actual < dummy)
            m_result[ID] = locField;
    }
    break;
 
    case OverlapMethod::MIN :
    {
        double actual = norm2(m_result[ID]);
        double dummy = norm2(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
ReconstructVector::buildPorts(){
 
    bool built = true;
 
    //input
    built = (built && createPortIn<mimmo::MimmoSharedPointer<MimmoObject>, ReconstructVector>(&m_geometry, M_GEOM, true));
    built = (built && createPortIn<dmpvecarr3E*, ReconstructVector>(this, &mimmo::ReconstructVector::addData, M_VECTORFIELD, true));
 
    //output
    built = (built && createPortOut<dmpvecarr3E*, ReconstructVector>(this, &ReconstructVector::getResultField, M_VECTORFIELD));
    built = (built && createPortOut<mimmo::MimmoSharedPointer<MimmoObject>, ReconstructVector>(&m_geometry, M_GEOM));
    built = (built && createPortOut<std::vector<dmpvecarr3E*>, ReconstructVector>(this, &mimmo::ReconstructVector::getResultFields, M_VECVECTORFIELDS));
    m_arePortsBuilt = built;
};
 
void ReconstructVector::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 = 4;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
            value = std::min(std::max(1, value),4);
        }
        setOverlapCriterium(value);
    }
};
 
void ReconstructVector::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 ReconstructVector::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);
}
 
 
}