VTKUnstructured.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{


void VTKUnstructuredGrid::HomogeneousInfoStreamer::setElementType( VTKElementType type){
 
    m_type = type;
 
}

void VTKUnstructuredGrid::HomogeneousInfoStreamer::setCellCount( uint64_t n){
 
    m_nCells = n ;
 
}

void VTKUnstructuredGrid::HomogeneousInfoStreamer::setConnectivityField( const VTKField *connectivity){
 
    m_connectivity = connectivity ;
 
}

void VTKUnstructuredGrid::HomogeneousInfoStreamer::flushData( std::fstream &str, const std::string &name, VTKFormat format){
 
    assert( m_type != VTKElementType::UNDEFINED ) ;
 
    VTKDataType connectivityDataType = m_connectivity->getDataType();
 
    if(name == "types" ){
        uint8_t type = (uint8_t) m_type ;
        for( uint64_t i=0; i<m_nCells; ++i){
            flushValue(str, format, type );
        }
 
    } else if(name == "offsets" ){
        uint8_t     n = vtk::getElementNodeCount(m_type) ;
        if (connectivityDataType == VTKDataType::Int64 || connectivityDataType == VTKDataType::UInt64) {
            uint64_t    offset(0) ;
            for( uint64_t i=0; i<m_nCells; ++i){
                offset += n ;
                flushValue(str, format, offset );
            }
        } else if (connectivityDataType == VTKDataType::Int32 || connectivityDataType == VTKDataType::UInt32) {
            uint32_t    offset(0) ;
            for( uint64_t i=0; i<m_nCells; ++i){
                offset += n ;
                flushValue(str, format, offset );
            }
        } else if (connectivityDataType == VTKDataType::Int16 || connectivityDataType == VTKDataType::UInt16) {
            uint16_t    offset(0) ;
            for( uint64_t i=0; i<m_nCells; ++i){
                offset += n ;
                flushValue(str, format, offset );
            }
        } else if (connectivityDataType == VTKDataType::Int8 || connectivityDataType == VTKDataType::UInt8) {
            uint8_t    offset(0) ;
            for( uint64_t i=0; i<m_nCells; ++i){
                offset += n ;
                flushValue(str, format, offset );
            }
        }
 
    }
 
}


VTKUnstructuredGrid::VTKUnstructuredGrid( VTKElementType elementType ) :VTK() {
 
    m_fh.setAppendix(getExtension());
 
    int nGeomFields = 6;
 
    m_geometry.resize(nGeomFields);
 
    m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)]       = VTKField("Points") ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::OFFSETS)]      = VTKField("offsets") ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::TYPES)]        = VTKField("types") ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::CONNECTIVITY)] = VTKField("connectivity") ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)] = VTKField("faces") ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)] = VTKField("faceoffsets") ;
 
    for( auto & field : m_geometry ){
        field.setLocation( VTKLocation::CELL ) ;
        field.setFieldType( VTKFieldType::SCALAR ) ;
        field.setDataType( VTKDataType::Int32 ) ;
        field.setCodification(m_geomCodex);
    }
 
    m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)].setLocation( VTKLocation::POINT ) ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)].setFieldType( VTKFieldType::VECTOR ) ;
    m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)].setDataType( VTKDataType::Float64 ) ;
 
    m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)].disable();
    m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)].disable();
 
    m_nConnectivityEntries = 0;
    m_nFaceStreamEntries   = 0;
 
    setElementType(elementType);
 
}

VTKUnstructuredGrid::VTKUnstructuredGrid( const std::string &dir, const std::string &name, VTKElementType elementType ):VTKUnstructuredGrid( elementType ){
 
    setNames( dir, name ) ; 
 
}

void VTKUnstructuredGrid::setElementType( VTKElementType elementType ){
 
    // Set element type
    m_elementType = elementType ;
 
    // Set homogeneous info streamer properties
    m_homogeneousInfoStreamer.setElementType( m_elementType );
    if ( m_elementType == VTKElementType::UNDEFINED ) {
        m_homogeneousInfoStreamer.setConnectivityField( nullptr );
        m_homogeneousInfoStreamer.setCellCount( 0 );
    } else {
        m_homogeneousInfoStreamer.setConnectivityField( &(m_geometry[getFieldGeomId(VTKUnstructuredField::CONNECTIVITY)]) );
        m_homogeneousInfoStreamer.setCellCount( m_cells );
    }
 
    // Early return if the grid is not made homogeously of one element type
    if ( m_elementType == VTKElementType::UNDEFINED ) {
        return;
    }
 
    // Types
    int types_gid = getFieldGeomId(VTKUnstructuredField::TYPES);
    m_geometry[types_gid].setDataType( VTKDataType::UInt8) ;
    m_geometry[types_gid].setStreamer(m_homogeneousInfoStreamer) ;
 
    // Offsets
    if ( m_elementType != VTKElementType::POLYGON && m_elementType != VTKElementType::POLYHEDRON) {
        int offsets_gid = getFieldGeomId(VTKUnstructuredField::OFFSETS);
        m_geometry[offsets_gid].setStreamer(m_homogeneousInfoStreamer) ;
    }
 
}

void VTKUnstructuredGrid::setDimensions( uint64_t ncells, uint64_t npoints, uint64_t nconn, uint64_t nfacestream ){
 
    // Grid size
    m_points = npoints ;
 
    m_cells = ncells ;
    if ( m_elementType != VTKElementType::UNDEFINED ) {
        m_homogeneousInfoStreamer.setCellCount( ncells );
    }
 
    // Connectivity
    bool unevenConnectivity = false;
    unevenConnectivity |= ( m_elementType == VTKElementType::UNDEFINED ) ;
    unevenConnectivity |= ( m_elementType == VTKElementType::POLYGON ) ;
    unevenConnectivity |= ( m_elementType == VTKElementType::POLYHEDRON ) ;
 
    if ( unevenConnectivity ) {
        m_nConnectivityEntries = nconn ;
    } else {
        m_nConnectivityEntries = m_cells *vtk::getElementNodeCount( m_elementType ) ;
    }
 
    assert( (m_cells != 0 && m_nConnectivityEntries != 0) || (m_cells == 0 && m_nConnectivityEntries == 0) );
 
    // Face stream
    bool faceSteamEnabled = false;
    faceSteamEnabled |= ( m_elementType == VTKElementType::POLYGON ) ;
    faceSteamEnabled |= ( m_elementType == VTKElementType::POLYHEDRON ) ;
    faceSteamEnabled |= ( m_elementType == VTKElementType::UNDEFINED && nfacestream > 0 ) ;
 
    if ( faceSteamEnabled ) {
        m_nFaceStreamEntries = nfacestream ;
 
        m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)].enable();
        m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)].enable();
    } else {
        m_nFaceStreamEntries = 0 ;
 
        m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)].disable();
        m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)].disable();
    }
 
}

void VTKUnstructuredGrid::setGeomData( VTKUnstructuredField fieldEnum, VTKBaseStreamer *streamer ){
 
    int      index = getFieldGeomId(fieldEnum) ;
    VTKField& field = m_geometry[index] ;
 
    field.setStreamer( *streamer ) ;
 
}

uint64_t VTKUnstructuredGrid::readConnectivityEntries( ){
 
    uint64_t                 nconn(0) ;
 
    std::fstream             str  ;
    std::fstream::pos_type   position_appended;
    std::string              line;
    char                     c_ ;
    uint32_t                 nbytes32 ;
    uint64_t                 nbytes64 ;
 
    // Geometry id of the connectivity
    int connectivity_gid = getFieldGeomId(VTKUnstructuredField::CONNECTIVITY);
    if(connectivity_gid < 0)    return nconn;
 
    if(!m_geometry[connectivity_gid].hasAllMetaData() || !m_geometry[connectivity_gid].isEnabled()){
        throw std::runtime_error("VTKUnstructuredGrid::readConnectivityEntries. Connectivity field is missing meta information or disabled");
    }
 
    if( m_geometry[connectivity_gid].getCodification() == VTKFormat::APPENDED ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary ) ;
        //Go to the initial position of the appended section
        while( getline(str, line) && (! bitpit::utils::string::keywordInString( line, "<AppendedData")) ){}
        str >> c_;
        while( c_ != '_') str >> c_;
        position_appended = str.tellg();
        str.close();
    }
 
    //Open in binary for read
    if( m_geometry[connectivity_gid].getCodification() == VTKFormat::APPENDED ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary);
        str.seekg( position_appended) ;
        str.seekg( m_geometry[connectivity_gid].getOffset(), std::ios::cur) ;
 
        int dataSize = VTKTypes::sizeOfType( m_geometry[connectivity_gid].getDataType() ) ;
        assert(dataSize > 0) ;
 
        if( m_headerType== "UInt32") {
            genericIO::absorbBINARY( str, nbytes32 ) ;
            nconn = nbytes32 / dataSize ;
        }
 
        if( m_headerType== "UInt64") {
            genericIO::absorbBINARY( str, nbytes64 ) ;
            nconn = nbytes64 / dataSize ;
        }
        str.close();
    }
 
 
 
    //Open in ascii for read
    if(  m_geometry[connectivity_gid].getCodification() == VTKFormat::ASCII ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary);
        str.seekg( m_geometry[connectivity_gid].getPosition() ) ;
 
        std::string              line ;
        std::vector<uint64_t>    temp;
 
        nconn = 0 ;
 
        getline( str, line) ;
        while( ! bitpit::utils::string::keywordInString(line,"/DataArray") ) {
            temp.clear() ;
            bitpit::utils::string::convertString( line, temp) ;
            nconn += temp.size() ;
            getline( str, line) ;
        }
 
        str.close();
    }
 
    return nconn ;
 
}

uint64_t VTKUnstructuredGrid::readFaceStreamEntries( ){
    uint64_t                 nface(0) ;
    std::fstream             str  ;
    std::fstream::pos_type   position_appended;
    std::string              line;
    char                     c_ ;
    uint32_t                 nbytes32 ;
    uint64_t                 nbytes64 ;
 
    // Geometry id of the facestream
    int facestream_gid = getFieldGeomId(VTKUnstructuredField::FACE_STREAMS);
    if(facestream_gid < 0 ) return nface;
    if(!m_geometry[facestream_gid].hasAllMetaData() || !m_geometry[facestream_gid].isEnabled()) return nface;
 
    if( m_geometry[facestream_gid].getCodification() == VTKFormat::APPENDED ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary ) ;
        //Go to the initial position of the appended section
        while( getline(str, line) && (! bitpit::utils::string::keywordInString( line, "<AppendedData")) ){}
        str >> c_;
        while( c_ != '_') str >> c_;
        position_appended = str.tellg();
        str.close();
    }
 
    //Open in binary for read
    if( m_geometry[facestream_gid].getCodification() == VTKFormat::APPENDED ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary);
        str.seekg( position_appended) ;
        str.seekg( m_geometry[facestream_gid].getOffset(), std::ios::cur) ;
 
        int dataSize = VTKTypes::sizeOfType( m_geometry[facestream_gid].getDataType() ) ;
        assert(dataSize > 0) ;
 
        if( m_headerType== "UInt32") {
            genericIO::absorbBINARY( str, nbytes32 ) ;
            nface = nbytes32 / dataSize ;
        }
 
        if( m_headerType== "UInt64") {
            genericIO::absorbBINARY( str, nbytes64 ) ;
            nface = nbytes64 / dataSize ;
        }
    }
 
    //Open in ASCII for read
    if(  m_geometry[facestream_gid].getCodification() == VTKFormat::ASCII ){
        str.open( m_fh.getPath( ), std::ios::in | std::ios::binary );
        str.seekg( m_geometry[facestream_gid].getPosition() ) ;
 
        std::string              line ;
        std::vector<uint64_t>    temp;
 
        nface = 0 ;
 
        getline( str, line) ;
        while( ! bitpit::utils::string::keywordInString(line,"/DataArray") ) {
            temp.clear() ;
            bitpit::utils::string::convertString( line, temp) ;
            nface += temp.size() ;
            getline( str, line) ;
        }
 
        str.close();
    }
 
    return nface ;
}

void VTKUnstructuredGrid::writeMetaInformation( ) const {
 
    std::fstream str ;
 
    //Since we are writing only ASCII information (without retrieving stream
    //position with seekg/tellg) it is safe to open it in text/ASCII mode.
    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=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"" << m_headerType << "\">" << std::endl;
    str << "  <UnstructuredGrid>"  << std::endl;;
    str << "    <Piece  NumberOfPoints=\"" << m_points << "\" NumberOfCells=\"" << m_cells << "\">" << std::endl;
 
    //Header for Data
    writeDataHeader( str, false );
 
    //Wring Geometry Information
    str << "      <Points>" << std::endl ;;
    writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)] ) ;
    str << "      </Points>" << std::endl;
 
    str << "      <Cells>" << std::endl ;;
    writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::OFFSETS)] ) ;
    writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::TYPES)] ) ;
    writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::CONNECTIVITY)] ) ;
    if (m_nFaceStreamEntries) {
        writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)] ) ;
        writeDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)] ) ;
    }
    str << "      </Cells>" << std::endl;
 
    //Closing Piece
    str << "    </Piece>" << std::endl;
    str << "  </UnstructuredGrid>"  << std::endl;
 
    //Appended Section
    if (isAppendedActive()) {
        str << "  <AppendedData encoding=\"raw\">" << std::endl;
        str << "_" << std::endl;
    }
 
    // Closing XML
    str << "</VTKFile>" << std::endl;
 
    str.close() ;
 
}

void VTKUnstructuredGrid::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=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
    str << "  <PUnstructuredGrid GhostLevel=\"0\">"  << std::endl;;
 
    //Header for Data
    writeDataHeader( str, true );
 
    //Wring Geometry Information
    str << "      <PPoints>" << std::endl;
    writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::POINTS)] ) ;
    str << "      </PPoints>" << std::endl;
 
    str << "      <PCells>" << std::endl;
    writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::OFFSETS)] ) ;
    writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::TYPES)] ) ;
    writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::CONNECTIVITY)] ) ;
    if (m_nFaceStreamEntries) {
        writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_STREAMS)] ) ;
        writePDataArray( str, m_geometry[getFieldGeomId(VTKUnstructuredField::FACE_OFFSETS)] ) ;
    }
    str << "      </PCells>" << 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) ;
        str << "    <Piece  Source=\"" << fho.getPath() <<  "\"/>" << std::endl;
    }
 
    str << "  </PUnstructuredGrid>"  << std::endl;
    str << "</VTKFile>" << std::endl;
 
    str.close() ;
 
 
}

void VTKUnstructuredGrid::readMetaInformation( ){
 
    std::fstream str;
    std::string line, temp;
 
    std::fstream::pos_type        position;
 
    //Although we are not explicit calling seekg/tellg on the stream, we still
    //need to open the file in binary mode because this mode is needed by the
    //for reading data array information.
    str.open( m_fh.getPath( ), std::ios::in | std::ios::binary ) ;
 
    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, "NumberOfPoints", '\"', temp) ;
    bitpit::utils::string::convertString( temp, m_points );
 
    bitpit::utils::string::getAfterKeyword( line, "NumberOfCells", '\"', temp) ;
    bitpit::utils::string::convertString( temp, m_cells );
 
 
    position = str.tellg() ;
    readDataHeader( str ) ;
 
    // Read metadata information
    for( auto &field : m_geometry ){ 
        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();
        }
    }
 
 
    str.close() ;
 
    if( m_elementType == VTKElementType::UNDEFINED) {
        setDimensions( m_cells, m_points, readConnectivityEntries(), readFaceStreamEntries() ) ;
    } else {
        // Metadata information read form file may not match the information
        // set in our own streamer. If the grid is homogeneous, we need to
        // reset all metadata that can't be overwritten.
        setElementType(m_elementType) ;
 
        // Set the dimension of the grid
        setDimensions( m_cells, m_points ) ;
    }
 
 
}

uint64_t VTKUnstructuredGrid::calcConnectivityEntries( ) const {
 
    return calcFieldEntries( m_geometry[getFieldGeomId(VTKUnstructuredField::CONNECTIVITY)] ) ;
}

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

uint64_t VTKUnstructuredGrid::calcFieldEntries( const VTKField &field ) const {
 
    uint64_t entries(0) ;
    const std::string &name = field.getName() ;
 
    if( name == "Points" ){
        entries = m_points * VTKField::getComponentCount(VTKFieldType::VECTOR) ;
 
    } else if( name == "offsets" ){
        entries = m_cells ;
 
    } else if( name == "types" ){
        entries = m_cells ;
 
    } else if( name == "connectivity"){
        entries = m_nConnectivityEntries ;
 
    } else if( name == "faces"){
        entries = m_nFaceStreamEntries ;
 
    } else if( name == "faceoffsets"){
        entries = m_cells ;
 
    } 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 VTKUnstructuredGrid::calcFieldComponents( const VTKField &field ) const {
 
    uint8_t comp ;
    const std::string &name = field.getName() ;
 
    if( name == "Points" ){
        comp = VTKField::getComponentCount(VTKFieldType::VECTOR) ;
 
    } else if( name == "offsets" ){
        comp = 1 ;
 
    } else if( name == "types" ){
        comp = 1 ;
 
    } else if( name == "connectivity" ){
       if( m_elementType != VTKElementType::UNDEFINED){
            comp = vtk::getElementNodeCount( m_elementType ) ;
 
       } else {
           comp = 1;
 
       }
 
    } else if( name == "faces" ){
        comp = 1 ;
 
    } else if( name == "faceoffsets" ){
        comp = 1 ;
 
    } else{
        comp = field.getComponentCount();
 
    }
 
    return comp ;
 
}

int VTKUnstructuredGrid::getFieldGeomId( VTKUnstructuredField field ) const {
 
    return static_cast<std::underlying_type<VTKElementType>::type>(field) ;
}

std::string VTKUnstructuredGrid::getExtension() const {
    return "vtu";
}
 
}