MeditMeshReader.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 "MeditMeshReader.hpp"
 
MeditMeshReader::MeditMeshReader(const CepsString &fileName, const CepsUInt &dim, CepsBool quads) :
  MeshReader(fileName, dim, quads),
  m_indexOffset(0),
  m_boundaryIndexOffset(0)
{
  m_nbAttrPerCell     = 1;
  m_nbAttrPerBdryCell = 1;
  m_extension         = "mesh";
  checkExtension ();
}
 
MeditMeshReader::~MeditMeshReader()
{
  close();
}
 
void
MeditMeshReader::readNodes(
  CepsUInt                   indexStart,
  CepsUInt                   indexEnd,
  CepsUInt                   nodeOffset,
  CepsVector<CepsReal>&      coords,
  CepsVector<CepsAttribute>& attr,
  CepsVector<CepsUInt>&      attrPtr
)
{
  checkCanReadNodes(indexStart,indexEnd);
 
  // Skip to the correct line
  this->close();
  this->open ();
  CepsInt    nbToRead = indexEnd-indexStart;
  CepsString dummy;
  for (CepsUInt i=0U; i<indexStart+m_nodesBegin; i++)
    std::getline(m_file,dummy);
 
  std::stringstream h;
  h << "Medit mesh reader: in file " << m_fileName << "," << std::endl 
    << "   unable to read ";
 
  for (CepsInt i=0; i<nbToRead; i++)
  {
    CepsUInt   nodeId = nodeOffset+i;
    CepsString n      = ceps::toString(i+1);
    for (CepsUInt j=0; j<3;j++)
      coords[3*nodeId+j] = ceps::readReal(m_file,h.str()+"coordinates of point #"+n);
  
    for (CepsUInt j=0; j<(CepsUInt)m_nbAttrPerNode; j++)
      attr[attrPtr[nodeId]+j] = ceps::readInt(m_file,h.str()+"attributes of point #"+n);
  
    attrPtr[nodeId+1] = attrPtr[nodeId] + m_nbAttrPerNode;
  }
 
}
 
void
MeditMeshReader::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);
 
  // 0. Init
  this->reset();
  CepsInt nbToRead = indexEnd-indexStart;
 
  CepsInt nbCellsToRead    = 0;
  CepsInt nbBoundaryToRead = 0;
 
  // 1. Reach starting line:
 
  // If only boundary cells to read
  if (indexStart >= m_nbCells)
  {
    nbBoundaryToRead = nbToRead;
    CepsString dummy;
    for (CepsInt i=0; i<(CepsInt)(m_bdryCellsBegin+(indexStart-m_nbCells)); i++)
      std::getline(m_file,dummy);
  }
  // Both kinds of cells to read
  else
  {
    nbCellsToRead    = (CepsInt)(m_nbCells-indexStart);
    nbBoundaryToRead = nbToRead-nbCellsToRead;
    CepsString dummy;
    for (CepsInt i = 0; i < (CepsInt) (m_cellsBegin+indexStart); i++)
      std::getline(m_file,dummy);
  }
 
  // 2. Read cells (if any to read)
  CepsInt i;
  for (i=0; i<nbCellsToRead; i++)
    getNextCell(
      cellOffset+i,i,nodeOffset,nodeIndices,nodeIndicesPtr,cellAttr,cellAttrPtr,isBoundaryCell,false
    );
 
  // 2.5 Finished reading cells. Now place stream on the boundary cells
  if (indexStart <= m_nbCells)
  {
    CepsString dummy;
    this->reset();
    for (CepsInt j=0; j<(CepsInt)m_bdryCellsBegin; j++)
      std::getline(m_file,dummy);
  }
 
  // 3. Read boundary cells
  for (; i < nbBoundaryToRead+nbCellsToRead; i++)
    getNextCell(
      cellOffset+i,i,nodeOffset,nodeIndices,nodeIndicesPtr,cellAttr,cellAttrPtr,isBoundaryCell,true
    );
 
  return;
}
 
void
MeditMeshReader::getNextCell(
  CepsUInt                         cellId,
  CepsUInt                         fileCellId,
  CepsUInt                         nodeOffset,
  CepsVector<CepsNodeGlobalIndex>& nodeIndices,
  CepsVector<CepsUInt>&            nodeIndicesPtr,
  CepsVector<CepsAttribute>&       cellAttr,
  CepsVector<CepsUInt>&            cellAttrPtr,
  CepsVector<CepsChar>&            isBoundaryCell,
  CepsBool                         boundaryFlag
)
{
  std::stringstream hh;
  hh << "Medit mesh reader: in file " << m_fileName << "," << std::endl << "   unable to read ";
  CepsString c = ceps::toString(fileCellId+1);
  CepsString h = hh.str();
 
  CepsUInt nbPts = boundaryFlag ? m_nbNodesPerBdryCell : m_nbNodesPerCell;
  CepsUInt nAttr = boundaryFlag ? m_nbAttrPerBdryCell  : m_nbAttrPerCell ;
 
  for (CepsUInt i=0; i<nbPts; i++)
    nodeIndices[nodeIndicesPtr[cellId]+i] = ceps::readInt(m_file,h+"node IDs of cell #"+c) + nodeOffset - m_indexOffset;
 
  for (CepsUInt i=0;i<nAttr;i++) 
    cellAttr[cellAttrPtr[cellId]] = ceps::readInt(m_file,h+"attributes of cell #"+c);
 
  cellAttrPtr   [cellId+1] = cellAttrPtr   [cellId] + nAttr;
  nodeIndicesPtr[cellId+1] = nodeIndicesPtr[cellId] + nbPts;
 
  isBoundaryCell[cellId] = boundaryFlag;
 
  if (boundaryFlag)
    m_bdryCellsRead++;
  else
    m_cellsRead++;
}
 
 
void
MeditMeshReader::initialize()
{
  CEPS_SAYS ("      Reading " << m_fileName << " for dimension " << m_expectedDim);
 
  this->reset ();
 
  // find the vertices
  CepsString keyword ("Vertices");
  CepsInt    line = this->lineIndex (keyword);
  CEPS_ABORT_IF(line==-1, "Failed to find any nodes in mesh file.");
 
  // next line is nb nodes
  m_nbNodes = ceps::readInt(m_file,"Medit mesh reader: could not read number of nodes");
  // nodes begin at next line
  m_nodesBegin    = line + 2;
  m_nbAttrPerNode = 1;
 
  // try to find cells first
  CepsString kwV = m_quads ? m_keywordsQ[m_expectedDim] : m_keywordsT[m_expectedDim];
  CepsString kwB = m_quads ? m_keywordsQ[m_expectedDim - 1] : m_keywordsT[m_expectedDim - 1];
 
  this->reset ();
  line = this->lineIndex (kwV);
  CEPS_ABORT_IF (line==-1,
    "Could not find any cells of type \"" << kwV << "\" in current Medit mesh file." << m_fileName
  );
 
  // next line is nb cells
  m_nbCells = ceps::readInt(m_file,"Medit mesh reader: could not read number of cells");
  m_cellsBegin = line + 2;
 
  // look for boundary cell now
  this->reset ();
  line = this->lineIndex (kwB);
  if (line == -1)
  {
    // no boundary cells
    CEPS_WARNS ("No boundary cells found for mesh " << m_fileName);
    m_bdryCellsBegin = 0;
    m_nbBdryCells    = 0;
  }
  else
  {
    m_nbBdryCells = ceps::readInt(m_file,"Medit mesh reader: could not read number of bdry cells");
    m_bdryCellsBegin = line + 2;
  }
 
  m_indexOffset         = 1;
  m_boundaryIndexOffset = 1;
 
  this->reset ();
  m_nbCellsTotal = m_nbCells + m_nbBdryCells;
  m_initialized  = true;
}