IOCloudPoints.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 "IOCloudPoints.hpp"
 
namespace mimmo {
 
IOCloudPoints::IOCloudPoints(bool readMode){
    m_name         = "mimmo.IOCloudPoints";
    m_read         = readMode;
    m_template     = false;
    m_dir       = ".";
    m_filename     = m_name+"_source.dat";
};
 
IOCloudPoints::IOCloudPoints(const bitpit::Config::Section & rootXML){
 
    m_name         = "mimmo.IOCloudPoints";
    m_template     = false;
    m_dir       = ".";
    m_filename     = m_name+"_source.dat";
    m_read         = true;
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
 
    std::string fallback_name2 = "1";
    std::string input2 = rootXML.get("IOmode", fallback_name2);
    input2 = bitpit::utils::string::trim(input2);
 
    m_read = bool(std::atoi(input2.c_str()));
 
    if(input == "mimmo.IOCloudPoints"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
 
IOCloudPoints::~IOCloudPoints(){};
 
IOCloudPoints::IOCloudPoints(const IOCloudPoints & other):BaseManipulation(other){
    m_read         = other.m_read;
    m_dir       = other.m_dir;
    m_filename     = other.m_filename;
    m_template = other.m_template;
};
 
IOCloudPoints & IOCloudPoints::operator=(IOCloudPoints other){
    swap(other);
    return *this;
};
 
void IOCloudPoints::swap(IOCloudPoints & x) noexcept
{
    std::swap(m_read        , x.m_read);
    std::swap(m_dir         , x.m_dir);
    std::swap(m_filename    , x.m_filename);
    std::swap(m_template    , x.m_template);
    std::swap(m_labels      , x.m_labels);
    std::swap(m_points      , x.m_points);
    m_scalarfield.swap(x.m_scalarfield);
    m_vectorfield.swap(x.m_vectorfield);
    BaseManipulation::swap(x);
}
 
void
IOCloudPoints::buildPorts(){
    bool built = true;
 
    built = (built && createPortIn<dmpvector1D*, IOCloudPoints>(this, &mimmo::IOCloudPoints::setScalarField, M_SCALARFIELD));
    built = (built && createPortIn<dmpvecarr3E*, IOCloudPoints>(this, &mimmo::IOCloudPoints::setVectorField, M_VECTORFIELD));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, IOCloudPoints>(this, &mimmo::IOCloudPoints::setGeometry, M_GEOM));
 
    built = (built && createPortOut<MimmoSharedPointer<MimmoObject>, IOCloudPoints>(this, &mimmo::IOCloudPoints::getGeometry, M_GEOM));
    built = (built && createPortOut<dmpvector1D*, IOCloudPoints>(this, &mimmo::IOCloudPoints::getScalarField, M_SCALARFIELD));
    built = (built && createPortOut<dmpvecarr3E*, IOCloudPoints>(this, &mimmo::IOCloudPoints::getVectorField, M_VECTORFIELD));
 
    m_arePortsBuilt = built;
}
 
dmpvector1D*
IOCloudPoints::getScalarField(){
    return &m_scalarfield;
};
 
dmpvecarr3E*
IOCloudPoints::getVectorField(){
    return &m_vectorfield;
};
 
bool
IOCloudPoints::isTemplate(){
    return m_template;
};
 
void
IOCloudPoints::setReadDir(std::string dir){
    if(!m_read)    return;
    m_dir = dir;
};
 
void
IOCloudPoints::setReadFilename(std::string filename){
    if(!m_read)    return;
    m_filename = filename;
};
 
void
IOCloudPoints::setWriteDir(std::string dir){
    if(m_read)    return;
    m_dir = dir;
};
 
void
IOCloudPoints::setGeometry(MimmoSharedPointer<MimmoObject> geometry){
    // Check if geometry is a point cloud
    if (geometry->getType() != 3) return;
    m_geometry = geometry;
};
 
void
IOCloudPoints::setWriteFilename(std::string filename){
    if(m_read)    return;
    m_filename = filename;
};
 
void
IOCloudPoints::setScalarField(dmpvector1D* scalarfield){
    if(m_read) return;
    m_scalarfield = *scalarfield;
};
 
void
IOCloudPoints::setVectorField(dmpvecarr3E* vectorfield){
    if(m_read) return;
    m_vectorfield = *vectorfield;
};
 
 
void
IOCloudPoints::setTemplate(bool flag){
    if(m_read) return;
    m_template = flag;
};
 
void
IOCloudPoints::clear(){
    m_vectorfield.clear();
    m_scalarfield.clear();
    m_template     = false;
    m_filename     = m_name+"_source.dat";
    BaseManipulation::clear();
}
 
void
IOCloudPoints::execute(){
    if(m_read){
        read();
    }else{
        write();
    }
};
 
void
IOCloudPoints::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    std::string input;
 
    //checking IOmode
 
    if(slotXML.hasOption("IOmode")){
        input = slotXML.get("IOmode");
        bool value = true;
        if(!input.empty()){
            std::stringstream ss(bitpit::utils::string::trim(input));
            ss>>value;
        }
        if (m_read != value){
            (*m_log)<< "warning in class "<<m_name<<": cannot absorb xml parameters for class IOmode mismatching"<<std::endl;
            throw std::runtime_error (m_name + " : xml absorbing failed ");
        }
    };
 
    BaseManipulation::absorbSectionXML(slotXML, name);
 
    if(m_read){
        if(slotXML.hasOption("ReadDir")){
            std::string input = slotXML.get("ReadDir");
            input = bitpit::utils::string::trim(input);
            setReadDir(input);
        };
 
        if(slotXML.hasOption("ReadFilename")){
            std::string input = slotXML.get("ReadFilename");
            input = bitpit::utils::string::trim(input);
            setReadFilename(input);
        };
    }else{
        if(slotXML.hasOption("WriteDir")){
            std::string input = slotXML.get("WriteDir");
            input = bitpit::utils::string::trim(input);
            setWriteDir(input);
        };
 
        if(slotXML.hasOption("WriteFilename")){
            std::string input = slotXML.get("WriteFilename");
            input = bitpit::utils::string::trim(input);
            setWriteFilename(input);
        };
    }
 
 
    if(slotXML.hasOption("Template")){
        std::string input = slotXML.get("Template");
        input = bitpit::utils::string::trim(input);
        bool temp = false;
        if(!input.empty()){
            std::stringstream ss(input);
            ss>>temp;
        }
        setTemplate(temp);
    };
 
}
 
void
IOCloudPoints::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::flushSectionXML(slotXML, name);
 
    slotXML.set("IOmode", std::to_string(int(m_read)));
 
    if(m_read){
        slotXML.set("ReadDir", m_dir);
        slotXML.set("ReadFilename", m_filename);
    }else{
        slotXML.set("WriteDir", m_dir);
        slotXML.set("WriteFilename", m_filename);
    }
 
    slotXML.set("Template", std::to_string(int(m_template)));
};
 
 
void
IOCloudPoints::plotOptionalResults(){
 
    if (getGeometry() == nullptr) return;
    //Assessing data;
    MimmoPiercedVector<double> f1 = m_scalarfield;
    MimmoPiercedVector<std::array<double,3>> f2 = m_vectorfield;
 
    if (!f1.getGeometry()) f1.setGeometry(getGeometry());
    if (!f2.getGeometry()) f2.setGeometry(getGeometry());
    //force both data to have location on POINT, otherwise all these class is pointless
    f1.setDataLocation(MPVLocation::POINT);
    f2.setDataLocation(MPVLocation::POINT);
    f1.completeMissingData(0.0);
    f2.completeMissingData({{0.0,0.0,0.0}});
 
    if(f1.getName() == f2.getName()){
        std::string temp = f1.getName();
        f1.setName(temp+"_scalar");
 
        temp = f2.getName();
        f2.setName(temp+"_vector");
 
    }
    //Write
    BaseManipulation::write(getGeometry(), f1, f2);
}
 
void
IOCloudPoints::read(){
 
    m_points.clear();
    m_labels.clear();
    m_scalarfield.clear();
    m_vectorfield.clear();
 
    // Instantiate a new MimmoObject of type point cloud
    m_geometry.reset(new MimmoObject(3));
 
#if MIMMO_ENABLE_MPI
    // Leave the reading only to the master rank
    if(getRank() == 0)
#endif
    {
 
        std::unordered_map<long, int> mapP;
        std::ifstream reading;
        std::string source = m_dir+"/"+m_filename;
        reading.open(source.c_str());
        if(reading.is_open()){
 
            std::string keyword, line;
            long label;
            darray3E dtrial;
            double temp;
 
            do{
                line.clear();
                keyword.clear();
                std::getline(reading,line);
                line = bitpit::utils::string::trim(line);
                std::stringstream ss(line);
                ss>>keyword;
                keyword = bitpit::utils::string::trim(keyword);
                if(keyword == "$POINT")    {
 
                    ss>>label;
                    m_labels.push_back(label);
 
                    dtrial.fill(0.0);
                    ss>>dtrial[0]>>dtrial[1]>>dtrial[2];
                    m_points.push_back(dtrial);
                }
 
            }while(!reading.eof());
 
            int counter = 0;
            for(auto &lab :m_labels){
                mapP[lab] = counter;
                ++counter;
            }
 
            m_scalarfield.reserve(m_points.size());
            m_vectorfield.reserve(m_points.size());
 
            reading.clear();
            reading.seekg(0, std::ios::beg);
 
            do{
                line.clear();
                keyword.clear();
                std::getline(reading,line);
                line = bitpit::utils::string::trim(line);
                std::stringstream ss(line);
                ss>>keyword;
                keyword = bitpit::utils::string::trim(keyword);
                if(keyword == "$SCALARF")    {
 
                    ss>>label>>temp;
                    m_scalarfield.insert(mapP[label],temp);
                }
                if(keyword == "$VECTORF")    {
 
                    dtrial.fill(0.0);
                    ss>>label>>dtrial[0]>>dtrial[1]>>dtrial[2];
                    m_vectorfield.insert(mapP[label],dtrial);
                }
 
            }while(!reading.eof());
 
        }else{
            (*m_log)<<"error of "<<m_name<<" : cannot open "<<m_filename<< " requested. Exiting... "<<std::endl;
            throw std::runtime_error (m_name + " : cannot open " + m_filename + " requested. Exiting... ");
        }
 
        reading.close();
 
    }
 
    // Fill geometry
 
    // Fill geometry with point cloud list
    // Only rank 0 read points, so only rank 0 enters the loop and fills its partition
    std::size_t ind = 0;
    for (darray3E & coords : m_points){
        long label = m_labels[ind];
        getGeometry()->addVertex(coords, label);
        getGeometry()->addConnectedCell(livector1D(1,label), bitpit::ElementType::VERTEX, label);
        ind++;
    }
 
    getGeometry()->update();
 
    // Clear points and label structure
    dvecarr3E().swap(m_points);
    livector1D().swap(m_labels);
 
    // Set the mimmo pierced vectors
    m_scalarfield.setGeometry(getGeometry());
    m_scalarfield.setDataLocation(MPVLocation::POINT);
    m_scalarfield.setName("scalarfield");
    m_vectorfield.setGeometry(getGeometry());
    m_vectorfield.setDataLocation(MPVLocation::POINT);
    m_vectorfield.setName("vectorfield");
 
};
 
void
IOCloudPoints::write(){
 
    //assessing data;
    if (!m_scalarfield.getGeometry()) m_scalarfield.setGeometry(getGeometry());
    if (!m_vectorfield.getGeometry()) m_vectorfield.setGeometry(getGeometry());
    //force both data to have location on POINT, otherwise all these class is pointless
    m_scalarfield.setDataLocation(MPVLocation::POINT);
    m_vectorfield.setDataLocation(MPVLocation::POINT);
 
    if(!m_scalarfield.empty())  m_scalarfield.completeMissingData(0.0);
    if(!m_vectorfield.empty()) m_vectorfield.completeMissingData({{0.0,0.0,0.0}});
 
    std::string filename;
    if(m_template){
        filename = "TEMPLATE_"+m_filename;
    }else{
        filename = m_filename;
    }
 
    std::ofstream writing;
    std::string source = m_dir+"/"+m_filename;
    std::string keyT1 = "{", keyT2 = "}";
 
#if MIMMO_ENABLE_MPI
    // Write on the file in a sequential way
    for (int irank = 0; irank < getGeometry()->getProcessorCount(); irank++){
        if(getRank() == irank){
 
            // If parallel and not the first process open file in appending mode
            if (getRank() == 0)
                writing.open(source.c_str());
            else
                writing.open(source.c_str(), std::ofstream::out | std::ofstream::app);
#else
            writing.open(source.c_str());
#endif
 
            if(writing.is_open()){
 
                MimmoSharedPointer<MimmoObject> geometry = getGeometry();
                darray3E coords;
                for(const long & label : geometry->getVerticesIds()){
                    coords = geometry->getVertexCoords(label);
                    writing<<"$POINT"<<'\t'<<label<<'\t'<<coords[0]<<'\t'<<coords[1]<<'\t'<<coords[2]<<std::endl;
                }
                writing<<""<<std::endl;
 
                for(const long & label : geometry->getVerticesIds()){
                    if(m_template){
                        std::string str1 = keyT1+"s"+std::to_string(label)+keyT2;
                        writing<<"$SCALARF"<<'\t'<<label<<'\t'<<str1<<std::endl;
                    }else{
                        if(!m_scalarfield.exists(label)) continue;
                        double val = m_scalarfield[label];
                        writing<<"$SCALARF"<<'\t'<<label<<'\t'<<val<<std::endl;
                    }
                }
                writing<<""<<std::endl;
 
                for(const long & label : geometry->getVerticesIds()){
                    if(m_template){
                        std::string str1 = keyT1+"x"+std::to_string(label)+keyT2;
                        std::string str2 = keyT1+"y"+std::to_string(label)+keyT2;
                        std::string str3 = keyT1+"z"+std::to_string(label)+keyT2;
 
                        writing<<"$VECTORF"<<'\t'<<label<<'\t'<<str1<<'\t'<<str2<<'\t'<<str3<<std::endl;
                    }else{
                        if(!m_vectorfield.exists(label))    continue;
                        std::array<double,3> val = m_vectorfield[label];
                        writing<<"$VECTORF"<<'\t'<<label<<'\t'<<val[0]<<'\t'<<val[1]<<'\t'<<val[2]<<std::endl;
                    }
                }
                writing<<""<<std::endl;
            }else{
                (*m_log)<<"error of "<<m_name<<" : cannot open "<<m_filename<< " requested. Exiting... "<<std::endl;
                throw std::runtime_error (m_name + " : cannot open " + m_filename + " requested. Exiting... ");
            }
 
            writing.close();
#if MIMMO_ENABLE_MPI
        }     // End if rank is the current one that has to write and wait all
 
        MPI_Barrier(getCommunicator());
 
    }// End loop on processors to write sequentially
#endif
 
};
 
 
}