Geometry.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/Geometry.hpp"
 
#include "geometry/partitioning/GeometryPartitioner.hpp"
#include "geometry/partitioning/PtscotchPartitioner.hpp"
 
// ---------------------------------------------------------------------------
// PUBLIC
// ---------------------------------------------------------------------------
 
// ---------------------------------------------------------------------------
// Create/destroy instance
 
Geometry::Geometry():
  CepsObject()
{
  initialize();
}
 
Geometry::Geometry(InputParameters *params) :
  CepsObject()
{
  initialize();
  setupWithParameters(params);
}
 
Geometry::~Geometry()
{
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim(dim))
      ceps::destroyObject(getMeshOfDim(dim));
}
 
void
Geometry::setupWithParameters(InputParameters *params)
{
  m_parameters = params;
 
  // Set meshes
  CepsVector<CepsString> meshes, keys = {"1dMesh", "2dMesh", "3dMesh"};
 
  for (CepsUInt dim=1; dim<=3; dim++)
    if (params->hasKey(keys[dim-1]))
    {
      meshes = ceps::split(params->getString(keys[dim-1]),",");
      for (CepsString &mesh : meshes)
        setMeshOfDim(dim,mesh);
    }
 
  CEPS_ABORT_IF(not has3dMesh() and not has2dMesh() and not has1dMesh(),
    "No meshes were provided in input file !"
  );
 
  // Connections between meshes
  if (params->hasKey("coupledNodesFile"))
    setCoupledNodesFile(params->getString("coupledNodesFile"));
 
  if (params->hasKey("unknowns per region"))
    setPartitionWeightsFromString(params->getString("unknowns per region"));
 
  // Compute domain partition
  setPartitioningMethodFromString(params->getString("partitioningMethod","PTScotchNode"));
  computePartition();
 
  CepsReal s = params->getReal("geometryScale",1.e0);
  if (not ceps::equals(s,1.e0))
    scale(s);
}
 
// ---------------------------------------------------------------------------
// Meshes management
 
void
Geometry::setMeshOfDim(CepsUInt dim, const CepsString &meshFileName)
{
  CEPS_ABORT_IF(dim==0 or dim>3,
    "Geometry::setMeshOfDim dim must be 1<=d<=3, got dim=" << dim
  );
  CepsString sd = ceps::toString(dim)+"D";
 
  if (not m_hasMeshOfDim[dim-1])
  {
    m_hasMeshOfDim[dim-1] = true;
    m_meshes      [dim-1] = ceps::getNew<Mesh>(dim);
    m_meshes      [dim-1]->setGeometry(this);
  }
 
  m_meshFilenames[dim-1].push_back(meshFileName);
}
 
void
Geometry::set3dMesh(const CepsString &meshFileName)
{
  setMeshOfDim(3,meshFileName);
}
 
void
Geometry::set2dMesh(const CepsString &meshFileName)
{
  setMeshOfDim(2,meshFileName);
}
 
void
Geometry::set1dMesh(const CepsString &meshFileName)
{
  setMeshOfDim(1,meshFileName);
}
 
Mesh*
Geometry::getMeshOfDim(CepsUInt dim) const
{
 
  CEPS_ABORT_IF(dim==0 or dim>3,
    "Geometry::getMeshOfDim dim must be 1<=d<=3, got dim=" << dim
  );
  return m_meshes.at(dim-1);
}
 
Mesh*
Geometry::get3dMesh() const
{
  return getMeshOfDim(3);
}
 
Mesh*
Geometry::get2dMesh() const
{
  return getMeshOfDim(2);
}
 
Mesh*
Geometry::get1dMesh() const
{
  return getMeshOfDim(1);
}
 
CepsBool
Geometry::hasMeshOfDim(CepsUInt dim) const
{
  CEPS_ABORT_IF(dim==0 or dim>3,
    "Geometry::hasMeshOfDim dim must be 1<=d<=3, got dim=" << dim
  );
  return m_hasMeshOfDim.at(dim-1);
}
 
CepsBool
Geometry::has3dMesh() const
{
  return hasMeshOfDim(3);
}
 
CepsBool
Geometry::has2dMesh() const
{
  return hasMeshOfDim(2);
}
 
CepsBool
Geometry::has1dMesh() const
{
  return hasMeshOfDim(1);
}
 
const CepsVector<CepsString>&
Geometry::getFileNamesForDim(CepsUInt dim) const
{
  return m_meshFilenames.at (dim-1);
}
 
const CepsVector<CepsString>&
Geometry::getVolumicFileNames() const
{
  return getFileNamesForDim(3);
}
 
const CepsVector<CepsString>&
Geometry::getSurfacicFileNames() const
{
  return getFileNamesForDim(2);
}
 
const CepsVector<CepsString>&
Geometry::getCableFileNames() const
{
  return getFileNamesForDim(1);
}
 
const CepsString&
Geometry::getCoupledNodesFileName() const
{
  return m_coupledNodesFileName;
}
 
void
Geometry::setCoupledNodesFile(const CepsString& coupledNodesFile)
{
  m_coupledNodesFileName = coupledNodesFile;
  m_hasCoupledNodes      = true;
}
 
// ---------------------------------------------------------------------------
// Partitioning
 
void
Geometry::setPartitioningMethodFromString(const CepsString& s)
{
  CepsString S = ceps::toUpper(s);
 
  if (S == "PTSCOTCHNODE")
    return setPartitioningMethod(PartitioningMethod::PTScotchNode);
 
  // Incorrect partitioning method was provided
  std::stringstream ss;
  ss << "Geometry Parameters: given partitioning method does not exist in CEPS." << std::endl
     << "   Valid keywords for partitioning methods are:" << std::endl;
  ss << "    - PTScotchNode" << std::endl;
  CEPS_ABORT(ss.str());
}
 
void
Geometry::setPartitioningMethod(PartitioningMethod method)
{
  m_partitioningMethod = method;
}
 
PartitioningMethod
Geometry::getPartitioningMethod() const
{
  return m_partitioningMethod;
}
 
void
Geometry::setPartitionWeights(const CepsMap<CepsAttribute,CepsUInt>& weights)
{
  m_partitionWeights = weights;
}
 
const CepsMap<CepsAttribute,CepsUInt>&
Geometry::getPartitionWeights() const 
{
  return m_partitionWeights;
}
 
CepsUInt
Geometry::computePartition()
{
  // Create a GeometryPartitioner that will do the job
  GeometryPartitioner* part = nullptr;
  switch(m_partitioningMethod)
  {
    case PartitioningMethod::PTScotchNode:
      CEPS_SAYS("Start reading and partitioning meshes with PtScotch");
      part = ceps::getNew<PtscotchPartitioner>(this);
      break;
    default:
      CEPS_ABORT("Unknown partitioning method");
  }
  CepsUInt err = part->computePartition();
 
  // Now get some data and share it if needed
 
  // Update max connectivity
  m_maxNodeConnectivity = this->computeMaxNodeConnectivity();
  m_maxCellDiameter     = this->computeMaxCellDiameter();
 
  // Update the number of halo nodes
  CepsUInt gridSize = ceps::getGridSize();
  if ((gridSize>1) && (m_nbHaloNodes==0))
  {
    // each process gets its local nb halo nodes
    CepsUInt localHaloNodes = this->getNbHaloNodesLocal();
 
    // exchange data with MPI
    CepsUInt globalHaloNodes = 0U;
 
    MPI_Allreduce(&localHaloNodes,&globalHaloNodes,1,CEPS_MPI_UINT,MPI_SUM, ceps::getCommunicator());
 
    m_nbHaloNodes = globalHaloNodes;
  }
 
  CepsUInt* cells     = ceps::newArray<CepsUInt>(gridSize);
  CepsUInt* boundary  = ceps::newArray<CepsUInt>(gridSize);
  CepsUInt* nodes     = ceps::newArray<CepsUInt>(gridSize);
  CepsUInt* haloNodes = ceps::newArray<CepsUInt>(gridSize);
  if (ceps::isProfiling())
  {
    // Print repartition of data
    this->getGlobalDistribution(cells,boundary,nodes,haloNodes);
 
    if (ceps::isDebug())
      this->printDataDistribution(cells,boundary,nodes,haloNodes,ceps::debugLog());
    if (ceps::isVerbose())
      this->printDataDistribution(cells,boundary,nodes,haloNodes,std::cout);
    this->printDataDistribution(cells,boundary,nodes,haloNodes,ceps::profilingLog());
  }
 
  #ifdef CEPS_DEBUG_ENABLED
    // Safety checks
    // this->getGlobalDistribution(cells,boundary,nodes,haloNodes);
    // CepsUInt nbCells = this->getNbCells();
    // CepsUInt nbNodes = this->getNbNodes();
    // CepsUInt i;
    // if (ceps::isMaster())
    // {
    //   CepsUInt sum = 0;
    //   // Sum of all nodes should be equal to total nb nodes
    //   for (i = 0; i<gridSize; i++)
    //     sum += nodes[i];
    //   CEPS_ABORT_IF (sum!=nbNodes,
    //     "Sum of local nodes on all process ("<< sum << ")" << 
    //     " differs from total nodes in geometry (" << nbNodes << ")."
    //   );
    //   // Check overlaping:
    //   sum = 0;
    //   for (i = 0; i < gridSize; i++)
    //     sum += cells[i];
    //   CEPS_ABORT_IF(sum<nbCells,
    //     "Invalid Geometry partitioning: incorrect nb cells per process."
    //   );
    // }
  #endif  // CEPS_DEBUG_ENABLED
  ceps::destroyTabular(cells    );
  ceps::destroyTabular(boundary );
  ceps::destroyTabular(nodes    );
  ceps::destroyTabular(haloNodes);
  ceps::destroyObject (part     );
 
  CEPS_SAYS("Information");
  CEPS_SAYS("   Number of nodes " << m_nbNodes );
  CEPS_SAYS("   Total number of cells " << m_nbCells + m_nbBoundaryCells);
  CEPS_SAYS("   Number of interior cells " << m_nbCells);
  CEPS_SAYS("   Number of boundary cells " << m_nbBoundaryCells);
  CEPS_SAYS("Geometry and partition set");
 
  return err;
}
 
// ---------------------------------------------------------------------------
// Number of nodes/cells
 
CepsUInt
Geometry::getNbCells() const
{
  return m_nbCells;
}
 
CepsUInt
Geometry::getNbBoundaryCells() const
{
  return m_nbBoundaryCells;
}
 
CepsUInt
Geometry::getNbNodes() const
{
  return m_nbNodes;
}
 
CepsUInt
Geometry::getNbHaloNodes() const
{
  return m_nbHaloNodes;
}
 
CepsUInt
Geometry::getNbOwnedNodesLocal() const
{
  CepsUInt nbNodes = 0;
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim (dim))
      nbNodes += m_meshes[dim-1]->getLocalNbNodes();
  return nbNodes;
}
 
CepsUInt 
Geometry::getNbHaloNodesLocal() const
{
  CepsUInt nbHaloNodes = 0;
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim(dim))
      nbHaloNodes += m_meshes[dim-1]->getLocalNbHaloNodes();
 
  return nbHaloNodes;
}
 
CepsUInt
Geometry::getNbNodesLocal() const
{
  return getNbOwnedNodesLocal() + getNbHaloNodesLocal();
}
 
CepsUInt
Geometry::getNbCellsLocal() const
{
  CepsUInt nbCells = 0;
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim (dim))
      nbCells += m_meshes[dim-1]->getLocalNbCells();
 
  return nbCells;
}
 
CepsUInt
Geometry::getNbBoundaryCellsLocal() const
{
  CepsUInt nbBoundary = 0;
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim(dim))
      nbBoundary += m_meshes[dim-1]->getLocalNbBoundaryCells();
 
  return nbBoundary;
}
 
CepsUInt
Geometry::getMaxNodeConnectivity() const
{
  return m_maxNodeConnectivity;
}
 
CepsReal
Geometry::getMaxCellDiameter() const
{
  return m_maxCellDiameter;
}
 
// ---------------------------------------------------------------------------
// Access to cellary entities
 
GeomNode*
Geometry::getNode(CepsNodeGlobalIndex globalID, CepsBool ownedOnly)
{
  GeomNode* node = nullptr;
  CepsUInt dim=1u;
  while (ceps::isNullPtr(node) and dim<=3)
  {
    if (hasMeshOfDim(dim))
    {
      if (ownedOnly)
        node = getMeshOfDim(dim)->getNode(globalID);
      else
        node = getMeshOfDim(dim)->getNodeOrHaloNode(globalID);
    }
    dim++;
  }
  return node;
}
 
GeomCell*
Geometry::getCell(CepsCellGlobalIndex globalID)
{
  GeomCell* elem = nullptr;
  for (CepsUInt dim=1; dim<=3; dim++)
    if (hasMeshOfDim(dim) and ceps::isNullPtr(elem))
    {
      elem = getMeshOfDim(dim)->getCell(globalID);
      if (ceps::isNullPtr(elem))
        elem = getMeshOfDim(dim)->getBoundaryCell(globalID);
 
      if (ceps::isValidPtr(elem))
        break;
    }
  return elem;
}
 
void
Geometry::print(std::ostream &os) const
{
  if (ceps::isMaster())
  {
    os << " Geometry" << std::endl;
    os << " ~~~~~~~~" << std::endl;
  }
 
  ceps::beginSequential();
  {
    os << "Process " << ceps::getRank () << " :" << std::endl;
 
    for (CepsUInt dim = 1; dim <= 3; dim++)
      if (hasMeshOfDim (dim))
      {
        os << "  " << dim << "D Mesh" << std::endl;
        os << "  -------" << std::endl;
        os << *getMeshOfDim (dim) << std::endl;
        os << std::endl;
      }
  }
  ceps::endSequential();
 
  if (ceps::isMaster())
    os << CEPS_STRING_SEPARATOR << std::endl;
  return;
}
 
std::ostream& operator<<
(std::ostream &os, const Geometry &g)
{
  g.print(os);
  return os;
}
 
// ---------------------------------------------------------------------------
// Data distribution
 
void 
Geometry::printDataDistribution(CepsUInt*     cells,
                                CepsUInt*     boundary,
                                CepsUInt*     nodes,
                                CepsUInt*     haloNodes,
                                std::ostream& os) const
{
  if (ceps::isMaster())
  {
    CepsUInt gridSize = ceps::getGridSize();
    CepsUInt i;
 
    CepsUInt totalCells     = this->getNbCells();
    CepsUInt totalBoundary  = this->getNbBoundaryCells();
    CepsUInt totalNodes     = this->getNbNodes();
    CepsUInt totalHaloNodes = this->getNbHaloNodes();
 
    // compute global overlapping
    CepsUInt sumCells = 0.0;
    for (i = 0; i < gridSize; i++)
      sumCells += cells[i];
 
    CepsReal ratHaloNodes = 100. * totalHaloNodes / (CepsReal) totalNodes;
    CepsReal overlapping  = 100 * (sumCells - totalCells) / (CepsReal) totalCells;
 
    os << "# " << CEPS_STRING_SEPARATOR << std::endl
       << "# Data Distribution" << std::endl
       << "# " << CEPS_STRING_SEPARATOR << std::endl
       << "# Global data:" << std::endl
       << "#\t- Number of nodes: " << totalNodes << std::endl
       << "#\t- Halo nodes: " << totalHaloNodes << " (" << ratHaloNodes << "% of nodes)" << std::endl
       << "#\t- Cells: " << totalCells << std::endl
       << "#\t- Boundary cells: " << totalBoundary << std::endl
       << "#\t- Theoretical optimum nodes per process: " << 100.0 / gridSize << "%" << std::endl
       << "#\t- Cell overlapping: " << overlapping << "%" << std::endl;
 
    // Avoid 0 division
    if (totalHaloNodes == 0)
      totalHaloNodes = 1;
 
    // CPU-wise info
    for (i = 0; i < gridSize; i++)
    {
      os << "# Process [" << i << "]" << std::endl;
      os << "#\t- Owned Nodes: " << nodes[i] << " (" << (nodes[i] / (CepsReal) totalNodes) * 100 << "%)" << std::endl;
      os << "#\t- Halo Nodes: " << haloNodes[i] << " (" << (haloNodes[i] / (CepsReal) totalHaloNodes) * 100 << "%)"
         << std::endl;
      os << "#\t- Cells: " << cells[i] << " (" << (cells[i] / (CepsReal) totalCells) * 100 << "%)"
         << std::endl;
      os << "#\t- Boundary Cells: " << boundary[i] << " (" << (boundary[i] / (CepsReal) totalBoundary) * 100 << "%)"
         << std::endl;
    }
    os << "# " << CEPS_STRING_SEPARATOR << std::endl;
  }
  return;
}
 
void
Geometry::getLocalDistribution(CepsUInt* nbCells,
                               CepsUInt* nbBoundaryCells,
                               CepsUInt* nbNodes,
                               CepsUInt* nbHaloNodes) const
{
  *nbCells         = this->getNbCellsLocal        ();
  *nbBoundaryCells = this->getNbBoundaryCellsLocal();
  *nbNodes         = this->getNbOwnedNodesLocal   ();
  *nbHaloNodes     = this->getNbHaloNodesLocal    ();
  return;
}
 
void 
Geometry::getGlobalDistribution(CepsUInt* cells,
                                CepsUInt* boundary,
                                CepsUInt* nodes,
                                CepsUInt* haloNodes) const
{
  // Fill array with local information
  CepsUInt nbCells, nbBdryCells, nbNodes, nbHaloNodes;
  getLocalDistribution(&nbCells, &nbBdryCells, &nbNodes, &nbHaloNodes);
 
  CepsUInt rank   = ceps::getRank();
  cells    [rank] = nbCells;
  boundary [rank] = nbBdryCells;
  nodes    [rank] = nbNodes;
  haloNodes[rank] = nbHaloNodes;
 
  // Share
  MPI_Gather(&nbCells    ,1,CEPS_MPI_UINT,cells    ,1,CEPS_MPI_UINT,0,ceps::getCommunicator());
  MPI_Gather(&nbBdryCells,1,CEPS_MPI_UINT,boundary ,1,CEPS_MPI_UINT,0,ceps::getCommunicator());
  MPI_Gather(&nbNodes    ,1,CEPS_MPI_UINT,nodes    ,1,CEPS_MPI_UINT,0,ceps::getCommunicator());
  MPI_Gather(&nbHaloNodes,1,CEPS_MPI_UINT,haloNodes,1,CEPS_MPI_UINT,0,ceps::getCommunicator());
  return;
}
 
// ---------------------------------------------------------------------------
// Indexing
 
const CepsVector<CepsSet<CepsNodeGlobalIndex>>&
Geometry::getNodesToSend()
{
  return m_nodesToSend;
}
 
const CepsVector<CepsSet<CepsNodeGlobalIndex>>&
Geometry::getNodesToReceive()
{
  return m_nodesToReceive;
}
 
const CepsVector<CepsNodeGlobalIndex>&
Geometry::getNodePartitionMap()
{
  return m_nodePartitionMap;
}
 
const CepsVector<CepsCellGlobalIndex>&
Geometry::getCellPartitionMap()
{
  return m_cellPartitionMap;
}
 
const CepsVector<CepsUInt>&
Geometry::getCellOffsetPerMesh()
{
  return m_cellOffsetPerMesh;
}
 
const CepsVector<CepsUInt>&
Geometry::getNodeOffsetPerMesh()
{
  return m_nodeOffsetPerMesh;
}
 
const CepsMultiMap<CepsNodeGlobalIndex,CepsNodeGlobalIndex>&
Geometry::getCoupledNodes()
{
  return m_coupledNodes;
}
 
// ---------------------------------------------------------------------------
// Coordinates
 
void
Geometry::scale(CepsReal scaleFactor)
{
  if (ceps::approxEquals (scaleFactor, 0.e0, 1.e-20))
    std::cout << "\033[1m\033[34mWARNING. \033[0mGeometry scale: scaling factor " << scaleFactor
              << " too small! Will not scale at all." << std::endl;
  else
  {
    // call scaling method on all available domains
    for (CepsUInt dim=1; dim<=3;dim++)
      if (hasMeshOfDim(dim))
        getMeshOfDim (dim)->scale(scaleFactor);
  }
}
 
// ---------------------------------------------------------------------------=
// PROTECTED
// ---------------------------------------------------------------------------=
 
CepsUInt
Geometry::computeMaxNodeConnectivity()
{
  if (not m_maxNodeConnectivity)
  {
    CepsUInt mymax = 0U;
    for (CepsUInt dim=1;dim<=3;dim++)
      if (hasMeshOfDim(dim))
        mymax = std::max(mymax,getMeshOfDim(dim)->getMaxNodeConnectivity());
    MPI_Allreduce(&mymax,&m_maxNodeConnectivity,1,CEPS_MPI_UINT,
                  MPI_MAX,ceps::getCommunicator());
  }
  return m_maxNodeConnectivity;
}
 
CepsReal
Geometry::computeMaxCellDiameter()
{
  if (m_maxCellDiameter<0.)
  {
    CepsReal mymax = -1.;
    for (CepsUInt dim=1;dim<=3;dim++)
      if (hasMeshOfDim(dim))
        mymax = std::max(mymax,getMeshOfDim(dim)->getMaxCellDiameter());
    MPI_Allreduce(&mymax,&m_maxCellDiameter,1,CEPS_MPI_REAL,
                  MPI_MAX,ceps::getCommunicator());
  }
  return m_maxCellDiameter;
}
 
// ---------------------------------------------------------------------------=
// PRIVATE
// ---------------------------------------------------------------------------=
 
void
Geometry::initialize()
{
  m_parameters          = nullptr;
  m_nbNodes             = 0;
  m_nbHaloNodes         = 0;
  m_nbCells             = 0;
  m_maxNodeConnectivity = 0;
  m_maxCellDiameter      = -1.;
  m_nbBoundaryCells     = 0;
  m_hasCoupledNodes     = false;
  setPartitioningMethod(PartitioningMethod::PTScotchNode);
 
  m_meshes       = {nullptr,nullptr,nullptr};
  m_hasMeshOfDim = {false  ,false  ,false  };
}
 
void
Geometry::setPartitionWeightsFromString(const CepsString& s)
{
  CepsVector<CepsString> spairs = ceps::split(s,",");
  for (CepsString spair: spairs)
  {
    std::istringstream iss(spair);
    CepsString em  = "Impossible to read number of unknowns in region with given attribute\n";
               em += "   The entry should be a comma separated list of pairs <attr,number>.\n";
               em += "   Got instead:\n   ";
               em += s; 
    CepsAttribute a = ceps::readInt(iss,em);
    CepsUInt      n = ceps::readInt(iss,em);
    m_partitionWeights.insert(std::make_pair(a,n));
  }
}
 
void
Geometry::updateMeshInfo(CepsUInt dim)
{
  // Get local info
  Mesh *mesh = getMeshOfDim(dim);
 
  if (mesh)
  {
    CepsUInt nNSCells    = mesh->getLocalNbCells        ();
    CepsUInt nNSBdCells  = mesh->getLocalNbBoundaryCells();
    CepsUInt nOwnedNodes = mesh->getLocalNbNodes        ();
 
    // Share
    CepsUInt gnNSCells(0U), gnNSBdCells(0U), gnOwnedNodes(0U);
    MPI_Allreduce(&nNSCells   ,&gnNSCells   ,1,CEPS_MPI_UINT,MPI_SUM,ceps::getCommunicator());
    MPI_Allreduce(&nNSBdCells ,&gnNSBdCells ,1,CEPS_MPI_UINT,MPI_SUM,ceps::getCommunicator());
    MPI_Allreduce(&nOwnedNodes,&gnOwnedNodes,1,CEPS_MPI_UINT,MPI_SUM,ceps::getCommunicator());
 
    // Udpate mesh info
    mesh->setNbCells        (gnNSCells   );
    mesh->setNbBoundaryCells(gnNSBdCells );
    mesh->setNbNodes        (gnOwnedNodes);
  }
  return;
}
 
void
Geometry::updateGeometryInfo()
{
  for (CepsUInt dim=1; dim<=3;dim++)
    if (hasMeshOfDim(dim))
      updateMeshInfo(dim);
  return;
}