VTKRectilinear.cpp

/*---------------------------------------------------------------------------*\
 *
 *  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 <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include "VTK.hpp"
 
#include "logger.hpp"
 
namespace bitpit{


VTKRectilinearGrid::VTKRectilinearGrid( ) :VTK() {
 
    m_fh.setAppendix( getExtension() );
 
    m_geometry.push_back( VTKField("x_Coord") ) ;
    m_geometry.push_back( VTKField("y_Coord") ) ;
    m_geometry.push_back( VTKField("z_Coord") ) ;
 
    for( auto & field : m_geometry ){
        field.setLocation( VTKLocation::POINT ) ;
        field.setFieldType( VTKFieldType::SCALAR ) ;
        field.setDataType( VTKDataType::Float64 ) ;
        field.setCodification(m_geomCodex);
    }
 
    for (int k = 0; k < 3; ++k) {
        m_localIndex[k]  = {{ -1, -1 }} ;
        m_globalIndex[k] = {{ -1, -1 }} ;
    }
 
    m_dimensions = -1 ;
 
}

VTKRectilinearGrid::VTKRectilinearGrid( const std::string &dir, const std::string &name ) :VTKRectilinearGrid( ) {
 
    setNames( dir, name ) ;
 
}

VTKRectilinearGrid::VTKRectilinearGrid( const std::string &dir, const std::string &name, VTKFormat codex, int n1, int n2, int m1, int m2, int l1, int l2 ) :VTKRectilinearGrid( ) {
 
    setNames( dir, name ) ;
 
    setDimensions( n1, n2, m1, m2, l1, l2) ;
    setGeomCodex( codex ) ;
 
}

VTKRectilinearGrid::VTKRectilinearGrid( const std::string &dir, const std::string &name, VTKFormat codex, int n, int m, int l ) : VTKRectilinearGrid( ) {
 
    setNames( dir, name ) ;
 
    setDimensions( 0, n-1, 0, m-1, 0, l-1) ;
    setGeomCodex( codex ) ;
 
}

VTKRectilinearGrid::VTKRectilinearGrid( const std::string &dir, const std::string &name, VTKFormat codex, int n1, int n2, int m1, int m2 ) :VTKRectilinearGrid( ) {
 
    setNames( dir, name ) ;
 
    setDimensions( n1, n2, m1, m2, 0, 0) ;
    setGeomCodex( codex ) ;
 
}

VTKRectilinearGrid::VTKRectilinearGrid( const std::string &dir, const std::string &name, VTKFormat codex, int n, int m ) : VTKRectilinearGrid( ) {
 
    setNames( dir, name ) ;
 
    setDimensions( 0, n-1, 0, m-1, 0, 0) ;
    setGeomCodex( codex ) ;
 
}

void VTKRectilinearGrid::readMetaInformation( ){
 
    std::fstream str;
    std::string line, temp;
 
    std::fstream::pos_type        position;
    std::array<int,6>             extensions ;
 
 
    str.open( m_fh.getPath( ), std::ios::in ) ;
 
    getline( str, line);
    while( ! bitpit::utils::string::keywordInString( line, "<VTKFile")){
        getline(str, line);
    }
 
    if( bitpit::utils::string::getAfterKeyword( line, "header_type", '\"', temp) ){
        setHeaderType( temp) ;
    }
 
    while( ! bitpit::utils::string::keywordInString( line, "<Piece")){
        getline(str, line);
    }
 
    bitpit::utils::string::getAfterKeyword( line, "Extent", '\"', temp) ;
    bitpit::utils::string::convertString( temp, extensions );
 
    m_localIndex[0][0] = extensions[0] ;
    m_localIndex[0][1] = extensions[1] ;
    m_localIndex[1][0] = extensions[2] ;
    m_localIndex[1][1] = extensions[3] ;
    m_localIndex[2][0] = extensions[4] ;
    m_localIndex[2][1] = extensions[5] ;
 
    position = str.tellg() ;
 
    readDataHeader( str ) ;
 
    // Read metadata information
    for( auto &field : m_geometry ){ //int i=0; i<geometry.size(); ++i){
        str.seekg( position) ;
        if( ! readDataArray( str, field ) ) {
#if BITPIT_ENABLE_DEBUG
            log::cout() <<"Geometry field " << field.getName() << " not found, it will be disabled" << std::endl ;
#endif
            field.disable();
        }
    }
 
 
    setDimensions( m_localIndex[0][0], m_localIndex[0][1], m_localIndex[1][0], m_localIndex[1][1], m_localIndex[2][0], m_localIndex[2][1] ) ;
    str.close() ; 
 
}

void VTKRectilinearGrid::writeMetaInformation( ) const {
 
    std::fstream str;
 
    str.open( m_fh.getPath( ), std::ios::out ) ;
    if (!str.is_open()) {
        throw std::runtime_error("Cannot create file \"" + m_fh.getName() + "\"" + " inside the directory \"" + m_fh.getDirectory() + "\"");
    }
 
    //Writing XML header
    str << "<?xml version=\"1.0\"?>" << std::endl;
 
    //Writing Piece Information
    str << "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\"  header_type=\"" << m_headerType << "\">" << std::endl; 
    str << "  <RectilinearGrid WholeExtent= \"" 
        << m_globalIndex[0][0] << " " << m_globalIndex[0][1]<< " "
        << m_globalIndex[1][0] << " " << m_globalIndex[1][1]<< " "
        << m_globalIndex[2][0] << " " << m_globalIndex[2][1]<< " "
        << "\" >" << std::endl;
 
    str << "    <Piece Extent= \" " 
        << m_localIndex[0][0] << " " << m_localIndex[0][1]<< " "
        << m_localIndex[1][0] << " " << m_localIndex[1][1]<< " "
        << m_localIndex[2][0] << " " << m_localIndex[2][1]<< " "
        << "\" >" << std::endl;
 
 
 
    //Header for Data
    writeDataHeader( str, false ) ;
 
    //Wring Geometry Information   
    str << "       <Coordinates>" << std::endl;
    writeDataArray( str, m_geometry[0] ) ;
    writeDataArray( str, m_geometry[1] ) ;
    writeDataArray( str, m_geometry[2] ) ;
    str << "       </Coordinates>" << std::endl;
 
    //Closing Piece
    str << "    </Piece>" << std::endl;
    str << "  </RectilinearGrid>" << std::endl;
 
    //Write Appended Section
    if (isAppendedActive()) {
        str << "  <AppendedData encoding=\"raw\">" << std::endl;
        str << "_" << std::endl;
    }
 
    //Closing XML
    str << "</VTKFile>" << std::endl;
 
    str.close() ;
 
}

void VTKRectilinearGrid::writeCollection( const std::string &outputName, const std::string &collectionName ) const {
 
    // Only one process in a parallel output should write the collection
    if (m_procs <= 1 || m_rank != 0) {
        return;
    }
 
    // Create file handler
    FileHandler fhp = createCollectionHandler(collectionName) ;
 
    // Initialize outout stream
    std::fstream str ;
    str.open( fhp.getPath( ), std::ios::out ) ;
    if (!str.is_open()) {
        throw std::runtime_error("Cannot create file \"" + fhp.getName() + "\"" + " inside the directory \"" + fhp.getDirectory() + "\"");
    }
 
    //Writing XML header
    str << "<?xml version=\"1.0\"?>" << std::endl;
 
    //Writing Piece Information
    str << "<VTKFile type=\"PRectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
    str << "  <PRectilinearGrid WholeExtent= \"" 
        << m_globalIndex[0][0] << " " << m_globalIndex[0][1]<< " "
        << m_globalIndex[1][0] << " " << m_globalIndex[1][1]<< " "
        << m_globalIndex[2][0] << " " << m_globalIndex[2][1]<< " "
        << "GhostLevel=\"0\">" << std::endl;
 
 
 
    //Header for Data
    writeDataHeader( str, true );
 
    //Wring Geometry Information
    str << "      <PCoordinates>" << std::endl;
    writePDataArray( str, m_geometry[0] ) ;
    writePDataArray( str, m_geometry[1] ) ;
    writePDataArray( str, m_geometry[2] ) ;
    str << "      </PCoordinates>" << std::endl;
 
 
    FileHandler fho(m_fh) ;
    fho.setSeries(false) ;
    fho.setDirectory(".") ;
    fho.setName(outputName) ;
    for( int i=0; i<m_procs; i++){
        fho.setBlock(i) ;
        const extension3D_t &index = m_procIndex[i] ;
 
        str << "    <Piece Extent= \" " 
            << index[0][0] << " " << index[0][1] << " "
            << index[1][0] << " " << index[1][1] << " "
            << index[2][0] << " " << index[2][1] << " "
            << "\" Source= \"" << fho.getPath() << "\"/>" << std::endl;
    }
 
    str << "  </PRectilinearGrid>"  << std::endl;
    str << "</VTKFile>" << std::endl;
 
    str.close() ;
 
}

void VTKRectilinearGrid::setDimensions( int n1, int n2, int m1, int m2, int l1, int l2 ){
 
    m_dimensions = ( l1 == l2 ) ? 2 :3 ;
 
    m_localIndex[0][0] = n1 ;
    m_localIndex[0][1] = n2 ;
    m_localIndex[1][0] = m1 ;
    m_localIndex[1][1] = m2 ;
    m_localIndex[2][0] = l1 ;
    m_localIndex[2][1] = l2 ;
 
    if(m_procs==1)
        m_globalIndex = m_localIndex ;
 
    m_cells = 1 ;
    m_points = 1 ;
 
    for(int d=0; d<m_dimensions; ++d){
        m_cells = m_cells *  ( m_localIndex[d][1] -m_localIndex[d][0] ) ;
        m_points = m_points *  ( m_localIndex[d][1] -m_localIndex[d][0] +1 ) ;
    }
 
}

void VTKRectilinearGrid::setDimensions( int n, int m, int l ){
 
    this->setDimensions( 0, n-1, 0, m-1, 0, l-1 );
 
}

void VTKRectilinearGrid::setDimensions( int n1, int n2, int m1, int m2 ){
 
    this->setDimensions( n1, n2, m1, m2, 0, 0 );
 
}

void VTKRectilinearGrid::setDimensions( int n, int m ){
 
    this->setDimensions( 0, n-1, 0, m-1, 0, 0 );
 
}

void VTKRectilinearGrid::setGlobalDimensions( int I, int J, int K ){
 
    m_globalIndex[0][0] = 0;
    m_globalIndex[0][1] = I;
    m_globalIndex[1][0] = 0;
    m_globalIndex[1][1] = J;
    m_globalIndex[2][0] = 0;
    m_globalIndex[2][1] = K;
 
}

void VTKRectilinearGrid::setGeomData( VTKRectilinearField fieldEnum, VTKBaseStreamer *streamer ){
 
    int      index = static_cast<int>(fieldEnum) ;
    VTKField& field = m_geometry[index] ;
 
    field.setStreamer( *streamer ) ;
 
}

void VTKRectilinearGrid::setGlobalDimensions( int I, int J ){
 
    m_globalIndex[0][0] = 0;
    m_globalIndex[0][1] = I;
    m_globalIndex[1][0] = 0;
    m_globalIndex[1][1] = J;
    m_globalIndex[2][0] = 0;
    m_globalIndex[2][1] = 0;
 
}

void VTKRectilinearGrid::setGlobalIndex( const std::vector<extension3D_t> &loc ){
 
    if( loc.size() != m_procs ) 
        log::cout() << "Size of loc_ in VTKRectilinearGrid::setParallelIndex does not fit m_procs " << std::endl ;
 
    m_procIndex   = loc ;
 
}

void VTKRectilinearGrid::setGlobalIndex( const std::vector<extension2D_t> &loc ){
 
    if( loc.size() !=m_procs ) log::cout() << "Size of loc_ in VTKRectilinearGrid::setParallelIndex does not fit m_procs " << std::endl ;
 
    m_procIndex.resize(m_procs) ;
 
    for( int i=0; i< m_procs; ++i){
        m_procIndex[i][0]   = loc[i][0] ;
        m_procIndex[i][1]   = loc[i][1] ;
        m_procIndex[i][2]   = {{0,0}} ;
    }
 
}

uint64_t VTKRectilinearGrid::calcFieldSize( const VTKField &field ) const {
 
    uint64_t bytes = calcFieldEntries(field) ;
    bytes *= VTKTypes::sizeOfType( field.getDataType() ) ;
 
    return bytes ;
 
}

uint64_t VTKRectilinearGrid::calcFieldEntries( const VTKField &field ) const {
 
    uint64_t entries(0) ;
    const std::string &name = field.getName() ;
 
    if( name == "x_Coord" ){
        entries = m_localIndex[0][1] -m_localIndex[0][0] +1 ;
 
    } else if( name == "y_Coord" ){
        entries = m_localIndex[1][1] -m_localIndex[1][0] +1 ;
 
    } else if( name == "z_Coord" ){
        entries = m_localIndex[2][1] -m_localIndex[2][0] +1 ;
 
    } else{
 
        VTKLocation location( field.getLocation() ) ;
        assert( location != VTKLocation::UNDEFINED) ;
 
        if( location == VTKLocation::CELL ){
            entries = m_cells ;
 
        } else if( location == VTKLocation::POINT ){
            entries = m_points ;
 
        }
 
        entries *= field.getComponentCount() ;
 
    }
 
    return entries ;
 
}

uint8_t VTKRectilinearGrid::calcFieldComponents( const VTKField &field ) const {
 
    uint8_t comp ;
    const std::string &name = field.getName() ;
 
    if( name == "x_Coord" || name == "y_Cooord" || name == "z_Coord" ){
        comp = 1 ;
 
    } else{
        comp = field.getComponentCount();
 
    }
 
    return comp ;
 
}

std::string VTKRectilinearGrid::getExtension() const {
    return "vtr";
}
 
}