VtkWriter.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS 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 3 of the License, or
   (at your option) any later version.
 
   CEPS 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 CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "VtkWriter.hpp"
 
#include <vtkDataObjectWriter.h>
#include <vtkDataSetToDataObjectFilter.h>
 
#include "common/core/CepsTypeTools.hpp"
#include "pde/problem/AbstractPdeProblem.hpp"
 
VtkWriter::VtkWriter(AbstractDiscretization* discr,
                     const CepsString&       baseName,
                     CepsUInt                nOutputs,
                     CepsOutputFormat        outputType,
                     CepsBool                writeGeomIndices) :
  m_discr           (discr),
  m_writtenTimes    ({}),
  m_internalIndex   (0),
  m_outputType      (CepsOutputFormat::ParaviewSeries),
  m_updateMesh      (true),
  m_writeGeomIndices(writeGeomIndices),
  m_unstructuredGrid()
{
  
  setOutputType(outputType);
 
  CEPS_ABORT_IF(ceps::isNullPtr(discr),
    "Vtk writer: nullptr pointer to discretization passed"
  );
  m_geom = discr->getProblem()->getGeometry();
 
  // Creating a subdirectory where to put vtu and pvtu files
  m_outputDir   = ceps::getDir(baseName);
  m_baseName    = ceps::getBaseName(baseName);
  m_outputNb    = 0;
  m_outputLgMax = static_cast<CepsUInt>(floor(log10(nOutputs))+1);
  m_writtenTimes.reserve(nOutputs+1);
 
  CEPS_ABORT_IF(not ceps::isPathInOneWord(m_outputDir),
    "Provided output directory is not a valid path name (do not forget \\ before spaces).\n"
    << "Given path: \n  " << m_outputDir
  );
 
  m_unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
 
  // Tools to handle the writing of pvtu file on multiple processors
  m_vtkMpiComm    = vtkSmartPointer<vtkMPICommunicator>::New();
  m_opaqueCommVTK = ceps::getPtrCommunicator();
  m_vtkMpiComm->InitializeExternal(&m_opaqueCommVTK);
 
  m_vtkMpiController = vtkSmartPointer<vtkMPIController>::New();
  m_vtkMpiController->SetCommunicator(m_vtkMpiComm);
 
  if (ceps::isMaster())
  {
    CepsString command = "mkdir -p " + m_outputDir + "/snapshots";
    ceps::execute(command);
  }
}
 
VtkWriter::~VtkWriter()
{
}
 
void
VtkWriter::addScalarData(const DHVecPtr& v, const CepsString &fieldName, Unknown* u)
{
  addScalarData(v,CepsVector<CepsString>({fieldName}),
                  CepsVector<Unknown*>({u}));
}
 
void
VtkWriter::addScalarData(const DHVecPtr& solution, const CepsVector<CepsString>& fieldNames,
                         const CepsVector<Unknown*>& us)
{
  if (ceps::isProfiling())
    getProfiler()->start("write","writing solution");
 
  // check if the mesh must be updated
  if (m_updateMesh)
    updateMeshData();
 
  CEPS_ABORT_IF(fieldNames.size()!=us.size(),
    "Mismatch between number of field names and number of unknowns given to wtk writer"
  );
 
  for (CepsUInt iu = 0u; iu<us.size(); iu++)
  {
 
    Unknown* u = us[iu];
 
    // Few sanity checks
    CEPS_ABORT_IF(ceps::isNullPtr(u),
      "Pointer to unknown passed to addScalar is null"
    );
 
    // Get data for this unknown
    CepsVector<CepsReal > data       ;
    CepsVector<CepsIndex> dataIndices;
    CepsBool              reduceOnMaster = isEnabledOption(m_outputType,CepsOutputFormat::VTKLegacy) 
                                        or isEnabledOption(m_outputType,CepsOutputFormat::VTKCeps);
    m_discr->extractValuesForUnknown(
      u->getIdentifier(),solution,dataIndices,data,
      reduceOnMaster,u->getLocation()!=CepsLocationFlag::ZeroD
    );
 
    if (ceps::isMaster() or not reduceOnMaster) // For those who write
    {
      CepsLocationFlag loc      = u->getLocation();
      CepsUInt         nToParse = data.size();
 
      // There may be additional nodes (P2) ! 
      CepsUInt nToWrite =  loc == CepsLocationFlag::ZeroD ? 1
                        :( loc == CepsLocationFlag::Point ? m_geom->getNbNodes() 
                        :                                   m_geom->getNbCells());
 
      vtkDoubleArray* scalars = vtkDoubleArray::New();
      scalars->SetName(fieldNames[iu].c_str());
      scalars->SetNumberOfTuples(nToWrite);
      scalars->Fill(NAN);
      for (CepsUInt i=0; i<nToParse; i++)
      {
        CepsIndex vtkId =  loc == CepsLocationFlag::ZeroD ? 0 
                        :( loc == CepsLocationFlag::Point ? m_geomToVtkNodeIndex.at(dataIndices[i])
                        :                                   m_geomToVtkCellIndex.at(dataIndices[i]));
        if (vtkId<CepsIndex(nToWrite))
          scalars->SetTuple1(vtkId,data[i]);
      }
 
      // this replaces any existing scalar field with the same name
      if (loc==CepsLocationFlag::ZeroD)
        m_unstructuredGrid->GetFieldData()->AddArray(scalars);
      else if (loc == CepsLocationFlag::Point)
        m_unstructuredGrid->GetPointData()->AddArray(scalars);
      else
        m_unstructuredGrid->GetCellData()->AddArray(scalars);
      scalars->Delete();
 
    }
 
  }
 
  if (ceps::isProfiling ())
    getProfiler ()->stop ("write");
}
 
void
VtkWriter::write(CepsReal time, CepsBool writeXmlEntry)
{
  // check if the mesh must be updated
  CEPS_ABORT_IF(m_updateMesh,
    "No data has been given to the writer with addScalarData"
  );
 
  if (ceps::isProfiling())
    getProfiler()->start("write");
 
  m_writtenTimes.push_back(time);
  addTimeFieldToUGrid(time);
 
  // We suppose data has been regrouped and given to the master node.
  if (isEnabledOption(m_outputType, CepsOutputFormat::Music))
  {
    if (ceps::isMaster())
    {
      writeLegacyVTK();
      writeMetaDataSet();
    }
  }
  else if (isEnabledOption(m_outputType, CepsOutputFormat::ParaviewSeries))
  {
    writeVTU();
    if (ceps::isMaster())
      writeSeriesFile();
  }
  else if (isEnabledOption(m_outputType, CepsOutputFormat::VTKCeps))
  {
    if (ceps::isMaster())
      writeDataArray();
  }
 
  if (ceps::isProfiling())
    getProfiler()->stop("write");
 
  m_outputNb++;
}
 
CepsString
VtkWriter::getByeByeMessage ()
{
  CepsString ext = "";
 
  return "Results can be viewed by opening "
          + m_outputDir + "/" + m_baseName
          + getOuputExtension();
}
 
void
VtkWriter::setOutputType(CepsOutputFormat flag)
{
  m_outputType = flag;
}
 
CepsOutputFormat
VtkWriter::getOutputType() const
{
  return m_outputType;
}
 
CepsString 
VtkWriter::getOuputExtension() const
{
  if (isEnabledOption(m_outputType, CepsOutputFormat::Music))
    return ".xml";
  else if (isEnabledOption(m_outputType, CepsOutputFormat::ParaviewSeries))
    return ".pvtu.series";  
  else if (isEnabledOption(m_outputType, CepsOutputFormat::VTKCeps))
    return".vtk_ceps";
  return "";
}
 
// -----------------------------------------------------------------------------------------------====
// Protected
// -----------------------------------------------------------------------------------------------====
 
void  
VtkWriter::writeVTU()
{
  std::stringstream filename;
  filename << m_outputDir << "/snapshots/" << m_baseName << "-" << std::setfill('0') 
           << std::setw (m_outputLgMax) << m_outputNb << ".pvtu";
 
  auto writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
  writer->SetDataModeToBinary();
  writer->SetNumberOfPieces(ceps::getGridSize());
  writer->SetStartPiece(ceps::getRank());
  writer->SetEndPiece(ceps::getRank());
  writer->SetCompressor(nullptr);
  writer->SetInputData(m_unstructuredGrid);
  writer->SetController(m_vtkMpiController);
  writer->SetFileName(filename.str().c_str());
  writer->Update();
  writer->Write();
}
 
void
VtkWriter::writeLegacyVTK()
{
  std::stringstream filename;
 
  filename << m_outputDir << "/snapshots/" << m_baseName << "-" << std::setfill('0') 
            << std::setw(m_outputLgMax) << m_outputNb << ".vtk";
 
  auto writer = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
  writer->SetFileTypeToBinary();
  writer->SetInputData(m_unstructuredGrid);
  writer->SetFileName(filename.str().c_str());
  writer->Write();
}
 
void
VtkWriter::writeDataArray()
{
  std::stringstream filename;
  
  filename << m_outputDir << "/snapshots/" << m_baseName << "-" << std::setfill('0')
           << std::setw(m_outputLgMax) << m_outputNb << ".vtk_ceps";
 
  auto fields = vtkSmartPointer<vtkDataSetToDataObjectFilter>::New();
  fields->AddInputData(m_unstructuredGrid);
  fields->CellDataOn();
  fields->PointDataOn();
  fields->FieldDataOn();
  fields->GeometryOff();
  fields->TopologyOff();
 
  // fields->AddInputDataObject(m_unstructuredGrid);
  auto writer = vtkSmartPointer<vtkDataObjectWriter>::New();
  writer->SetFileTypeToBinary();
  writer->SetFieldDataName("data");
  writer->AddInputConnection(fields->GetOutputPort());
  writer->SetFileName(filename.str().c_str());
  writer->Write();
}
 
void 
VtkWriter::writeSeriesFile()
{
  CepsString series = m_outputDir + "/" + m_baseName + ".pvtu.series";
  std::ofstream file(series, std::ofstream::trunc);
  std::stringstream filename;
 
  file << "{" << std::endl;
  file << "  \"file-series-version\" : \"1.0\"," << std::endl;
  file << "  \"files\" : [" << std::endl;
 
  for (CepsUInt i = 0UL; i < m_writtenTimes.size(); ++i)
  {
    filename << "./snapshots/" << m_baseName << "-" << std::setfill('0') 
             << std::setw (m_outputLgMax) << i << ".pvtu";
 
    file << "    { ";
    file << "\"time\" : " << m_writtenTimes[i] << ", ";
    file << "\"name\" : \"" << filename.str() << "\", ";
    file << "}," << std::endl;
 
    filename.str("");
  }
  
  file << "  ]" << std::endl;
  file << "}" << std::endl;
  file.close();
}
 
void 
VtkWriter::writeMetaDataSet()
{
  CepsString meta = m_outputDir + "/" + m_baseName + ".xml";
  std::ofstream file(meta, std::ofstream::trunc);
  std::stringstream filename;
 
  file << "<?xml version=\"1.0\"?>" << endl;
  file << "<VTKFile type=\"vtkDataManager\" version=\"0.1\" byte_order=\"LittleEndian\" "
          "compressor=\"vtkZLibDataCompressor\">" << std::endl;
  file << "  <vtkDataManager>" << std::endl;
 
  for (CepsUInt i = 0UL; i < m_writtenTimes.size(); ++i)
  {
    filename << "./snapshots/" << m_baseName << "-" << std::setfill('0') 
             << std::setw (m_outputLgMax) << i << ".vtk";
 
    file << "    <vtkMetaDataSet type=\"1\" name=\"\" tag=\"\">" << std::endl;
    file << "      <metadata key=\"Time\" value=\"" << m_writtenTimes[i] << "\">" << std::endl;
    file << "      <DataSet group=\"0\" dataset=\"" << i << "\" file=\"" << filename.str() << "\"/>" << std::endl;
    file << "    </vtkMetaDataSet>" << std::endl;
 
    filename.str("");
  }
 
  file << "   </vtkDataManager>" << std::endl;
  file << "</VTKFile>" << std::endl;
  file.close();
}
 
void
VtkWriter::updateMeshData()
{
 
  // Reset if needed
  m_internalIndex = 0;
  m_geomToVtkNodeIndex.clear();
 
  insertPointData();
  insertElementData();
 
  // Geometry data is up to date
  m_updateMesh = false;
}
 
void
VtkWriter::insertPointData()
{
 
  // Vtk array for points coordinates
  vtkPoints* points = vtkPoints::New(VTK_DOUBLE);
  points->GetData()->SetName("Vertex coordinates");
 
  // Number of nodes (all of them if writing legacy vtk)
  CepsBool reduceOnMaster = isEnabledOption(m_outputType, CepsOutputFormat::Music) 
                            or isEnabledOption(m_outputType, CepsOutputFormat::VTKCeps);
  CepsUInt nToWrite = reduceOnMaster ? m_geom->getNbNodes() : m_geom->getNbNodesLocal();
 
  // Vtk array for global indices (after partitioning)
  vtkIdTypeArray* globalIndices = vtkIdTypeArray::New();
  if (m_writeGeomIndices)
  {
    globalIndices->SetName("globalIndices");
    globalIndices->SetNumberOfComponents(1);
    globalIndices->SetNumberOfTuples(nToWrite);
    globalIndices->FillComponent(0,-1);
  }
 
  // Now fill the arrays
  if (reduceOnMaster)
  {
    // Write whole mesh in one fle : we follow the original ordering of the geometry,
    // before partitioning
    CepsVector<CepsIndex> partMap = m_geom->getNodePartitionMap();
 
    for (CepsUInt i=0; i<nToWrite; i++)
    {
      CepsIndex gId = partMap[i];
      // We build the reverse node partition mapping (for cells later on) 
      if (ceps::isMaster())
      {
        m_geomToVtkNodeIndex[gId] = i;
        if (m_writeGeomIndices)
          globalIndices->SetTuple1(i,gId);
      }
 
      GeomNode* n = m_geom->getNode(gId,true); // Fetch only owned nodes
      CepsVector<CepsReal> xyz(3,0.e0);
 
      if (ceps::isParallel())
      {
        // Send coordinates to everyone (incl master)
        CepsVector<CepsReal> myxyz(3,0.e0);
        if (ceps::isValidPtr(n))
          for (CepsUInt j=0;j<3;j++)
            myxyz[j] = n->getCoordinate(j);
        MPI_Allreduce(myxyz.data(),xyz.data(),3,CEPS_MPI_REAL,MPI_SUM,ceps::getCommunicator());
      }
      else
        // Just get coordinates
        for (CepsUInt j=0;j<3;j++)
          xyz[j] = n->getCoordinate(j);
 
      if (ceps::isMaster())
        points->InsertNextPoint(xyz[0],xyz[1],xyz[2]);
    }
  }
  else
  {
    // Write blocks.
    // In each block, we follow the same order for nodes:
    // 3D owned, 3D halo, 2D owned, 2D halo, 1D owned, 1D halo
    // We build a map between geomID and vtkIndex in this block along the way
    for (CepsUInt dim=3;dim>0;dim--)
      if (m_geom->hasMeshOfDim(dim))
      {
        for (const CepsVector<GeomNode*>& nodeVector: {m_geom->getMeshOfDim(dim)->getOwnedNodes(),
                                                       m_geom->getMeshOfDim(dim)->getHaloNodes ()})
          for (GeomNode* n : nodeVector)
          {
            m_geomToVtkNodeIndex[n->getGlobalIndex()] = m_internalIndex;
            points->InsertNextPoint(n->x(),n->y(),n->z());
            if (m_writeGeomIndices)
              globalIndices->SetTuple1(m_internalIndex,n->getGlobalIndex());
            m_internalIndex++;
          }
      }
  }
 
  m_unstructuredGrid->SetPoints(points);
  if (m_writeGeomIndices)
    m_unstructuredGrid->GetPointData()->AddArray(globalIndices);
 
  points       ->Delete();
  globalIndices->Delete();
 
}
 
void
VtkWriter::insertElementData()
{
  CepsBool reduceOnMaster = isEnabledOption(m_outputType, CepsOutputFormat::Music) 
    or isEnabledOption(m_outputType, CepsOutputFormat::VTKCeps);
 
  if (reduceOnMaster)
  {
    // We write the whole mesh, so we follow the original ordering of cells
    CepsUInt              nToWrite = m_geom->getNbCells();
    CepsVector<CepsIndex> partMap  = m_geom->getCellPartitionMap();
 
    for (CepsUInt i=0; i<nToWrite; i++)
    {
      CepsIndex gId = partMap[i];
      if (ceps::isMaster())
        m_geomToVtkCellIndex[gId] = i;
 
      GeomCell* c   = m_geom->getCell(gId);
      // Get cell dimension
      CepsUInt  dim;
      if (ceps::isParallel())
      {
        CepsUInt mydim = ceps::isValidPtr(c) ? c->getDimension() : 0u;
        MPI_Allreduce(&mydim,&dim,1,CEPS_MPI_UINT,MPI_SUM,ceps::getCommunicator());
      }
      else
        dim = c->getDimension();
 
      auto vtkCell = allocateCellOfDim(dim);
 
      // Get node IDs
      vtkIdList *nodeIds = vtkCell->GetPointIds();
      for (vtkIdType j=0; j<vtkCell->GetNumberOfPoints(); j++)
      {
        CepsIndex nodeId;
        if (ceps::isParallel())
        {
          CepsIndex mynodeId = ceps::isValidPtr(c) ? c->getNodeAt(j)->getGlobalIndex() : 0u;
          MPI_Allreduce(&mynodeId,&nodeId,1,CEPS_MPI_UINT,MPI_SUM,ceps::getCommunicator());
        }
        else
          nodeId = c->getNodeAt(j)->getGlobalIndex();
 
        nodeId = m_geomToVtkNodeIndex[nodeId];
        nodeIds->SetId(j,nodeId);
      }
 
      m_unstructuredGrid->InsertNextCell(vtkCell->GetCellType(),nodeIds);
      vtkCell->Delete();
 
    }
 
  }
  else
  {
    // We write blocks, so we follow the order 3D, 2D, 1D of cells in the block
    m_internalIndex = 0;
    if (m_geom->has3dMesh())
      insertElementsOfDim(3);
    if (m_geom->has2dMesh())
      insertElementsOfDim(2);
    if (m_geom->has1dMesh())
      insertElementsOfDim(1);
  }
}
 
void
VtkWriter::insertElementsOfDim(CepsUInt dim)
{
  Mesh* m = m_geom->getMeshOfDim(dim);
  for (GeomCell* cell : m->getCells())
    insertCell(cell);
  for (GeomCell* cell : m->getBoundaryCells())
    insertCell(cell);
}
 
void
VtkWriter::insertCell(GeomCell* cell)
{
  auto       vtkCell = allocateCellOfDim(cell->getDimension());
  vtkIdList *nodeIds = vtkCell->GetPointIds();
 
  for (vtkIdType j=0; j<vtkCell->GetNumberOfPoints(); j++)
  {
    CepsIndex nodeId = cell->getNodeAt(j)->getGlobalIndex();
    nodeId = m_geomToVtkNodeIndex[nodeId];
    nodeIds->SetId(j,nodeId);
  }
 
  m_geomToVtkCellIndex[cell->getGlobalIndex()] = m_internalIndex;
  m_internalIndex++;
 
  m_unstructuredGrid->InsertNextCell(vtkCell->GetCellType(),nodeIds);
  vtkCell->Delete();
}
 
vtkCell*
VtkWriter::allocateCellOfDim(CepsUInt dim)
{
  if (dim== 3U) 
    return vtkTetra::New();
  else if (dim== 2U)
    return vtkTriangle::New();
  else if (dim== 1U) 
    return vtkLine::New();
 
  return vtkVertex::New();
}
 
void
VtkWriter::addTimeFieldToUGrid(CepsReal time)
{
  auto timef = vtkSmartPointer<vtkDoubleArray>::New();
  timef->SetNumberOfValues(1);
  timef->SetValue(0,time);
  timef->SetName("TIME_STEP");
  m_unstructuredGrid->GetFieldData()->AddArray(timef);
}