vtkmultiwriter.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright (C) 2008-2013 by Andreas Lauser
  Copyright (C) 2009 by Anneli Schoeniger

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM 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 General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef EWOMS_VTK_MULTI_WRITER_HH
#define EWOMS_VTK_MULTI_WRITER_HH

#include "vtkscalarfunction.hh"
#include "vtkvectorfunction.hh"
#include "vtktensorfunction.hh"

#include <ewoms/io/baseoutputwriter.hh>

#include <opm/material/common/Valgrind.hpp>

#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

#if HAVE_MPI
#include <mpi.h>
#endif

#include <list>
#include <string>
#include <limits>
#include <sstream>
#include <fstream>

namespace Ewoms {
template <class GridView, int vtkFormat>
class VtkMultiWriter : public BaseOutputWriter
{
    enum { dim = GridView::dimension };

    typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGVertexLayout> VertexMapper;
    typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;

public:
    typedef BaseOutputWriter::Scalar Scalar;
    typedef BaseOutputWriter::Vector Vector;
    typedef BaseOutputWriter::Tensor Tensor;
    typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
    typedef BaseOutputWriter::VectorBuffer VectorBuffer;
    typedef BaseOutputWriter::TensorBuffer TensorBuffer;

    typedef Dune::VTKWriter<GridView> VtkWriter;

    VtkMultiWriter(const GridView &gridView,
                   const std::string &simName = "",
                   std::string multiFileName = "")
        : gridView_(gridView)
        , elementMapper_(gridView)
        , vertexMapper_(gridView)
    {
        simName_ = (simName.empty()) ? "sim" : simName;
        multiFileName_ = multiFileName;
        if (multiFileName_.empty()) {
            multiFileName_ = simName_;
            multiFileName_ += ".pvd";
        }

        curWriterNum_ = 0;

        commRank_ = gridView.comm().rank();
        commSize_ = gridView.comm().size();
    }

    ~VtkMultiWriter()
    {
        finishMultiFile_();

        if (commRank_ == 0)
            multiFile_.close();
    }

    int curWriterNum() const
    { return curWriterNum_; }

    void gridChanged()
    {
        elementMapper_.update();
        vertexMapper_.update();
    }

    void beginWrite(double t)
    {
        if (!multiFile_.is_open()) {
            startMultiFile_(multiFileName_);
        }

        curWriter_ = new VtkWriter(gridView_, Dune::VTK::conforming);
        ++curWriterNum_;

        curTime_ = t;
        curOutFileName_ = fileName_();
    }

    ScalarBuffer *allocateManagedScalarBuffer(int numEntities)
    {
        ScalarBuffer *buf = new ScalarBuffer(numEntities);
        managedScalarBuffers_.push_back(buf);
        return buf;
    }

    VectorBuffer *allocateManagedVectorBuffer(int numOuter, int numInner)
    {
        VectorBuffer *buf = new VectorBuffer(numOuter);
        for (int i = 0; i < numOuter; ++ i)
            (*buf)[i].resize(numInner);

        managedVectorBuffers_.push_back(buf);
        return buf;
    }

    void attachScalarVertexData(ScalarBuffer &buf, std::string name)
    {
        sanitizeScalarBuffer_(buf);

        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkScalarFunction<GridView, VertexMapper> VtkFn;
        FunctionPtr fnPtr(new VtkFn(name,
                                    gridView_,
                                    vertexMapper_,
                                    buf,
                                    /*codim=*/dim));
        curWriter_->addVertexData(fnPtr);
    }

    void attachScalarElementData(ScalarBuffer &buf, std::string name)
    {
        sanitizeScalarBuffer_(buf);

        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkScalarFunction<GridView, ElementMapper> VtkFn;
        FunctionPtr fnPtr(new VtkFn(name,
                                    gridView_,
                                    elementMapper_,
                                    buf,
                                    /*codim=*/0));
        curWriter_->addCellData(fnPtr);
    }

    void attachVectorVertexData(VectorBuffer &buf, std::string name)
    {
        sanitizeVectorBuffer_(buf);

        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkVectorFunction<GridView, VertexMapper> VtkFn;
        FunctionPtr fnPtr(new VtkFn(name,
                                    gridView_,
                                    vertexMapper_,
                                    buf,
                                    /*codim=*/dim));
        curWriter_->addVertexData(fnPtr);
    }

    void attachTensorVertexData(TensorBuffer &buf, std::string name)
    {
        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkTensorFunction<GridView, VertexMapper> VtkFn;

        for (size_t colIdx = 0; colIdx < buf[0].N(); ++colIdx) {
            std::ostringstream oss;
            oss << name <<  "[" << colIdx << "]";

            FunctionPtr fnPtr(new VtkFn(oss.str(),
                                        gridView_,
                                        vertexMapper_,
                                        buf,
                                        /*codim=*/dim,
                                        colIdx));
            curWriter_->addVertexData(fnPtr);
        }
    }

    void attachVectorElementData(VectorBuffer &buf, std::string name)
    {
        sanitizeVectorBuffer_(buf);

        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkVectorFunction<GridView, ElementMapper> VtkFn;
        FunctionPtr fnPtr(new VtkFn(name,
                                    gridView_,
                                    elementMapper_,
                                    buf,
                                    /*codim=*/0));
        curWriter_->addCellData(fnPtr);
    }

    void attachTensorElementData(TensorBuffer &buf, std::string name)
    {
        typedef typename VtkWriter::VTKFunctionPtr FunctionPtr;
        typedef Ewoms::VtkTensorFunction<GridView, ElementMapper> VtkFn;

        for (size_t colIdx = 0; colIdx < buf[0].N(); ++colIdx) {
            std::ostringstream oss;
            oss << name <<  "[" << colIdx << "]";

            FunctionPtr fnPtr(new VtkFn(oss.str(),
                                        gridView_,
                                        elementMapper_,
                                        buf,
                                        /*codim=*/0,
                                        colIdx));
            curWriter_->addCellData(fnPtr);
        }
    }

    void endWrite(bool onlyDiscard = false)
    {
        if (!onlyDiscard) {
            std::string fileName;
            // write the actual data as vtu or vtp (plus the pieces file in the parallel case)
            fileName = curWriter_->write(/*name=*/curOutFileName_.c_str(),
                                         static_cast<Dune::VTK::OutputType>(vtkFormat));

            // determine name to write into the multi-file for the
            // current time step
            multiFile_.precision(16);
            multiFile_ << "   <DataSet timestep=\"" << curTime_ << "\" file=\""
                       << fileName << "\"/>\n";
        }
        else
            --curWriterNum_;

        // discard managed objects and the current VTK writer
        delete curWriter_;
        while (managedScalarBuffers_.begin() != managedScalarBuffers_.end()) {
            delete managedScalarBuffers_.front();
            managedScalarBuffers_.pop_front();
        }
        while (managedVectorBuffers_.begin() != managedVectorBuffers_.end()) {
            delete managedVectorBuffers_.front();
            managedVectorBuffers_.pop_front();
        }

        // temporarily write the closing XML mumbo-jumbo to the mashup
        // file so that the data set can be loaded even if the
        // simulation is aborted (or not yet finished)
        finishMultiFile_();
    }

    template <class Restarter>
    void serialize(Restarter &res)
    {
        res.serializeSectionBegin("VTKMultiWriter");
        res.serializeStream() << curWriterNum_ << "\n";

        if (commRank_ == 0) {
            size_t fileLen = 0;
            size_t filePos = 0;
            if (multiFile_.is_open()) {
                // write the meta file into the restart file
                filePos = multiFile_.tellp();
                multiFile_.seekp(0, std::ios::end);
                fileLen = multiFile_.tellp();
                multiFile_.seekp(filePos);
            }

            res.serializeStream() << fileLen << "  " << filePos << "\n";

            if (fileLen > 0) {
                std::ifstream multiFileIn(multiFileName_.c_str());
                char *tmp = new char[fileLen];
                multiFileIn.read(tmp, static_cast<long>(fileLen));
                res.serializeStream().write(tmp, fileLen);
                delete[] tmp;
            }
        }

        res.serializeSectionEnd();
    }

    template <class Restarter>
    void deserialize(Restarter &res)
    {
        res.deserializeSectionBegin("VTKMultiWriter");
        res.deserializeStream() >> curWriterNum_;

        if (commRank_ == 0) {
            std::string dummy;
            std::getline(res.deserializeStream(), dummy);

            // recreate the meta file from the restart file
            size_t filePos, fileLen;
            res.deserializeStream() >> fileLen >> filePos;
            std::getline(res.deserializeStream(), dummy);
            if (multiFile_.is_open())
                multiFile_.close();

            if (fileLen > 0) {
                multiFile_.open(multiFileName_.c_str());

                char *tmp = new char[fileLen];
                res.deserializeStream().read(tmp, fileLen);
                multiFile_.write(tmp, fileLen);
                delete[] tmp;
            }

            multiFile_.seekp(filePos);
        }
        else {
            std::string tmp;
            std::getline(res.deserializeStream(), tmp);
        }
        res.deserializeSectionEnd();
    }

private:
    std::string fileName_()
    {
        // use a new file name for each time step
        std::ostringstream oss;
        oss << simName_ << "-" << std::setw(5) << std::setfill('0')
            << curWriterNum_;
        return oss.str();
    }

    std::string fileSuffix_()
    { return (GridView::dimension == 1) ? "vtp" : "vtu"; }

    void startMultiFile_(const std::string &multiFileName)
    {
        // only the first process writes to the multi-file
        if (commRank_ == 0) {
            // generate one meta vtk-file holding the individual time steps
            multiFile_.open(multiFileName.c_str());
            multiFile_ << "<?xml version=\"1.0\"?>\n"
                          "<VTKFile type=\"Collection\"\n"
                          "         version=\"0.1\"\n"
                          "         byte_order=\"LittleEndian\"\n"
                          "         compressor=\"vtkZLibDataCompressor\">\n"
                          " <Collection>\n";
        }
    }

    void finishMultiFile_()
    {
        // only the first process writes to the multi-file
        if (commRank_ == 0) {
            // make sure that we always have a working meta file
            std::ofstream::pos_type pos = multiFile_.tellp();
            multiFile_ << " </Collection>\n"
                          "</VTKFile>\n";
            multiFile_.seekp(pos);
            multiFile_.flush();
        }
    }

    // make sure the field is well defined if running under valgrind
    // and make sure that all values can be displayed by paraview
    void sanitizeScalarBuffer_(ScalarBuffer &b)
    {
        for (unsigned j = 0; j < b.size(); ++j) {
            Valgrind::CheckDefined(b[j]);

            // set values which are too small to 0 to avoid
            // problems with paraview
            if (std::abs(b[j]) < std::numeric_limits<float>::min()) {
                b[j] = 0.0;
            }
        }
    }

    void sanitizeVectorBuffer_(VectorBuffer &b)
    {
        size_t nComps = b[0].size();
        for (size_t i = 0; i < b.size(); ++i) {
            for (size_t j = 0; j < nComps; ++j) {
                Valgrind::CheckDefined(b[i][j]);

                // set values which are too small to 0 to avoid
                // problems with paraview
                if (std::abs(b[i][j]) < std::numeric_limits<float>::min()) {
                    b[i][j] = 0.0;
                }
            }
        }
    }

    const GridView gridView_;
    ElementMapper elementMapper_;
    VertexMapper vertexMapper_;

    std::string simName_;
    std::ofstream multiFile_;
    std::string multiFileName_;

    int commSize_; // number of processes in the communicator
    int commRank_; // rank of the current process in the communicator

    VtkWriter *curWriter_;
    double curTime_;
    std::string curOutFileName_;
    int curWriterNum_;

    std::list<ScalarBuffer *> managedScalarBuffers_;
    std::list<VectorBuffer *> managedVectorBuffers_;
};
} // namespace Ewoms

#endif