TetgenMeshReader.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 "TetgenMeshReader.hpp"
 
TetgenMeshReader::TetgenMeshReader(
  const CepsString &fileName,
  const CepsUInt   &dim,
  CepsBool          quads
) :
    MeshReader(fileName,dim,quads),
    m_indexOffset(0),
    m_bdryIndexOffset(0),
    m_bdryMarker(0)
{
  m_extension   = "ele";
  m_initialized = false;
  checkExtension();
  this->open ();
}
 
TetgenMeshReader::~TetgenMeshReader()
{
  if (m_file.is_open())
    m_file.close();
  if (m_bdryFile.is_open())
    m_bdryFile.close();
  if (m_nodesFile.is_open())
    m_nodesFile.close();
}
 
// redefinition of the FileReader open method
CepsBool
TetgenMeshReader::open()
{
  CEPS_SAYS("      Reading " << m_fileName << " and associated files");
  // The cell file
  m_file.open(m_fileName.c_str(), std::ifstream::in | std::ifstream::binary);
  // Check has already been done in FileReader::open()
 
  // The node file
  m_nodesFileName = ceps::changeExtension(m_fileName, "node");
  m_nodesFile.open (m_nodesFileName, std::fstream::in | std::ifstream::binary);
  CEPS_ABORT_IF(not m_nodesFile.good(),"Could not open mesh file " << m_nodesFileName);
 
  // Are there boundary cell files ?
  // check for faces first
  m_bdryFileName = ceps::changeExtension (m_fileName, "face");
  m_bdryFile.open (m_bdryFileName, std::fstream::in | std::ifstream::binary);
  if (!m_bdryFile.good())
  {
    // Did not find face file, maybe there is an edge file ?
    m_bdryFileName = ceps::changeExtension (m_fileName, "edge");
 
    m_bdryFile.open(m_bdryFileName, std::fstream::in | std::ifstream::binary);
    if (!m_bdryFile.good())
    {
      // No edge file, so no boundary cells
      CEPS_WARNS (
        "Tetgen mesh : found no boundary cell file " << std::quoted("*.face") << " or "
                                                     << std::quoted("*.edge") << " associated to "
                                                     << m_fileName
      );
    }
  }
 
  return true;
}
 
void
TetgenMeshReader::initialize()
{
  CepsString dum;
 
  // Nodes
 
  m_nodesFile.seekg(0, std::ios::beg);
  std::stringstream hn;
  hn << "Tetgen mesh reader: in file " << m_nodesFileName << "," << std::endl << "   ";
 
  m_nbNodes       = ceps::readInt(m_nodesFile,hn.str()+"unable to read number of nodes"          );
  m_nodesFile >> dum; // ignore dimensionality
  m_nbAttrPerNode = ceps::readInt(m_nodesFile,hn.str()+"unable to read nb of attributes per node");
  m_bdryMarker    = ceps::readInt(m_nodesFile,hn.str()+"unable to read boundary marker"          );
 
  m_nodesBegin = 1;
 
  // Cells
 
  m_file.seekg (0, std::ios::beg);
  std::stringstream hc;
  hc << "Tetgen mesh reader: in file " << m_fileName << "," << std::endl << "   ";
 
  m_nbCells        = ceps::readInt(m_file,hc.str()+"unable to read number of cells"              );
  m_nbNodesPerCell = ceps::readInt(m_file,hc.str()+"unable to read number of nodes per cell"     );
  m_nbAttrPerCell  = ceps::readInt(m_file,hc.str()+"unable to read number of attributes per cell");
  m_indexOffset    = ceps::readInt(m_file,hc.str()+"unable to read index offset"                 );
 
  // Check nbNodesPerCell with expected dimension
  CepsBool incorrectNNodes = (m_quads and m_nbNodesPerCell != pow (2, m_expectedDim)) 
                          or (not m_quads and m_nbNodesPerCell != m_expectedDim + 1);
  CEPS_ABORT_IF (incorrectNNodes,
    "Provided .ele file does not have elements with the expected number of nodes " << m_fileName
  );
 
  m_cellsBegin = 1;
 
  // Boundary cells
  if (m_bdryFile.is_open())
  {
    m_bdryFile.seekg (0, std::ios::beg);
    std::stringstream hb;
    hb << "Tetgen mesh reader: in file " << m_bdryFileName << "," << std::endl << "   ";
 
    m_nbBdryCells        = ceps::readInt(m_bdryFile,hb.str()+"unable to read number of boundary cells"              );
    m_nbNodesPerBdryCell = ceps::readInt(m_bdryFile,hb.str()+"unable to read number of nodes per boundary cell"     );
    m_nbAttrPerBdryCell  = ceps::readInt(m_bdryFile,hb.str()+"unable to read number of attributes per boundary cell");
    m_bdryIndexOffset    = ceps::readInt(m_bdryFile,hb.str()+"unable to read index offset of boundary cells"        );
    m_bdryCellsBegin = 1;
    // Check nbNodesPerCell with expected dimension
    incorrectNNodes = (m_quads and m_nbNodesPerBdryCell != pow (2, m_expectedDim - 1))
                   or (not m_quads and m_nbNodesPerBdryCell != m_expectedDim);
    CEPS_ABORT_IF (incorrectNNodes,
      "Provided .face or .edge file does not have elements with the expected number of nodes "
      << m_bdryFileName
    );
  }
  else
  {
    // no boundary
    m_nbBdryCells        = 0;
    m_nbNodesPerBdryCell = 0;
    m_nbAttrPerBdryCell  = 0;
    m_bdryIndexOffset    = 0;
    m_bdryCellsBegin     = 0;
  }
 
  m_nbCellsTotal = m_nbCells + m_nbBdryCells;
  m_initialized  = true;
  this->resetStreams();
}
 
CepsBool
TetgenMeshReader::hasBoundaryMarker() const
{
  return static_cast<CepsBool>(m_bdryMarker);
}
 
void
TetgenMeshReader::readNodes(
  CepsUInt                   indexStart,
  CepsUInt                   indexEnd,
  CepsUInt                   nodeOffset,
  CepsVector<CepsReal>&      coords,
  CepsVector<CepsAttribute>& attr,
  CepsVector<CepsUInt>&      attrPtr
)
{
  checkCanReadNodes(indexStart,indexEnd);
 
  CepsUInt nbToRead = indexEnd-indexStart;
 
  // place stream
  CepsString dummy;
  for (CepsUInt i=0U; i<indexStart; i++)
    getline(m_nodesFile,dummy);
 
  std::stringstream hh;
  hh << "Tetgen mesh reader: in file " << m_nodesFileName << "," 
     << std::endl << "   unable to read ";
  CepsString h = hh.str();
 
  for (CepsUInt i=0; i<nbToRead; i++)
  {
    CepsString n      = ceps::toString(i+1);
    CepsUInt   nodeId = nodeOffset+i;
    // Ignore index of node
    ceps::readInt(m_nodesFile,h+"expected node index "+n);
    for (CepsUInt j=0; j<3;j++)
      coords[nodeId*3+j] = ceps::readReal(m_nodesFile,h+"coordinates of point #"+n);
 
    for (CepsUInt j=0; j<(CepsUInt)m_nbAttrPerNode;j++)
      attr[attrPtr[nodeId]+j] = ceps::readInt(m_nodesFile,h+"attributes of point #"+n);
 
    attrPtr[nodeId+1] = attrPtr[nodeId] + m_nbAttrPerNode;
 
    // ignore boundary marker
    if (m_bdryMarker)
      ceps::readInt(m_nodesFile,h+"boundary marker");
  }
 
}
 
void
TetgenMeshReader::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
 
  CepsInt nbToRead = indexEnd-indexStart;
 
  // Place streams
  CepsInt nbCellsToRead    = 0;
  CepsInt nbBoundaryToRead = 0;
 
  // No cells to read, only boundary
  if (indexStart >= m_nbCells)
  {
    nbBoundaryToRead = nbToRead;
    CepsString dummy;
    for (CepsUInt i=0U; i < indexStart-m_nbCells; i++)
      getline(m_bdryFile, dummy);
  }
  else
  {
    nbCellsToRead    = m_nbCells-indexStart;
    nbBoundaryToRead = nbToRead-nbCellsToRead;
    CepsString dummy;
    for (CepsUInt i = 0U; i < indexStart; i++)
      getline (m_file, dummy);
  }
 
  // Actually read
  CepsInt i;
  for (i = 0; i < nbCellsToRead; i++)
    getNextCell(
      cellOffset+i,i,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,false
    );
 
  for (; i < nbBoundaryToRead+nbCellsToRead; i++)
    getNextCell(
      cellOffset+i,i-m_cellsRead,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,true
    );
 
  return;
}
 
void
TetgenMeshReader::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 << "Tetgen mesh reader: in file " << m_fileName << "," << std::endl << "   unable to read ";
  CepsString c = ceps::toString(fileCellId+1);
  CepsString h = hh.str();
 
  std::ifstream* file  = boundaryFlag ? &m_bdryFile : &m_file; 
  CepsUInt       nbPts = boundaryFlag ? m_nbNodesPerBdryCell : m_nbNodesPerCell;
  CepsUInt       nAttr = boundaryFlag ? m_nbAttrPerBdryCell  : m_nbAttrPerCell ;
 
  // ignore the index that comes first
  ceps::readInt(*file,h+"expected cell index "+c);
 
  for(CepsUInt i=0;i<nbPts;i++)
    nodeIndices[nodeIndicesPtr[cellId]+i] = ceps::readInt(*file,h+"node IDs of cell #"+c) 
                                              + nodeOffset - m_indexOffset;
 
  for (CepsUInt i=0;i<nAttr;i++) 
    cellAttr[cellAttrPtr[cellId]+i] = ceps::readInt(*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
TetgenMeshReader::resetStreams()
{
  CepsString temp;
  m_file.seekg(0,std::ios::beg);
  for (CepsInt i=0U; i<m_cellsBegin; i++)
    std::getline(m_file,temp);
 
  m_nodesFile.seekg(0,std::ios::beg);
  for (CepsInt i=0U; i<m_nodesBegin; i++)
    std::getline(m_nodesFile,temp);
 
  if (m_bdryFile.is_open())
  {
    m_bdryFile.seekg(0,std::ios::beg);
    for (CepsInt i=0U; i<m_bdryCellsBegin; i++)
      std::getline(m_bdryFile,temp);
  }
  return;
}