CoeffReader.hpp

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
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#ifndef _COEFF_READER_HPP_
#define _COEFF_READER_HPP_
 
#include "common/CepsCommon.hpp"
#include "geometry/Geometry.hpp"
#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkDataSetReader.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
 
 
class CoeffReader : public FileReader
{
 
public:
 
  CoeffReader (const CepsString& fileName, const CepsLocationFlag& loc, Geometry* geom);
 
  ~CoeffReader () override;
 
  CepsBool
  isValidExtension(CepsString ext);
 
  void
  getData(const CepsString& label,DistributedInfos<CepsMathScalar>& dest);
 
  void
  getData(const CepsString& label,DistributedInfos<CepsMathVertex>& dest);
 
  void
  getData(const CepsString& label,DistributedInfos<CepsMathTensor>& dest);
 
  CepsBool
  hasArray(const CepsString &label, CepsInt dim=-1) const;
 
protected:
 
  void
  buildArrayNamesList();
 
  template<typename _Result>
  void
  read(const CepsString &label, DistributedInfos<_Result>& dest);
 
  template<typename _Result>
  void
  readVtk(const CepsString &label, DistributedInfos<_Result>& dest);
 
  template<typename _Result>
  void
  readMsh(const CepsString &label, DistributedInfos<_Result>& dest);
 
  // /// @brief Share data that has been read on master
  // void
  // broadCastData (const CepsString &label);
 
protected:
 
  Geometry*        m_geom; 
  CepsLocationFlag m_loc;  
 
  CepsArray<CepsVector<CepsString>,3> m_cellArrayNames;
  CepsArray<CepsVector<CepsString>,3> m_nodeArrayNames;
 
  vtkSmartPointer<vtkDataSet> m_vtkData;  
};
 
template<typename _Result>
void
CoeffReader::read(const CepsString &label, DistributedInfos<_Result>& dest)
{  
  CepsString ext = ceps::getExtension (m_fileName);
  CEPS_SAYS("  Reading coefficient \"" << label << "\" from file " << m_fileName);
  if (ext == "vtk")
    readVtk(label,dest);
  else 
    readMsh(label,dest);
 
  return;
}
 
template<typename _Result>
void
CoeffReader::readVtk (const CepsString &label, DistributedInfos<_Result>& dest)
{
 
  vtkIdType     nbData;
  vtkDataArray *array = nullptr;
 
  if (m_loc == CepsLocationFlag::Cell)
  {
    nbData = m_vtkData->GetNumberOfCells();
    array  = m_vtkData->GetCellData()->GetArray(label.c_str());
  }
  else
  {
    nbData = m_vtkData->GetNumberOfPoints();
    array  = m_vtkData->GetPointData()->GetArray(label.c_str());
  }
  // array should be found, as its presence was checked beforehand
 
  CepsReal* d = ceps::newArray<CepsReal>(9);
 
  for (vtkIdType i = 0; i < nbData; ++i)
  {
    array->GetTuple(i,d);
    // Node data, can be owned or halo
    // Even if the node is halo, we store the data in the "owned"
    // section of the distributed info. For a much faster access later on
    if (m_loc==CepsLocationFlag::Point)
    {
      CepsGlobalIndex gId  = m_geom->getNodePartitionMap()[i];
      GeomNode*       node = m_geom->getNode(gId);
      if (ceps::isValidPtr(node))
        dest.add(ceps::convertReal<_Result>(d),gId,gId,ceps::getRank());
    }
    // Cell data, owned only
    else 
    {
      CepsGlobalIndex gId  = m_geom->getCellPartitionMap()[i];
      GeomCell*       cell = m_geom->getCell(gId);
      if (ceps::isValidPtr(cell))
        dest.add(ceps::convertReal<_Result>(d),gId,gId,ceps::getRank());
    }
  }
  ceps::destroyTabular(d);
 
}
 
template<typename _Result>
void
CoeffReader::readMsh (const CepsString &label, DistributedInfos<_Result>& dest)
{
  // Check that the file contains an array defined where we want
  // with the correct name. We build a list of all arrays
  CepsInt    nRealTags,nIntTags,nComp,nToRead;
  CepsString e = "in file " + m_fileName + ":\n   ";
 
  // Reader "header" information of data field
  reset();
  std::stringstream keyword;
  keyword << "\"" << label << "\"";
  lineIndex(keyword.str());
 
  nRealTags   = ceps::readInt (m_file,e+"cannot read number of real tags");
  for (CepsInt i = 0U; i < nRealTags; i++)
    ceps::readReal(m_file,e+"cannot read real tag");
  nIntTags    = ceps::readInt (m_file,e+"cannot read number of int tags");
  ceps::readInt (m_file,e+"cannot read time step index");
  nComp       = ceps::readInt (m_file,"cannot read number of components of coefficients");
  nToRead     = ceps::readInt (m_file,"cannot read number of data to read");
  for (CepsInt i = 3U; i < nIntTags; i++)
    ceps::readInt (m_file,e+"cannot read int tag");
 
  // Read the data
  CepsReal* d = ceps::newArray<CepsReal>(9);
  CepsInt   eId;
  for (CepsInt n = 0U; n < nToRead; n++)
  {
    eId = ceps::readInt(m_file,e+"cannot read entity ID");
    for (CepsInt i = 0; i < nComp; i++)
    {
      d[i] = ceps::readReal(m_file,e+"could not read " + ceps::toString(i+1) + "-th value" 
                                    +" on node/cell with ID " + ceps::toString(eId));
 
    }
    // Node can be owned or halo
    if (m_loc==CepsLocationFlag::Point)
    {
      CepsGlobalIndex gId  = m_geom->getNodePartitionMap()[n];
      GeomNode*       node = m_geom->getNode(gId);
      if (ceps::isValidPtr(node))
      {
        dest.add(ceps::convertReal<_Result>(d),gId,gId,node->getOwner());
      }
    }
    // Cells are only owned
    else 
    {
      CepsGlobalIndex gId  = m_geom->getCellPartitionMap()[n];
      GeomCell*       cell = m_geom->getCell(gId);
      if (ceps::isValidPtr(cell))
      {
        dest.add(ceps::convertReal<_Result>(d),gId,gId,ceps::getRank());
      }
    }
  }
  ceps::destroyTabular(d);
 
}
 
#endif  //_COEFF_READER_HPP_