CoeffReader.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 "pde/common/readers/CoeffReader.hpp"
 
CoeffReader::CoeffReader(const CepsString& fileName, const CepsLocationFlag& loc, Geometry* geom) :
    FileReader(fileName),
    m_geom(geom),
    m_loc(loc)
{
 
  CEPS_ABORT_IF(loc!=CepsLocationFlag::Cell and loc !=CepsLocationFlag::Point,
    "Given location of data must be either Cell or Node, not ZeroD"
  );
 
  CEPS_ABORT_IF(ceps::isNullPtr(geom), "Pointer on Geometry is null");
 
  // Coefficients over multiple meshes are not implemented
  // We would need to get the list of meshes from the geometry
  // and the offsets of nodes/cells to get the correct globalID
  CepsUInt count = 0;
  for (CepsUInt d = 1U; d <= 3U; d++)
    if (geom->hasMeshOfDim(d))
      count++;
 
  if (count != 1)
    CEPS_ABORT(
      "Importing coefficients defined on geometry is not implemented" << std::endl << 
      "   when geometry is made from multiple meshes"
    );
 
  m_fileName = fileName;
 
  CepsString ext = ceps::getExtension(fileName);
  CEPS_ABORT_IF(not isValidExtension(ext),
    "CoeffReader does not support files with extension " << ext << std::endl << 
    "   Supported extensions are .vtk and .msh"
  );
 
  buildArrayNamesList();
}
 
CoeffReader::~CoeffReader()
{
}
 
CepsBool
CoeffReader::isValidExtension(CepsString ext)
{
  return ext=="msh" or ext=="vtk";
}
 
void
CoeffReader::getData(const CepsString& label,DistributedInfos<CepsMathScalar>& dest)
{
  CEPS_ABORT_IF(not hasArray(label,0),
    "Reader has no scalar array with name " << label << " in file " << m_fileName
  );
  read(label,dest);
  return;
}
 
void
CoeffReader::getData(const CepsString& label,DistributedInfos<CepsMathVertex>& dest)
{
  CEPS_ABORT_IF(not hasArray(label,1),
    "Reader has no vector array with name " << label << " in file " << m_fileName
  );
  read(label,dest);
  return;
}
 
void
CoeffReader::getData(const CepsString& label,DistributedInfos<CepsMathTensor>& dest)
{
  CEPS_ABORT_IF(not hasArray(label,2),
    "Reader has no tensor array with name " << label << " in file " << m_fileName
  );
  read(label,dest);
  return;
}
 
CepsBool
CoeffReader::hasArray (const CepsString &label, CepsInt dim) const
{
  CepsBool found = false;
  CEPS_ABORT_IF(dim<-1 or dim>2,
    "Provided dimension must be either -1 or 0 1 2"
  );
  CepsUInt i   = dim==-1 ? 0 : dim;
  CepsUInt end = dim==-1 ? 3 : dim+1;
  while (i<end and not found)
  {
    if (m_loc==CepsLocationFlag::Cell)
      found = ceps::contains(m_cellArrayNames[i],label);
    if (m_loc==CepsLocationFlag::Point)
      found = ceps::contains(m_nodeArrayNames[i],label);
    i++;
  }
  return found;
}
 
void
CoeffReader::buildArrayNamesList ()
{
  CepsString ext = ceps::getExtension (m_fileName);
  if (ext == "vtk")
  {
    vtkDataSetReader *reader = vtkDataSetReader::New ();
    reader->SetFileName(m_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_vtkData = vtkDataSet::SafeDownCast(reader->GetOutput());
    reader->Delete();
 
    if (m_loc == CepsLocationFlag::Cell)
      for (CepsInt i = 0; i < m_vtkData->GetCellData()->GetNumberOfArrays(); i++)
        switch(m_vtkData->GetCellData()->GetArray(i)->GetNumberOfComponents())
        {
          case 1:
            m_cellArrayNames[0].push_back(m_vtkData->GetCellData()->GetArrayName(i));
            break;
          case 3:
            m_cellArrayNames[1].push_back(m_vtkData->GetCellData()->GetArrayName(i));
            break;
          case 9:
            m_cellArrayNames[2].push_back(m_vtkData->GetCellData()->GetArrayName(i));
            break;
        }
    else
      for (CepsInt i = 0; i < m_vtkData->GetPointData()->GetNumberOfArrays(); i++)
        switch(m_vtkData->GetPointData()->GetArray(i)->GetNumberOfComponents())
        {
          case 1:
            m_nodeArrayNames[0].push_back(m_vtkData->GetPointData()->GetArrayName(i));
            break;
          case 3:
            m_nodeArrayNames[1].push_back(m_vtkData->GetPointData()->GetArrayName(i));
            break;
          case 9:
            m_nodeArrayNames[2].push_back(m_vtkData->GetPointData()->GetArrayName(i));
            break;
        }
  }
  else  // msh
  {
    CepsArray<CepsString,2> keywords = {"$ElementData","$NodeData"};
    CepsArray<CepsArray<CepsVector<CepsString>,3>*,2> nameArrays = {&m_cellArrayNames,&m_nodeArrayNames};
 
    for (CepsUInt ii=0;ii<2;ii++)
    {
      CepsInt iLine = -1;
      FileReader::reset();
      do
      {
        iLine = lineIndex(keywords[ii]);
        if (iLine != -1)
        {
          CepsString line; 
          CepsString e = "in file " + m_fileName + ":\n   ";
          CepsInt    nStringTags,nRealTags,nComp;
 
          nStringTags = ceps::readInt(m_file,e+"cannot read number of string tags after " + keywords[ii]);
          for(CepsInt i=0;i<nStringTags+1;i++)
            std::getline(m_file,line);
          CepsString label = line.substr(line.find_first_of('\"')+1,line.find_last_of('\"')-1);
          // Get the nb of components
          skipLines(0);
          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");
          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");
          switch (nComp)
          {
            case 1:
              (*nameArrays[ii])[0].push_back(label);
              break;
            case 3:
              (*nameArrays[ii])[1].push_back(label);
              break;
            case 9:
              (*nameArrays[ii])[2].push_back(label);
              break;
          }
        }
      } while (iLine != -1);
    }
  }
}