GeometryReader.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 "geometry/readers/GeometryReader.hpp"
 
GeometryReader::GeometryReader (
  const CepsVector<CepsString> &volumicMeshes,
  const CepsVector<CepsString> &surfacicMeshes,
  const CepsVector<CepsString> &cableMeshes,
  const CepsString             &coupledNodesFile
) :
    m_totalNbNodes (0u),
    m_totalNbCells (0u),
    m_totalNbBdryCells (0u)
{
  // How many meshes do we have ?
  CepsUInt nVolumic  = volumicMeshes.size ();
  CepsUInt nSurfacic = surfacicMeshes.size ();
  CepsUInt nCable    = cableMeshes.size ();
 
  m_nbMeshes = nVolumic + nSurfacic + nCable;
 
  // Allocate structures
  m_meshFileNames.resize (m_nbMeshes);
  m_readers         = ceps::newArray<MeshReader *> (m_nbMeshes, nullptr);
  m_readerType      = ceps::newArray<CepsUInt> (m_nbMeshes, 0u);
  m_nodeOffsets     = ceps::newArray<CepsUInt> (m_nbMeshes + 1, 0u);
  m_cellOffsets     = ceps::newArray<CepsUInt> (m_nbMeshes + 1, 0u);
  m_meshNbNodes     = ceps::newArray<CepsUInt> (m_nbMeshes, 0u);
  m_meshNbCells     = ceps::newArray<CepsUInt> (m_nbMeshes, 0u);
  m_meshNbBdryCells = ceps::newArray<CepsUInt> (m_nbMeshes, 0u);
 
  CepsUInt index = 0U;
  // save mesh reader type, and build the reader
  for (CepsUInt i = 0U; i < nVolumic; i++, index++)
  {
    m_readerType[index] = MeshReaderType::VOLUMIC;
    m_readers[index]    = this->getMeshReaderFor (volumicMeshes[i], index, 3);
  }
  for (CepsUInt i = 0U; i < nSurfacic; i++, index++)
  {
    m_readerType[index] = MeshReaderType::SURFACIC;
    m_readers[index]    = this->getMeshReaderFor (surfacicMeshes[i], index, 2);
  }
  for (CepsUInt i = 0U; i < nCable; i++, index++)
  {
    m_readerType[index] = MeshReaderType::CABLE;
    m_readers[index]    = this->getMeshReaderFor (cableMeshes[i], index, 1);
  }
 
  // Read meshes
  initializeDataStructures();
 
  // Deal with coupled nodes and/or junctions
  m_hasCoupledNodes = not coupledNodesFile.empty();
  if (m_hasCoupledNodes)
    readCoupledNodes(coupledNodesFile);
 
  // Update cell counts
  for (CepsUInt i = 0U; i < m_nbMeshes; i++)
  {
    m_totalNbNodes     += m_meshNbNodes    [i];
    m_totalNbCells     += m_meshNbCells    [i];
    m_totalNbBdryCells += m_meshNbBdryCells[i];
  }
}
 
GeometryReader::~GeometryReader ()
{
  if (m_nbMeshes)
  {
    for (CepsUInt i=0U; i<m_nbMeshes; i++)
      ceps::destroyObject(m_readers[i]);
    ceps::destroyTabular(m_readers);
    ceps::destroyTabular(m_readerType);
    ceps::destroyTabular(m_nodeOffsets);
    ceps::destroyTabular(m_cellOffsets);
    ceps::destroyTabular(m_meshNbNodes);
    ceps::destroyTabular(m_meshNbCells);
    ceps::destroyTabular(m_meshNbBdryCells);
  }
}
 
void
GeometryReader::getCellOffsetsPerMeshType(CepsUInt cellOffsets[4])
{
  memset (cellOffsets, 0, 4*sizeof(CepsUInt));
  for (CepsUInt i=0; i<m_nbMeshes; i++)
  {
    // Fill next types with mesh nb of cells
    CepsUInt nbCells = m_readers[i]->getNbCells() + m_readers[i]->getNbBdryCells();
    for (CepsUInt j = m_readerType[i]+1; j<4; j++)
      cellOffsets[j] += nbCells;
  }
  return;
}
 
void
GeometryReader::getCellOffsetPerMesh(CepsVector<CepsUInt> &offsets)
{
  offsets.clear();
  offsets.resize(m_nbMeshes + 1);
  offsets[0] = 0U;
  for (CepsUInt i=0U; i<m_nbMeshes; i++)
    offsets[i+1] = offsets[i] 
                 + m_readers[i]->getNbCells() 
                 + m_readers[i]->getNbBdryCells();
  return;
}
 
void
GeometryReader::getNodeOffsetsPerMeshType(CepsUInt offsets[4])
{
  memset(offsets, 0, 4*sizeof(CepsUInt));
  for (CepsUInt i=0U; i<m_nbMeshes; i++)
  {
    CepsUInt nbNodes = m_readers[i]->getNbNodes ();
    for (CepsUInt j = m_readerType[i]+1; j<4; j++)
      offsets[j] += nbNodes;
  }
  return;
}
 
void
GeometryReader::getNodeOffsetPerMesh(CepsVector<CepsUInt> &offsets)
{
  offsets.clear();
  offsets.resize(m_nbMeshes+1);
  offsets[0] = 0U;
  for (CepsUInt i=0U; i<m_nbMeshes; i++)
    offsets[i+1] = offsets[i] + m_readers[i]->getNbNodes();
}
 
CepsUInt
GeometryReader::getNbNodes()
{
  return m_totalNbNodes;
}
 
CepsUInt
GeometryReader::getNbCells()
{
  return m_totalNbCells;
}
 
CepsUInt
GeometryReader::getNbBdryCells()
{
  return m_totalNbBdryCells;
}
 
CepsUInt
GeometryReader::readNodes(
  CepsUInt                   indexStart,
  CepsUInt                   nbNodesToRead,
  CepsVector<CepsReal>&      coords,
  CepsVector<CepsAttribute>& attr,
  CepsVector<CepsUInt>&      attrPtr,
  CepsUInt&                  attrSize
)
{
  // 0. No nodes to read. Abort
  if (nbNodesToRead < 1)
  {
    attrSize = 0;
    return 1;
  }
 
  // 1. -----------------------------------------------------
  // We must determine in which meshes we have to read nodes,
  // and for each mesh, where to start, where to end, and number
  // of node attributes
 
  CepsVector<CepsUInt> meshesToRead;
  CepsVector<CepsInt>  meshesStart;
  CepsVector<CepsInt>  meshesEnd;
  CepsVector<CepsUInt> meshesNbAttr;
 
  // 1.1 Where do we start ?
  CepsUInt iMesh;
  for (iMesh = 0; ((CepsInt)indexStart-(CepsInt)m_nodeOffsets[iMesh+1]) >= 0; iMesh++)
  {}
 
  CEPS_ABORT_IF(iMesh >= m_nbMeshes,
    "GeometryReader: starting node index is larger than number of nodes in all meshes"
  );
  
  meshesToRead.push_back(iMesh);
  meshesStart .push_back((CepsInt)(indexStart-m_nodeOffsets[iMesh]));
  meshesNbAttr.push_back(m_readers[iMesh]->getNbAttributesPerNode());
 
  // 1.2 Register next meshes
  CepsBool finished = false;
  CepsUInt indexEnd = indexStart + nbNodesToRead;
 
  while (not finished)
  {
    if ((CepsInt)(indexEnd-m_nodeOffsets[iMesh+1]) <= 0)
    {
      // We end with this mesh, read remaining nodes
      meshesEnd.push_back((CepsInt) (indexEnd-m_nodeOffsets[iMesh]));
      finished = true;
    }
    else
    {
      if (iMesh+1 == m_nbMeshes)
      {
        CEPS_ABORT("Number of nodes to read is too large!");
        return 0;
      }
 
      // Must read all this mesh remaining cells and boundary cells
      meshesEnd.push_back((CepsInt)(m_nodeOffsets[iMesh + 1] - m_nodeOffsets[iMesh]));
 
      // Onto next mesh
      meshesToRead.push_back(iMesh+1);
      meshesStart .push_back(0);
      meshesNbAttr.push_back(m_readers[iMesh+1]->getNbAttributesPerNode());
 
      // since we read contiguously:
      iMesh++;
    }
  }  // We now know what is going to be read
 
  // 2. -----------------------------------------------------
  // Prepare data arrays
 
  // Memory size for coords
  coords.resize(3 * nbNodesToRead);
 
  // Memory size of attr
  CepsUInt memorySize = 0U;
  for (CepsUInt i=0U; i<meshesToRead.size(); i++)
    memorySize += meshesNbAttr[i]*(meshesEnd[i]-meshesStart[i]);
 
  if (memorySize>0)
    attr.resize(memorySize, 0);
 
  attrPtr.resize(nbNodesToRead+1, 0);
  // Make sure its zero, just in case.
  attrPtr[0] = 0;
 
  // 3. -----------------------------------------------------
  // Now fill the arrays
 
  CepsUInt nbNodesRead(0u), start(0u), end(0u), toRead(0u);
  // Loop on meshes
  for (CepsUInt i=0U; i<meshesToRead.size(); i++)
  {
    toRead = meshesToRead[i];
    start  = meshesStart [i];
    end    = meshesEnd   [i];
    m_readers[toRead]->readNodes(start,end,nbNodesRead,coords,attr,attrPtr);
    nbNodesRead += end-start;
  }
  attrSize = memorySize;
  return 0;
}
 
CepsUInt
GeometryReader::readCells(
  CepsUInt                         indexStart,
  CepsUInt                         nbCellsToRead,
  CepsVector<CepsNodeGlobalIndex> &nodeIndices,
  CepsVector<CepsUInt>            &nodeIndicesPtr,
  CepsVector<CepsAttribute>       &attr,
  CepsVector<CepsUInt>            &attrPtr,
  CepsVector<CepsChar>            &isBoundary,
  CepsUInt                        &nNodeIndices,
  CepsUInt                        &nAttr
)
{
  // 0. ------------------------------------------------------------
  // Nothing to read
  if (nbCellsToRead < 1)
  {
    nAttr = 0;
    return 1;
  }
 
  // 1. ------------------------------------------------------------
  // Determine which mesh will be read
  CepsUInt indexEnd = indexStart + nbCellsToRead;
  CepsInt  offset   = indexStart;
 
  // build 2 arrays that will respectively hold the number of cells
  // and the number of boundary cells to read for each type
  CepsUInt             nbCells[3]     = {0, 0, 0};
  CepsUInt             nbBdryCells[3] = {0, 0, 0};
  CepsVector<CepsUInt> meshesToRead;
  CepsVector<CepsInt>  meshesStart;
  CepsVector<CepsInt>  meshesEnd;
 
  // 1.1 what is the first mesh we must read ?
  CepsUInt iMesh;
  for (iMesh=0U; ((CepsInt)indexStart-(CepsInt)m_cellOffsets[iMesh+1]) >= 0; iMesh++)
  {}
  CEPS_ABORT_IF (iMesh >= m_nbMeshes,
    "GeometryReader: starting cell index is larger than # of cells in all meshes"
  );
  offset -= m_cellOffsets[iMesh];
 
  meshesToRead.push_back (iMesh);
  meshesStart.push_back (offset);
 
  // 1.2 Register next meshes
  CepsBool finished     = false;
  CepsUInt cellAttrSize = 0U;
  while (not finished)
  {
    // Reading will be finished with this mesh
    if ((CepsInt) (indexEnd - m_cellOffsets[iMesh + 1]) <= 0)
    {
      CepsInt totalRemainingCells = indexEnd - m_cellOffsets[iMesh] - offset;
      CepsInt remainingBdry       = 0;
      CepsInt remainingCells      = 0;
      // Only boundary cells to read
      if (offset >= (CepsInt) m_meshNbCells[iMesh])
      {
        if (m_meshNbBdryCells[iMesh] == 0)
        {
          CEPS_ABORT ("Cannot read last boundary cells in mesh " << m_meshFileNames[iMesh]);
          return 0;
        }
        remainingBdry = totalRemainingCells;
      }
      // Cells and maybe boundary cells to read
      else
      {
        CepsInt temp = m_meshNbCells[iMesh] - offset;
        if (temp < totalRemainingCells)
        {
          // will read boundary
          remainingCells = temp;
          remainingBdry  = totalRemainingCells - temp;
        }
        else
        {
          // only cells
          remainingCells = totalRemainingCells;
        }
      }
 
      nbCells    [m_readerType[iMesh]] += remainingCells;
      nbBdryCells[m_readerType[iMesh]] += remainingBdry;
 
      cellAttrSize += remainingCells * m_readers[iMesh]->getNbAttributesPerCell    () +
                      remainingBdry  * m_readers[iMesh]->getNbAttributesPerBdryCell();
 
      meshesEnd.push_back(offset + totalRemainingCells);
      // We're done with this mesh
      finished = true;
    }
    else
    {
      // We will have to read all cells and boundary cells in this mesh,
      // starting from offset.
 
      CEPS_ABORT_IF(iMesh+1 == m_nbMeshes,
        "Number of nodes to read is too large!"
      );
 
      CepsInt remainingCells = 0;
      if (offset <= (CepsInt)m_meshNbCells[iMesh])
        remainingCells = m_meshNbCells[iMesh]-offset;
      CepsInt remainingBdry = m_cellOffsets[iMesh+1]-m_cellOffsets[iMesh]-offset-remainingCells;
 
      // Last cell (or boundary) of current mesh
      meshesEnd.push_back((CepsInt)m_cellOffsets[iMesh+1]-m_cellOffsets[iMesh]);
      // Next mesh
      meshesToRead.push_back(iMesh+1);
      meshesStart .push_back(0);
 
      // update cells to read
      nbCells    [m_readerType[iMesh]] += remainingCells;
      nbBdryCells[m_readerType[iMesh]] += remainingBdry;
 
      cellAttrSize += remainingCells*m_readers[iMesh]->getNbAttributesPerCell    () +
                      remainingBdry *m_readers[iMesh]->getNbAttributesPerBdryCell();
      // since we read contiguously:
      iMesh++;
      offset = 0;
    }
  }
 
  // 2. ------------------------------------------------------------
  // We now know what is going to be read
 
  // Memory size for node indices
  CepsUInt niSize = computeNumberOfNodeIndices(nbCells,nbBdryCells);
  nodeIndices   .resize(niSize);
  nodeIndicesPtr.resize(nbCellsToRead+1);
  isBoundary    .resize(nbCellsToRead);
 
  // Memory size of attr
  if (cellAttrSize > 0)
    attr.resize(cellAttrSize,0);
  attrPtr.resize(nbCellsToRead+1,0);
 
  // 3. ------------------------------------------------------------
  // Fill these arrays !
  CepsInt  nbCellsRead = 0;
  CepsUInt toRead(0u),start(0u),end(0u);
  // Loop on meshes
  for (CepsUInt i=0U; i<meshesToRead.size(); i++)
  {
    toRead = meshesToRead[i];
    start  = meshesStart [i];
    end    = meshesEnd   [i];
    m_readers[toRead]->readCells(
      start,end,m_nodeOffsets[toRead],m_cellOffsets[toRead],
      nodeIndices,nodeIndicesPtr,attr,attrPtr,isBoundary
    );
    nbCellsRead += end-start;
  }
 
  nNodeIndices = nodeIndices.size();
  nAttr        = cellAttrSize;
  return 0;
}
 
const CepsMultiMap<CepsNodeGlobalIndex, CepsNodeGlobalIndex> &
GeometryReader::getCoupledNodes ()
{
  return m_coupledNodes;
}
 
MeshReader *
GeometryReader::getMeshReaderFor (const CepsString &meshName, CepsInt index, CepsUInt dim)
{
  // Store string given by user
  m_meshFileNames.at (index) = meshName;
 
  CepsString dir      = ceps::getDir (meshName);
  CepsString fileName = ceps::getFilename (meshName);
 
  // Test missing extensions
  if (ceps::fileExists (fileName + ".mesh", dir))
    fileName.append (".mesh");
  else if (ceps::fileExists (fileName + ".ele", dir))
    fileName.append (".ele");
  else if (ceps::fileExists (fileName + ".msh", dir))
    fileName.append (".msh");
  else if (ceps::fileExists (fileName + ".pvtu", dir))
    fileName.append (".pvtu");
  else if (ceps::fileExists (fileName + ".vtk", dir))
    fileName.append (".vtk");
 
  // Resolve relative paths
  fileName = dir + "/" + fileName;
 
  CepsString ext = ceps::getExtension (fileName);
  if (ext == "mesh")
    return ceps::getNew<MeditMeshReader>(fileName, dim);
  // Tetgen format file
  else if (ext == "ele")
    return ceps::getNew<TetgenMeshReader>(fileName, dim);
  else if (ext == "msh")
    return ceps::getNew<GmshMeshReader>(fileName, dim);
  else if (ext == "vtk" or ext == "pvtu")
    return ceps::getNew<VtkMeshReader>(fileName, dim);
  else
    CEPS_ABORT ("Unsupported mesh format for file: " << fileName);
 
  return nullptr;
}
 
void
GeometryReader::initializeDataStructures ()
{
  for (CepsUInt i = 0; i < m_nbMeshes; i++)
  {
    m_readers[i]->initialize ();
    m_meshNbNodes[i]     = m_readers[i]->getNbNodes ();
    m_meshNbCells[i]     = m_readers[i]->getNbCells ();
    m_meshNbBdryCells[i] = m_readers[i]->getNbBdryCells ();
 
    // Update node, cell and boundary offsets
    m_nodeOffsets[i + 1] = m_nodeOffsets[i] + m_meshNbNodes[i];
    m_cellOffsets[i + 1] = m_cellOffsets[i] + m_meshNbCells[i] + m_meshNbBdryCells[i];
  }
}
 
void
GeometryReader::readCoupledNodes (const CepsString &fileName)
{
  CoupledNodesReader *r = ceps::getNew<CoupledNodesReader>(fileName);
  using couple_t        = std::pair<CepsNodeGlobalIndex, CepsNodeGlobalIndex>;
  for (NodeCoupling &nc : r->getCoupledNodesBetween (m_meshFileNames))
  {
    CepsInt  meshIndexA (ceps::argVector (m_meshFileNames, nc.getMesh (0)));
    CepsInt  meshIndexB (ceps::argVector (m_meshFileNames, nc.getMesh (1)));
    CepsUInt offsetA (m_nodeOffsets[meshIndexA]);
    CepsUInt offsetB (m_nodeOffsets[meshIndexB]);
    for (couple_t c : nc.getNodeCouples ())
    {
      // We store a bi-directional link
      m_coupledNodes.insert (couple_t (c.first + offsetA, c.second + offsetB));
      m_coupledNodes.insert (couple_t (c.second + offsetB, c.first + offsetA));
    }
  }
}
 
CepsUInt
GeometryReader::computeNumberOfNodeIndices (CepsUInt nbCells[3], CepsUInt nbBoundary[3])
{
  GeomCell *c3 = ceps::getNew<GeomSimplex<3u>>();
  GeomCell *c2 = ceps::getNew<GeomSimplex<2u>>();
  GeomCell *c1 = ceps::getNew<GeomSimplex<1u>>();
  GeomCell *c0 = ceps::getNew<GeomSimplex<0u>>();
 
  CepsUInt res = nbCells[0] * c3->getNumberOfNodes () +
         (nbCells[1] + nbBoundary[0]) * c2->getNumberOfNodes () +
         (nbCells[2] + nbBoundary[1]) * c1->getNumberOfNodes () +
         nbBoundary[2] * c0->getNumberOfNodes ();
 
  ceps::destroyObject(c3);
  ceps::destroyObject(c2);
  ceps::destroyObject(c1);
  ceps::destroyObject(c0);
 
  return res;
}