GmshMeshReader.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 "GmshMeshReader.hpp"
 
GmshMeshReader::GmshMeshReader(const CepsString &fileName, const CepsUInt &dim, CepsBool quads) :
    MeshReader (fileName,dim,quads)
{
 
  checkFormatVersion();
  m_extension = "msh";
  checkExtension();
 
  if (m_quads)
  {
    m_cellType     = m_keywordsQ[dim];
    m_bdryCellType = m_keywordsQ[dim-1];
  }
  else
  {
    m_cellType     = m_keywordsT[dim];
    m_bdryCellType = m_keywordsT[dim-1];
  }
}
 
GmshMeshReader::~GmshMeshReader()
{
  this->close ();
}
 
void
GmshMeshReader::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 << "Gmsh 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);
    // Ignore index of node
    ceps::readInt(m_file,h.str()+"expected node index "+n);
    for (CepsUInt j=0; j<3;j++)
      coords[3*nodeId+j] = ceps::readReal(m_file,h.str()+"coordinates of point #"+n);
  
    // No attribute to read in this file format, we put -1
    attr   [nodeId  ] = -1;
    attrPtr[nodeId+1] = attrPtr[nodeId] + 1;
  }
}
 
void
GmshMeshReader::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();
 
  // Reach starting line
  CepsString             line;
  CepsVector<CepsString> words;
  for (CepsInt i = 0; i < m_cellsBegin; i++)
    std::getline(m_file,line);
 
  // Cells and boundary cells can be mixed in the file, so we loop
  // twice in the list of cells
 
  // 1. Non-boundary cells
  for (CepsUInt i=0; i<m_nbCellsTotal; ++i)
    getNextCell(
      cellOffset+m_cellsRead,i,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,false
    );
 
  // 2. Boundary cells
  // Reach starting line (same as cells)
  this->reset();
  for (CepsInt i=0U; i<m_cellsBegin; i++)
    std::getline(m_file,line);
 
  for (CepsUInt i=0; i<m_nbCellsTotal; ++i)
    getNextCell(
      cellOffset+m_cellsRead+m_bdryCellsRead,i,nodeOffset,nodeIndices,nodeIndicesPtr,
      cellAttr,cellAttrPtr,isBoundaryCell,true
    );
 
  return;
}
 
void
GmshMeshReader::getNextCell(
  CepsUInt                         cellId,
  CepsUInt                         fileCellId,
  CepsUInt                         nodeOffset,
  CepsVector<CepsNodeGlobalIndex>& nodeIndices,
  CepsVector<CepsUInt>&            nodeIndicesPtr,
  CepsVector<CepsAttribute>&       cellAttr,
  CepsVector<CepsUInt>&            cellAttrPtr,
  CepsVector<CepsChar>&            isBoundaryCell,
  CepsBool                         boundaryFlag
)
{
  CepsString line;
  std::getline(m_file,line);
  auto words = ceps::split(line);
 
  CepsUInt   nbPoints = boundaryFlag ? m_nbNodesPerBdryCell : m_nbNodesPerCell;
  CepsString cellType = boundaryFlag ? m_bdryCellType       : m_cellType;
 
  // Cell line format:
  // cellID cellType Ntags <tags> <nodes list>
  if (words.size()>1 and words[1] == cellType)
  {
    std::stringstream hh;
    hh << "Gmsh mesh reader: in file " << m_fileName << "," << std::endl << "   unable to read ";
    CepsString h = hh.str();
    CepsString c = ceps::toString(fileCellId+1);
 
    std::istringstream iss(line);
    ceps::readInt(iss,h+"expected cell #"+c);
    ceps::readInt(iss,h+"type of cell for cell #"+c);
    CepsUInt nbAttrToRead = ceps::readInt(iss,h+"number of attributes of cell #"+c);
    CepsUInt nbAttr = boundaryFlag ? m_nbAttrPerBdryCell : m_nbAttrPerCell;
 
    for (CepsUInt i=0u; i<nbAttr; i++)
      cellAttr[cellAttrPtr[cellId]+i] = i<nbAttrToRead ? ceps::readInt(iss,h+"attribute of cell #"+c) :-1;
      
    for (CepsUInt i=0U; i<nbPoints; i++)
      nodeIndices[nodeIndicesPtr[cellId]+i] = ceps::readInt(iss,h+"node indices of cell #"+c)+nodeOffset-1;
    
    cellAttrPtr   [cellId+1] = cellAttrPtr   [cellId] + nbAttr;
    nodeIndicesPtr[cellId+1] = nodeIndicesPtr[cellId] + nbPoints;
    isBoundaryCell[cellId  ] = boundaryFlag;
 
    if (boundaryFlag)
      m_bdryCellsRead++;
    else
      m_cellsRead++;
  }
 
}
 
void
GmshMeshReader::initialize()
{
  CEPS_SAYS ("      Reading " << m_fileName << " for dimension " << m_expectedDim);
  std::stringstream hh;
  hh << "Gmsh mesh reader: in file " << m_fileName << "," << std::endl << "   ";
  CepsString h = hh.str();
 
  // Find the nodes
  this->reset();
  m_nodesBegin = this->lineIndex("$Nodes");
  CEPS_ABORT_IF(m_nodesBegin==-1, h+"\"$Nodes\" keyword not found");
 
  m_nbNodes = ceps::readInt(m_file,h+"unable to read the number of nodes");
  m_nodesBegin += 2;
  m_nbAttrPerNode = 1; // nodes don't have any attribute in msh format v2, we will put -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();
  m_cellsBegin = this->lineIndex ("$Elements");
  CEPS_ABORT_IF(m_cellsBegin==-1,h+"\"$Elements\" keyword not found");
 
  // Read cell types (second number on each cell line) and count for each type
  CepsString             line;
  CepsVector<CepsString> words, readElemTypes;
  m_nbCellsTotal      = 0u;
  m_nbCells           = 0u;
  m_nbBdryCells       = 0u;
  m_nbAttrPerCell     = 0u;
  m_nbAttrPerBdryCell = 0u;
  m_cellsBegin        = m_cellsBegin + 2;
  m_bdryCellsBegin    = m_cellsBegin;
 
  m_nbCellsTotal = ceps::readInt(m_file,h+"unable to read the number of cells");
  std::getline(m_file, line);  // End the line with number of cells
  for (CepsUInt i=0; i<m_nbCellsTotal; i++)
  {
    std::getline (m_file, line);
    words = ceps::split (line);
 
    try
    {
      CepsUInt nAttr(std::stoi(words[2]));
      if (words[1]==m_cellType)
      {
        m_nbCells++;
        m_nbAttrPerCell = std::max(m_nbAttrPerCell, nAttr);
      }
      else if (words[1]==m_bdryCellType)
      {
        m_nbBdryCells++;
        m_nbAttrPerBdryCell = std::max(m_nbAttrPerBdryCell, nAttr);
      }
    }
    catch (...)
    {
      CEPS_ABORT (h+"error while reading cell types and number of attributes");
    }
  }
 
  CEPS_ABORT_IF(not m_nbCells and not m_nbBdryCells,
    h+"could not find any cells for dimension " << m_expectedDim
  );
 
  this->reset();
  m_initialized = true;
}
 
void
GmshMeshReader::checkFormatVersion()
{
  this->reset();
  CepsInt index = this->lineIndex ("$MeshFormat");
  CEPS_ABORT_IF(index==-1,"Unable to determine mesh format in " << m_fileName);
 
  CepsString line;
  std::getline(m_file, line);
  CEPS_ABORT_IF(not line.starts_with ("2"),
    "Gmsh files format other than 2 are not supported by CEPS."
    << "Please select \"version 2\" when exporting your mesh in gmsh." << std::endl
    << "File : " << m_fileName
  );
}