MeditGeometryWriter.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 "MeditGeometryWriter.hpp"
 
MeditGeometryWriter::
MeditGeometryWriter (const CepsString &fileName, Geometry *geom) :
  FileWriter (fileName),
  m_showGeometryPartitioning (false)
{
  m_geom = geom;
  CEPS_ABORT_IF(ceps::isNullPtr(m_geom),
    "Geometry passed to MeditGeometryWriter is nullptr"
  );
}
 
MeditGeometryWriter::
~MeditGeometryWriter()
{
}
 
void MeditGeometryWriter::
showGeometryPartitioning(CepsBool flag)
{
  m_showGeometryPartitioning = flag;
}
 
void MeditGeometryWriter::
write()
{
  if (ceps::isMaster ())
  {
    open();
    writeHeader();
    CEPS_SEPARATOR;
    CEPS_SAYS("Writing mesh in " << m_fileName);
  }
  writeNodes();
  writeCells();
  if (ceps::isMaster())
    close();
}
 
void MeditGeometryWriter::
writeHeader ()
{
  m_file << "MeshVersionFormatted 2" << std::endl;
  m_file << "Dimension" << std::endl << "3" << std::endl;
}
 
void MeditGeometryWriter::
writeNodes()
{
  CepsUInt nbNodes = m_geom->getNbNodes();
  MPI_Comm comm    = ceps::getCommunicator();
 
  CepsVector<CepsReal> coords(nbNodes*3,0.),myCoords(nbNodes*3,0.);
  CepsVector<CepsInt > attrs (nbNodes  ,0 ),myAttrs (nbNodes  ,0 );
 
  // Gather nodes data
  // FIXPERF
  // As the partitioning may not place the CPUs in order, it would be really cumbersome
  // to do an all gather. So for now, we can do ann AllReduce
  for (CepsUInt dim=1;dim<=3;dim++)
    if (m_geom->hasMeshOfDim(dim))
    {
      const CepsVector<GeomNode *> &nodes = m_geom->getMeshOfDim(dim)->getOwnedNodes();
      for (GeomNode* n: nodes)
      {
        CepsNodeGlobalIndex nID = n->getGlobalIndex();
        myCoords[3*nID  ] = n->x();
        myCoords[3*nID+1] = n->y();
        myCoords[3*nID+2] = n->z();
        if (m_showGeometryPartitioning)
          myAttrs[nID] = ceps::getRank();
        else
          // Medit format has only one attribute, so we use the first one only
          myAttrs[nID] = *(n->getAttributes().begin());
      }
    }
 
 
  MPI_Allreduce(myCoords.data(),coords.data(),nbNodes*3,CEPS_MPI_REAL,MPI_SUM,comm);
  MPI_Allreduce(myAttrs .data(),attrs .data(),nbNodes  ,CEPS_MPI_INT ,MPI_SUM,comm);
 
  // Now master writes the nodes
  // We do it with a buffer
  if (ceps::isMaster ())
  {
    m_file << std::endl << "Vertices" << std::endl << nbNodes << std::endl;
    std::stringstream toWrite (std::stringstream::out);
    toWrite.unsetf (std::ios::floatfield);
    toWrite.precision (15);
    CepsInt  counter = 0;
    for (CepsUInt i = 0; i < nbNodes; i++, counter++)
    {
      toWrite << coords[i*3  ] << " " 
              << coords[i*3+1] << " " 
              << coords[i*3+2] << " " 
              << attrs [i] << std::endl;
 
      if (counter > 4096)
      {
        m_file.write(toWrite.str().c_str(),toWrite.str().length());
        toWrite.str("");
        counter = 0;
      }
    }
    m_file.write(toWrite.str().c_str(),toWrite.str().length());
  }
}
 
void MeditGeometryWriter::
writeCells ()
{
  CepsBool surfacicCellsWritten = false;
  CepsBool cableCellsWritten    = false;
  CepsUInt nCells;
  Mesh*    mesh;
 
  // 3D cells
  if (m_geom->has3dMesh())
  {
    mesh   = m_geom->get3dMesh();
    nCells = mesh->getNbCells();
    if (ceps::isMaster())
      m_file << std::endl << "Tetrahedra" << std::endl << nCells << std::endl;
 
    this->writeCellsOf(mesh);
 
    // 2D boundary cells of 3d mesh
    nCells = mesh->getNbBoundaryCells();
    if (nCells)
    {
      // Maybe there is a 2d mesh, so get it's number of cells
      if (m_geom->has2dMesh())
        nCells += (m_geom->get2dMesh())->getNbCells();
 
      m_file << std::endl;
      m_file << "Triangles" << std::endl << nCells << std::endl;
      surfacicCellsWritten = true;
      this->writeCellsOf (mesh,true);
    }
  }
 
  // 2D cells of 2d mesh
  if (m_geom->has2dMesh ())
  {
    mesh   = m_geom->get2dMesh();
    nCells = mesh->getNbCells();
    if (!surfacicCellsWritten && ceps::isMaster())
      m_file << std::endl << "Triangles" << std::endl << nCells << std::endl;
 
    this->writeCellsOf(mesh);
 
    // 1D boundary cells of 2D mesh
    nCells = mesh->getNbBoundaryCells();
    if (nCells)
    {
      if (m_geom->has1dMesh())
        nCells += (m_geom->get1dMesh())->getNbCells();
 
      m_file << std::endl;
      m_file << "Edges" << std::endl << nCells << std::endl;
      cableCellsWritten = true;
      this->writeCellsOf(mesh,true);
    }
  }
 
  // 1D cells of 1d mesh
  if (m_geom->has1dMesh ())
  {
    mesh   = m_geom->get1dMesh();
    nCells = mesh->getNbCells();
    if (! cableCellsWritten && ceps::isMaster())
      m_file << std::endl << "Edges" << std::endl << nCells << std::endl;
 
    this->writeCellsOf(mesh);
  }
}
 
void MeditGeometryWriter::
writeCellsOf(Mesh *mesh, CepsBool boundary)
{
  CepsUInt  nbLocalCells    = boundary ? mesh->getLocalNbBoundaryCells() 
                                       : mesh->getLocalNbCells();
  auto      cells           = boundary ? mesh->getBoundaryCells()
                                       : mesh->getCells();
  CepsUInt  gridSize        = ceps::getGridSize();
  CepsBool  isMaster        = ceps::isMaster();
  CepsUInt  rank            = ceps::getRank();
  MPI_Comm  comm            = ceps::getCommunicator();
  
  CepsVector<CepsNodeGlobalIndex> cellNodeIndicesToSend;
  CepsVector<CepsInt>             cellAttrToSend;
  CepsVector<CepsCellGlobalIndex> cellIndicesToSend;
 
  // the array that will hold the number of local cells per process
  CepsUInt *nbCellsArray = new CepsUInt[gridSize];
  nbCellsArray[rank]     = nbLocalCells;
  // each process sends 1 CepsUInt to the root
  MPI_Gather(&nbLocalCells,1,CEPS_MPI_UINT,nbCellsArray,1,CEPS_MPI_UINT,0,comm);
 
  // Try to reserve some memory for better performance.
  // sizeHint is the number of nodes of the first cell
  CepsUInt sizeHint = 0;
  if (nbLocalCells)
    sizeHint = cells[0]->getNumberOfNodes();
 
  // If the meshes are very small (nb cells < nb process) then it is
  // possible that the root process does not own any cells. In that
  // case, sizeHint_t is zero, and we will have a problem when trying
  // to write the cells. To make sure that root process is aware of
  // the number of nodes per cell, exchange this information across
  // the computing grid.
  CepsUInt *sizeHintArray = new CepsUInt[gridSize];
  sizeHintArray[rank]   = sizeHint;
  MPI_Gather(&sizeHint,1,CEPS_MPI_UINT,sizeHintArray,1,CEPS_MPI_UINT,0,comm);
  if (isMaster)
  {
    CepsUInt max = sizeHintArray[0];
    for (CepsUInt i = 1U; i < gridSize; i++)
      if (sizeHintArray[i] > max)
        max = sizeHintArray[i];
    sizeHint = max;
  }
 
  cellNodeIndicesToSend.reserve (nbLocalCells*sizeHint);
  cellAttrToSend       .reserve (nbLocalCells);
  cellIndicesToSend    .reserve (nbLocalCells);
 
  CepsUInt size = cells.size ();
  GeomCell *cell;
  CepsVector<CepsNodeGlobalIndex> nodeIndices;
  for (CepsUInt i = 0U; i < size; i++)
  {
    nodeIndices.clear ();
    cell        = cells[i];
    nodeIndices = cell->getNodeIndices();
    // insert the node indices
    cellNodeIndicesToSend.insert(cellNodeIndicesToSend.end(),
                                 nodeIndices.begin(),nodeIndices.end());
    // push_back this cell's attribute. We can only take the first.
    cellAttrToSend.push_back(*(cell->getAttributes().begin()));
    // push_back this cell's globalIndex
    cellIndicesToSend.push_back((CepsCellGlobalIndex)cell->getGlobalIndex());
  }
 
  // Communicate the sizes of the vectors
 
  // Cell indices first
  int *cellIndicesToSendSizes  = new int[gridSize];
  int  localIndicesSize        = (int) cellIndicesToSend.size ();
  cellIndicesToSendSizes[rank] = localIndicesSize;
  MPI_Gather(&localIndicesSize,1,MPI_INT,cellIndicesToSendSizes,1,MPI_INT,0,comm);
 
  // Cell node indices
  int *cellNodeIndicesToSendSizes  = new int[gridSize];
  int  localNodeIndicesSize        = (int) cellNodeIndicesToSend.size ();
  cellNodeIndicesToSendSizes[rank] = localNodeIndicesSize;
  MPI_Gather(&localNodeIndicesSize,1,MPI_INT,cellNodeIndicesToSendSizes,1,MPI_INT,0,comm);
 
  // Allocate the space needed on master process
  // These variable will be relevant on root process only
  CepsUInt *cellIndices          = nullptr;
  CepsUInt *cellNodeIndices      = nullptr;
  CepsInt  *cellAttr             = nullptr;
  CepsUInt  cellIndicesTotalSize = 0, cellNodeIndicesTotalSize = 0;
  if (isMaster)
  {
    for (CepsUInt i = 0; i < gridSize; i++)
    {
      cellIndicesTotalSize += cellIndicesToSendSizes[i];
      cellNodeIndicesTotalSize += cellNodeIndicesToSendSizes[i];
    }
    cellIndices     = new CepsUInt[cellIndicesTotalSize];
    cellNodeIndices = new CepsUInt[cellNodeIndicesTotalSize];
    cellAttr        = new CepsInt[cellIndicesTotalSize];
  }
 
  // Now send the 'real' data: cell indices and cell nodes indices
 
  // Start with cell indices.
  // Create the displacement array
  int* displ = new int[gridSize];
  displ[0] = 0;
  if (isMaster)
    for (CepsUInt i = 1; i < gridSize; i++)
      displ[i] = displ[i-1] + cellIndicesToSendSizes[i-1];
 
  MPI_Gatherv(cellIndicesToSend.data(),cellIndicesToSendSizes[rank],CEPS_MPI_UINT,
              cellIndices,cellIndicesToSendSizes,displ,CEPS_MPI_UINT,0,comm);
 
  // similar for attributes
  MPI_Gatherv(cellAttrToSend.data(),cellIndicesToSendSizes[rank],CEPS_MPI_INT,
              cellAttr,cellIndicesToSendSizes,displ,CEPS_MPI_INT,0,comm);
 
  // Now, send the node indices.
  // Update displacement array
  if (isMaster)
    for (CepsUInt i = 1; i < gridSize; i++)
      displ[i] = displ[i-1] + cellNodeIndicesToSendSizes[i-1];
 
  MPI_Gatherv(cellNodeIndicesToSend.data(),cellNodeIndicesToSendSizes[rank],CEPS_MPI_UINT,
              cellNodeIndices,cellNodeIndicesToSendSizes,displ,CEPS_MPI_UINT,0,comm);
 
  if (isMaster)
  {
    // build the array of cell offsets
    CepsUInt *cellOffset = new CepsUInt[gridSize+1];
    cellOffset[0]      = 0;
    for (CepsUInt i = 1; i <= gridSize; i++) 
    {
      cellOffset[i] = cellOffset[i-1] + cellIndicesToSendSizes[i-1];
    }
 
    std::stringstream  toWrite(std::stringstream::out);
    std::set<CepsUInt> duplicates;
    CepsUInt           index, counter=0u;
    // iterate through our array, by parts
    for (CepsUInt j = 0; j < cellIndicesTotalSize; j++)
    {
      index = cellIndices[j];
      if (duplicates.find (index) == duplicates.end ())
      {
        // write it in new order
        for (CepsUInt k=0;k<sizeHint;k++)
        {
          CepsNodeGlobalIndex nID = cellNodeIndices[(j*sizeHint)+k];
          toWrite << nID+1 << " ";
        }
        // write a node attribute, either the owning proc index
        // or the true attribute
        if (m_showGeometryPartitioning)
          toWrite << this->determineCellOwner(cellOffset,j) << std::endl;
        else
          toWrite << cellAttr[j] << std::endl;
 
        duplicates.insert(index);
      }
      counter++;
 
      // Buffered writing: flush
      if (counter == 4096)
      {
        m_file.write(toWrite.str().c_str(),toWrite.str().length());
        toWrite.str ("");
        counter = 0;
      }
    }
    m_file.write(toWrite.str().c_str(),toWrite.str().length());
    ceps::destroyTabular(cellOffset);
    ceps::destroyTabular(displ);
  }
  // free memory
  ceps::destroyTabular(cellNodeIndicesToSendSizes);
  ceps::destroyTabular(cellIndicesToSendSizes);
 
  if (ceps::isMaster())
  {
    ceps::destroyTabular(cellIndices);
    ceps::destroyTabular(cellNodeIndices);
    ceps::destroyTabular(cellAttr);
  }
}
 
CepsUInt MeditGeometryWriter::
determineCellOwner(CepsUInt *cellOffset, CepsUInt index)
{
  CepsUInt gridSize = ceps::getGridSize ();
  CepsUInt proc = 0;
  while(proc<gridSize and index>= cellOffset[proc+1])
    proc++;
  return proc;
 
  //for (CepsUInt proc = 0; proc < gridSize; proc++)
  //  if (index < cellOffset[proc+1])
  //    // proc was the owner of the cell
  //    return proc;
 
  // NEVER REACHED
  //return 0;
}