ITHACAstream.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * In real Time Highly Advanced Computational Applications for Finite Volumes
 * Copyright (C) 2017 by the ITHACA-FV authors
-------------------------------------------------------------------------------
License
    This file is part of ITHACA-FV
    ITHACA-FV is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    ITHACA-FV 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 ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
Class
    ITHACAstream
Description
    Methods for input and outputs from and to files
SourceFiles
    ITHACAstream.C
    ITHACAstreamTemplates.C
\*---------------------------------------------------------------------------*/

 
#ifndef ITHACAstream_H
#define ITHACAstream_H
 
#include "fvCFD.H"
#include "IOmanip.H"
#include <stdio.h>
#include <sys/types.h>
#include <dirent.h>
#include <algorithm>
#include <fstream>
#include <string>
#include <stdexcept>
#include <sstream>
#include <vector>
#include <cstdio>
#include <typeinfo>
#include <iostream>
#include <cassert>
#include <zlib.h>
#include <map>
#include "ITHACAassert.H"
#include "ITHACAparameters.H"
#include "ITHACAutilities.H"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
#pragma GCC diagnostic ignored "-Wignored-attributes"
#include <Eigen/Eigen>
#include "EigenFunctions.H"
 
#include <unsupported/Eigen/CXX11/Tensor>
#pragma GCC diagnostic pop
 
#define MAXBUFSIZE (static_cast<int> (1e6))
#define PBSTR "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
#define PBWIDTH 60
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
/*---------------------------------------------------------------------------*\
                    Namespace ITHACAstream Declaration
\*---------------------------------------------------------------------------*/
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
#include "cnpy.H"
#pragma GCC diagnostic pop

namespace ITHACAstream
{
 
//--------------------------------------------------------------------------
template<typename Type>
void exportFvMatrix(fvMatrix<Type>& Matrix, word folder,
                    word MatrixName);
 
//--------------------------------------------------------------------------
template<typename T, int dim>
void exportMatrix(Eigen::Matrix < T, -1, dim > & matrix, word Name,
                  word type = "python", word folder = "./Matrices");
 
//--------------------------------------------------------------------------
void exportMatrix(List <Eigen::MatrixXd>& matrix, word Name,
                  word type = "python", word folder = "./Matrices");
 
//--------------------------------------------------------------------------
template<typename T>
void exportTensor(Eigen::Tensor<T, 3 > tensor, word Name,
                  word type = "python", word folder = "./Matrices");
 
 
//----------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void read_fields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                 Lfield, word Name,
                 fileName casename, int first_snap = 0, int n_snap = 0);
 
 
 
//----------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh >
void read_fields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                 Lfield,
                 GeometricField<Type, PatchField, GeoMesh>& field,
                 fileName casename, int first_snap = 0, int n_snap = 0);
 
//----------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void readMiddleFields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                      Lfield,
                      GeometricField<Type, PatchField, GeoMesh>& field,
                      fileName casename);
 
//----------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void readConvergedFields(
    PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
    GeometricField<Type, PatchField, GeoMesh>& field,
    fileName casename);
 
//----------------------------------------------------------------------

template<class Type, template<class> class PatchField, class GeoMesh>
void read_last_fields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                      Lfield,
                      const GeometricField<Type, PatchField, GeoMesh>& field,
                      const fileName casename);
 
//----------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void readLastFields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                    Lfield,
                    const GeometricField<Type, PatchField, GeoMesh>& field,
                    const fileName casename);
 
//--------------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void exportFields(PtrList<GeometricField<Type, PatchField, GeoMesh >> &
                  field,
                  word folder, word fieldname);
 
//--------------------------------------------------------------------------
/* One has to provide the complete filename with the absolute or relative path */
Eigen::MatrixXd readMatrix(word filename);
 
//--------------------------------------------------------------------------
/* One has to provide the folder containing the matrix files and the filename of the
 the matrix. Since it is stored as a List of matrices each matrix must be stored in a
 different file with the following format:
 matFileName(i)_mat.txt */
List <Eigen::MatrixXd> readMatrix(word folder, word mat_name);
 
//--------------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
                    fileName subfolder, fileName folder,
                    word fieldName);
 
//--------------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh>
void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
                    fileName subfolder, fileName folder);
 
//--------------------------------------------------------------------------
template<typename T>
void exportList(T& list, word folder, word filename);
 
 
//--------------------------------------------------------------------------
template <typename T, int Nrows, typename IND>
void SaveSparseMatrix(Eigen::SparseMatrix<T, Nrows, IND>& m,
                      word folder, word MatrixName);
 
//--------------------------------------------------------------------------
template <typename T, int Nrows, typename IND>
void ReadSparseMatrix(Eigen::SparseMatrix<T, Nrows, IND>& m,
                      word folder, word MatrixName);
 
//--------------------------------------------------------------------------
template <typename MatrixType>
void SaveDenseMatrix(MatrixType& Matrix, word folder, word MatrixName);
 
//--------------------------------------------------------------------------
template <typename MatrixType>
void ReadDenseMatrix(MatrixType& Matrix, word folder, word MatrixName);
 
 
//----------------------------------------------------------------------
template <typename TensorType>
void SaveDenseTensor(TensorType& Tensor, word folder, word MatrixName);
 
 
//----------------------------------------------------------------------
template <typename TensorType>
void ReadDenseTensor(TensorType& Tensor, word folder, word MatrixName);
 
//--------------------------------------------------------------------------
template <typename MatrixType>
void SaveSparseMatrixList(List<MatrixType>& MatrixList, word folder,
                          word MatrixName);
 
//--------------------------------------------------------------------------
template <typename T, int Nrows, typename IND>
void ReadSparseMatrixList(List<Eigen::SparseMatrix<T, Nrows, IND >> & m,
                          word folder, word MatrixName);
 
//--------------------------------------------------------------------------
template <typename MatrixType>
void SaveDenseMatrixList(List<MatrixType>& MatrixList, word folder,
                         word MatrixName);
 
 
//--------------------------------------------------------------------------
template <typename MatrixType>
void ReadDenseMatrixList(List<MatrixType>& m, word folder,
                         word MatrixName);
 
//--------------------------------------------------------------------------
int numberOfFiles(word folder, word MatrixName, word ext = "");
 
//--------------------------------------------------------------------------
void writePoints(pointField points, fileName folder, fileName subfolder);
 
//--------------------------------------------------------------------------
void printProgress(double percentage);
 
template<typename T>
void save(const List<Eigen::SparseMatrix<T >> & MatrixList, word folder,
          word MatrixName);
 
template<typename T>
void load(List<Eigen::SparseMatrix<T >> & MatrixList, word folder,
          word MatrixName);
 
//--------------------------------------------------------------------------
template<class Type, template<class> class PatchField, class GeoMesh >
GeometricField<Type, PatchField, GeoMesh> readFieldByIndex(
    const GeometricField<Type, PatchField, GeoMesh>& field,
    fileName casename,
    label index
);
 
};
 
namespace Foam
{
template<typename T>
Ostream& operator<< (Ostream& os, const Eigen::Matrix < T, -1, -1 > & mat)
{
    os << mat.rows() << mat.cols() << UList<T>(const_cast<Eigen::MatrixXd&>
            (mat).data(), mat.size());
    return os;
}
 
template<typename T>
Istream& operator>> (Istream& is, Eigen::Matrix < T, -1, -1 > & mat)
{
    int nrow, ncol;
    is >> nrow >> ncol;
    mat.resize(nrow, ncol);
    UList<T> list(mat.data(), nrow * ncol);
    is >> list;
    return is;
}
 
template<typename T>
Ostream& operator<< (Ostream& os, const Eigen::Matrix < T, -1, 1 > & mat)
{
    os << mat.rows() << mat.cols() << UList<T>(const_cast<Eigen::VectorXd&>
            (mat).data(), mat.size());
    return os;
}
 
template<typename T>
Istream& operator>> (Istream& is, Eigen::Matrix < T, -1, 1 > & mat)
{
    int nrow, ncol;
    is >> nrow >> ncol;
    mat.resize(nrow, ncol);
    UList<T> list(mat.data(), nrow * ncol);
    is >> list;
    return is;
}
 
template<typename T>
Ostream& operator<< (Ostream& os, const Eigen::Tensor<T, 3>& tens)
{
    os << tens.dimension(0) << tens.dimension(1) << tens.dimension(
           2) << UList<T>(const_cast<Eigen::Tensor<T, 3>&>(tens).data(),
                          tens.size());
    return os;
}
 
template<typename T>
Istream& operator>> (Istream& is, Eigen::Tensor<T, 3>& tens)
{
    int d1, d2, d3;
    is >> d1 >> d2 >> d3;
    tens.resize(d1, d2, d3);
    UList<T> list(tens.data(), d1 * d2 * d3);
    is >> list;
    return is;
}
 
}
 
template <typename T, int whatever, typename IND>
void ITHACAstream::SaveSparseMatrix(Eigen::SparseMatrix<T, whatever, IND>& m,
                                    word folder, word MatrixName)
{
    typedef Eigen::Triplet<int> Trip;
    std::vector<Trip> res;
    m.makeCompressed();
    mkDir(folder);
    std::fstream writeFile;
    writeFile.open(folder + MatrixName, std::ios::binary | std::ios::out);
 
    if (writeFile.is_open())
    {
        IND rows, cols, nnzs, outS, innS;
        rows = m.rows()     ;
        cols = m.cols()     ;
        nnzs = m.nonZeros() ;
        outS = m.outerSize();
        innS = m.innerSize();
        writeFile.write(reinterpret_cast<const char*> (& rows), sizeof(IND));
        writeFile.write(reinterpret_cast<const char*> (& cols), sizeof(IND));
        writeFile.write(reinterpret_cast<const char*> (& nnzs), sizeof(IND));
        writeFile.write(reinterpret_cast<const char*> (& outS), sizeof(IND));
        writeFile.write(reinterpret_cast<const char*> (& innS), sizeof(IND));
        writeFile.write(reinterpret_cast<const char*>(m.valuePtr()),
                        sizeof(T  ) * m.nonZeros());
        writeFile.write(reinterpret_cast<const char*>(m.outerIndexPtr()),
                        sizeof(IND) * m.outerSize());
        writeFile.write(reinterpret_cast<const char*>(m.innerIndexPtr()),
                        sizeof(IND) * m.nonZeros());
        writeFile.close();
    }
}
 
template <typename T, int whatever, typename IND>
void ITHACAstream::ReadSparseMatrix(Eigen::SparseMatrix<T, whatever, IND>& m,
                                    word folder, word MatrixName)
{
    std::fstream readFile;
    readFile.open(folder + MatrixName, std::ios::binary | std::ios::in);
 
    if (!readFile.good())
    {
        std::cout << folder + MatrixName <<
        " file does not exist, try to rerun the Offline Stage!" << std::endl;
        exit(EXIT_FAILURE);
    }
    else if (readFile.is_open())
    {
        IND rows, cols, nnz, inSz, outSz;
        readFile.read(reinterpret_cast<char*>(& rows ), sizeof(IND));
        readFile.read(reinterpret_cast<char*>(& cols ), sizeof(IND));
        readFile.read(reinterpret_cast<char*>(& nnz  ), sizeof(IND));
        readFile.read(reinterpret_cast<char*>(& inSz ), sizeof(IND));
        readFile.read(reinterpret_cast<char*>(& outSz), sizeof(IND));
        m.resize(rows, cols);
        m.makeCompressed();
        m.resizeNonZeros(nnz);
        readFile.read(reinterpret_cast<char*>(m.valuePtr()), sizeof(T  ) * nnz  );
        readFile.read(reinterpret_cast<char*>(m.outerIndexPtr()), sizeof(IND) * outSz);
        readFile.read(reinterpret_cast<char*>(m.innerIndexPtr()), sizeof(IND) * nnz );
        m.finalize();
        readFile.close();
    }
}
 
template <typename MatrixType>
void ITHACAstream::SaveDenseMatrix(MatrixType& Matrix, word folder,
                                   word MatrixName)
{
    mkDir(folder);
    std::ofstream out(folder + MatrixName,
                      std::ios::out | std::ios::binary | std::ios::trunc);
    typename MatrixType::Index rows = Matrix.rows(), cols = Matrix.cols();
    out.write(reinterpret_cast<char*> (& rows),
              sizeof(typename MatrixType::Index));
    out.write(reinterpret_cast<char*> (& cols),
              sizeof(typename MatrixType::Index));
    out.write(reinterpret_cast<char*> (Matrix.data()),
              rows * cols * sizeof(typename MatrixType::Scalar) );
    out.close();
}
 
template <typename MatrixType>
void ITHACAstream::ReadDenseMatrix(MatrixType& Matrix, word folder,
                                   word MatrixName)
{
    std::ifstream in;
    in.open(folder + MatrixName, std::ios::in | std::ios::binary);
 
    if (!in.good())
    {
        std::cout << folder + MatrixName <<
        " file does not exist, try to rerun the Offline Stage!" << std::endl;
        exit(EXIT_FAILURE);
    }
    else if (in.is_open())
    {
        typename MatrixType::Index rows = 0, cols = 0;
        in.read(reinterpret_cast<char*> (& rows), sizeof(typename MatrixType::Index));
        in.read(reinterpret_cast<char*> (& cols), sizeof(typename MatrixType::Index));
        Matrix.resize(rows, cols);
        in.read( reinterpret_cast<char*>(Matrix.data()),
                 rows * cols * sizeof(typename MatrixType::Scalar) );
        in.close();
    }
}
 
 
template <typename TensorType>
void ITHACAstream::SaveDenseTensor(TensorType& Tensor, word folder,
                                   word MatrixName)
{
    std::ofstream out(folder + MatrixName,
                      std::ios::out | std::ios::binary | std::ios::trunc);
    typename TensorType::Dimensions dim = Tensor.dimensions();
    int tot = 1;
 
    for (int k = 0; k < dim.size(); k++)
    {
        tot *= dim[k];
    }
 
    out.write(reinterpret_cast<char*> (& dim),
              sizeof(typename TensorType::Dimensions));
    out.write(reinterpret_cast<char*> (Tensor.data()),
              tot * sizeof(typename TensorType::Scalar) );
    out.close();
}
 
template <typename TensorType>
void ITHACAstream::ReadDenseTensor(TensorType& Tensor, word folder,
                                   word MatrixName)
{
    std::ifstream in;
    in.open(folder + MatrixName, std::ios::in | std::ios::binary);
    typename TensorType::Dimensions dim;
    in.read(reinterpret_cast<char*> (& dim),
            sizeof(typename TensorType::Dimensions));
    auto dims = Tensor.dimensions();
    M_Assert(dims.size() == dim.size(),
             "The rank of the tensor you want to fill does not coincide with the rank of the tensor you are reading");
    int tot = 1;
 
    for (int k = 0; k < dim.size(); k++)
    {
        tot *= dim[k];
    }
 
    Tensor.resize(dim);
    in.read( reinterpret_cast<char*>(Tensor.data()),
             tot * sizeof(typename TensorType::Scalar) );
    in.close();
}
 
template <typename MatrixType>
void ITHACAstream::SaveSparseMatrixList(List<MatrixType>& MatrixList,
                                        word folder, word MatrixName)
{
    for (int i = 0; i < MatrixList.size(); i++)
    {
        SaveSparseMatrix(MatrixList[i], folder, MatrixName + name(i));
    }
}
 
template <typename MatrixType>
void ITHACAstream::SaveDenseMatrixList(List<MatrixType>& MatrixList,
                                       word folder, word MatrixName)
{
    for (int i = 0; i < MatrixList.size(); i++)
    {
        SaveDenseMatrix(MatrixList[i], folder, MatrixName + name(i));
    }
}
 
template <typename T, int whatever, typename IND>
void ITHACAstream::ReadSparseMatrixList(
    List<Eigen::SparseMatrix<T, whatever, IND >> & m, word folder, word MatrixName)
{
    int number_of_files = numberOfFiles(folder, MatrixName);
    std::cout << "Reading the Matrix " + folder + MatrixName << std::endl;
    M_Assert(number_of_files != 0,
             "Check if the file you are trying to read exists" );
    Eigen::SparseMatrix<T, whatever, IND> A;
 
    for (int i = 0; i < number_of_files; i++)
    {
        ReadSparseMatrix(A, folder, MatrixName + name(i));
        m.append(A);
    }
}
 
 
template <typename MatrixType>
void ITHACAstream::ReadDenseMatrixList(List<MatrixType>& m, word folder,
                                       word MatrixName)
{
    int number_of_files = numberOfFiles(folder, MatrixName);
    std::cout << "Reading the Matrix " + folder + MatrixName << std::endl;
    M_Assert(number_of_files != 0,
             "Check if the file you are trying to read exists" );
    MatrixType A;
 
    for (int i = 0; i < number_of_files; i++)
    {
        ReadDenseMatrix(A, folder, MatrixName + name(i));
        m.append(A);
    }
}
 
 
 
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
 
#endif