VtkMeshReader.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 "geometry/readers/VtkMeshReader.hpp"
 
VtkMeshReader::VtkMeshReader(const CepsString &fileName, const CepsUInt &dim, CepsBool quads) :
    MeshReader (fileName, dim, quads),
    m_mesh(vtkSmartPointer<vtkDataSet>())
{
  CEPS_ABORT_IF(not isVtkFile(m_fileName),
    "Implemented vtk mesh reader is not compatible with extension "
    << ceps::getExtension (m_fileName) << "." << std::endl
    << "   Valid extenstions are vtk and pvtu."
  );
}
 
VtkMeshReader::~VtkMeshReader()
{
}
 
void
VtkMeshReader::initialize()
{
  CEPS_SAYS("      Reading " << m_fileName);
 
  // determine if reading a legacy file or a pvtu file
  CepsString file = ceps::getFilename (m_fileName);
  CepsString dir  = ceps::getDir (m_fileName);
  CepsString ext  = ceps::getExtension (m_fileName);
 
  if (ext=="vtk")
    initializeLegacy(m_fileName);
  else if (ext=="pvtu")
    initializePvtu(m_fileName);
 
  // vertices
  m_nbNodes = m_mesh->GetNumberOfPoints();
 
  m_nbCells     = 0u;
  m_nbBdryCells = 0u;
 
  // try to find cells first
  vtkIdType cellType, nbCellsTotal = m_mesh->GetNumberOfCells();
  m_nbCellsTotal = nbCellsTotal;
  vtkIdType  cellNbPoints;
  vtkIdList *cellPoints = vtkIdList::New();
  for (vtkIdType i = 0; i < nbCellsTotal; ++i)
  {
    cellType = m_mesh->GetCellType(i);
    switch (cellType)
    {
      // case VTK_HEXAHEDRON:
      // case VTK_VOXEL:
      //   if (m_quads and m_expectedDim == 3)
      //     m_nbCells++;
      //   break;
      case VTK_TETRA:
        if (not m_quads and m_expectedDim == 3)
          m_nbCells++;
        break;
      case VTK_TRIANGLE:
      {
        if (not m_quads)
        {
          if (m_expectedDim == 3)
            m_nbBdryCells++;
          else if (m_expectedDim == 2)
            m_nbCells++;
        }
        break;
      }
      // case VTK_QUAD:
      // {
      //   if (m_quads)
      //   {
      //     if (m_expectedDim == 3)
      //       m_nbBdryCells++;
      //     else if (m_expectedDim == 2)
      //       m_nbCells++;
      //   }
      //   break;
      // }
      case VTK_POLYGON:
      {
        m_mesh->GetCellPoints(i, cellPoints);
        cellNbPoints = cellPoints->GetNumberOfIds ();
        if ((cellNbPoints == 3 and not m_quads) or (cellNbPoints == 4 and m_quads))
        {
          if (m_expectedDim == 3)
            m_nbBdryCells++;
          else if (m_expectedDim == 2)
            m_nbCells++;
        }
        break;
      }
      case VTK_LINE:
      {
        if (m_expectedDim == 2)
          m_nbBdryCells++;
        else if (m_expectedDim == 1)
          m_nbCells++;
        break;
      }
    }
  }
 
  CEPS_ABORT_IF(m_nbCells == 0 and m_nbBdryCells == 0,
    "No cells or boundary cells found in file " << m_fileName << " for mesh of dimension "
                                                << m_expectedDim
  );
 
  m_nbAttrPerNode     = 0;
  m_nbAttrPerCell     = 0;
  m_nbAttrPerBdryCell = 0;
  // check for attributes
  determineCellDataToRead ();
  determinePointDataToRead();
 
  m_initialized = true;
}
 
void
VtkMeshReader::readNodes(
  CepsUInt                   indexStart,
  CepsUInt                   indexEnd,
  CepsUInt                   nodeOffset,
  CepsVector<CepsReal>&      coords,
  CepsVector<CepsAttribute>& attr,
  CepsVector<CepsUInt>&      attrPtr
)
{
  checkCanReadNodes(indexStart,indexEnd);
 
  double       pt[3];
  vtkPointSet *p = vtkPointSet::SafeDownCast(m_mesh);
 
  for (vtkIdType i=0; (CepsUInt) i<m_nbNodes; ++i)
  {
    CepsUInt nodeId = nodeOffset+i;
    p->GetPoint (i,pt);
    if (i>=indexStart and i<indexEnd)
    {
      coords[3*nodeId  ] = static_cast<CepsReal>(pt[0]);
      coords[3*nodeId+1] = static_cast<CepsReal>(pt[1]);
      coords[3*nodeId+2] = static_cast<CepsReal>(pt[2]);
 
      // read attributes
      for (CepsUInt j=0; j<m_pointArrayToRead.size(); ++j)
        attr[attrPtr[nodeId]+j] = m_pointArrayToRead[j]->GetTuple1(i);
 
      attrPtr[nodeId+1] = attrPtr[nodeId] + m_pointArrayToRead.size();
    }
  }
 
}
 
void
VtkMeshReader::readCells(
  CepsUInt                         indexStart,
  CepsUInt                         indexEnd,
  CepsUInt                         nodeOffset,
  CepsUInt                         cellOffset,
  CepsVector<CepsNodeGlobalIndex>& nodeIndices,
  CepsVector<CepsUInt>&            nodeIndicesPtr,
  CepsVector<CepsAttribute>&       cellAttr,
  CepsVector<CepsUInt>&            cellAttrPtr,
  CepsVector<CepsChar>&            isBoundaryCell
)
{
  checkCanReadCells(indexStart,indexEnd);
 
  // We do the loop twice, the first parses the volumic cells, the second the boundary cells
  for (vtkIdType i=0; i<m_nbCellsTotal; ++i)
    getNextCell(
      cellOffset+m_cellsRead,i,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,false
    );
  
  for (vtkIdType i=0; i<m_nbCellsTotal; ++i)
    getNextCell(
      cellOffset+m_cellsRead+m_bdryCellsRead,i,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,true
    );
 
  return;
}
 
CepsBool
VtkMeshReader::isVtkFile(const CepsString &fileName) const
{
  CepsString extension = ceps::getExtension (fileName);
  return (extension == "vtk" or extension == "pvtu");
}
 
void
VtkMeshReader::getNextCell(
  CepsUInt                         cellId,
  CepsUInt                         fileCellId,
  CepsUInt                         nodeOffset,
  CepsVector<CepsNodeGlobalIndex>& nodeIndices,
  CepsVector<CepsUInt>&            nodeIndicesPtr,
  CepsVector<CepsAttribute>&       cellAttr,
  CepsVector<CepsUInt>&            cellAttrPtr,
  CepsVector<CepsChar>&            isBoundaryCell,
  CepsBool                         boundaryFlag
)
{
  vtkIdList *ptIds = vtkIdList::New();
  m_mesh->GetCellPoints(fileCellId,ptIds);
 
  CepsUInt nbPoints = boundaryFlag ? m_nbNodesPerBdryCell : m_nbNodesPerCell;
 
  // Add the cell only if there is a match between the read Ids and the expected number of points
  if (nbPoints == ptIds->GetNumberOfIds())
  {
    for (CepsUInt i = 0U; i < nbPoints; i++)
      nodeIndices[nodeIndicesPtr[cellId]+i] = ptIds->GetId(i)+nodeOffset;
  
    // Read cellAttributes
    for (CepsUInt i=0; i<m_cellArrayToRead.size(); ++i)
      cellAttr[cellAttrPtr[cellId]+i] = m_cellArrayToRead[i]->GetTuple1(fileCellId);
 
    cellAttrPtr   [cellId+1] = cellAttrPtr   [cellId] + m_cellArrayToRead.size();
    nodeIndicesPtr[cellId+1] = nodeIndicesPtr[cellId] + nbPoints;
    isBoundaryCell[cellId  ] = boundaryFlag;
 
    if (boundaryFlag)
      m_bdryCellsRead++;
    else
      m_cellsRead++;
  }
 
  ptIds->Delete();
}
 
void
VtkMeshReader::initializePvtu(const CepsString &fileName)
{
  vtkXMLPUnstructuredGridReader *reader = vtkXMLPUnstructuredGridReader::New();
  reader->SetFileName (fileName.c_str());
  reader->Update();
 
  reader->GetOutput()->Register(reader);
  m_mesh = vtkDataSet::SafeDownCast(reader->GetOutput());
  reader->Delete ();
}
 
void
VtkMeshReader::initializeLegacy(const CepsString &fileName)
{
  vtkNew<vtkDataSetReader> reader;
 
  reader->SetFileName (fileName.c_str());
  reader->SetReadAllScalars(1);
  reader->SetReadAllVectors(1);
  reader->SetReadAllNormals(1);
  reader->SetReadAllTensors(1);
  reader->SetReadAllFields (1);
 
  reader->Update();
 
  reader->GetOutput()->Register(reader);
  m_mesh = vtkDataSet::SafeDownCast(reader->GetOutput());
 
}
 
void
VtkMeshReader::determineCellDataToRead()
{
  vtkCellData  *cellData = m_mesh->GetCellData();
  vtkDataArray *array    = nullptr;
 
  CepsInt i;
  CepsInt nbCells = m_mesh->GetNumberOfCells();
  for (i = 0; i < cellData->GetNumberOfArrays(); ++i)
  {
    array = cellData->GetArray(i);
    // check only arrays of integer type
    switch (array->GetDataType())
    {
      case VTK_UNSIGNED_SHORT:
      case VTK_UNSIGNED_CHAR:
      case VTK_UNSIGNED_LONG:
      case VTK_UNSIGNED_INT:
      case VTK_ID_TYPE:
      case VTK_SHORT:
      case VTK_LONG:
      case VTK_INT:
      case VTK_LONG_LONG:
      {
        if ((array->GetNumberOfComponents()==1) and (array->GetNumberOfTuples() == nbCells))
        {
          // use this array as an attribute field
          m_nbAttrPerCell++;
          m_nbAttrPerBdryCell++;
          m_cellArrayToRead.push_back (array);
          CEPS_SAYS("        Detected attributes array \"" << array->GetName() << "\" on cells");
        }
        break;
      }
      default:
        // ignore this array
        break;
    }
  }
}
 
void
VtkMeshReader::determinePointDataToRead()
{
  vtkPointData *data  = m_mesh->GetPointData();
  vtkDataArray *array = nullptr;
 
  CepsInt i;
  for (i = 0; i < data->GetNumberOfArrays(); ++i)
  {
    array = data->GetArray(i);
    // check only arrays of integer type
    switch (array->GetDataType())
    {
      case VTK_UNSIGNED_SHORT:
      case VTK_UNSIGNED_CHAR:
      case VTK_UNSIGNED_LONG:
      case VTK_UNSIGNED_INT:
      case VTK_ID_TYPE:
      case VTK_SHORT:
      case VTK_INT:
      {
        if ((array->GetNumberOfComponents() == 1) and ((CepsUInt)(array->GetNumberOfTuples()) == m_nbNodes))
        {
          // use this array as an attribute field
          m_nbAttrPerNode++;
          m_pointArrayToRead.push_back(array);
          CEPS_SAYS("        Detected attributes array \"" << array->GetName() << "\" on points");
        }
        break;
      }
      default: break;
    }
  }
}