FiniteElements.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 "pde/finiteElements/FiniteElements.hpp"
#include "pde/problem/AbstractPdeProblem.hpp"
#include "pde/common/DistributedFactory.hpp"
#include "pde/assemblers/FE/FEMassAssembler.hpp"
#include "pde/assemblers/FE/FEIntegrator.hpp"
#include "pde/assemblers/FE/FEDivKGradAssembler.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
 
// ----------------------------------------------------------------------------------------
// Public methods
// ----------------------------------------------------------------------------------------
 
FiniteElements::FiniteElements(AbstractPdeProblem* problem, CepsUInt order) :
    AbstractDiscretization(problem),
    m_order               (order),
    m_refElems            ({}),
    m_nodes               (nullptr),
    m_volElems            (nullptr),
    m_bdrElems            (nullptr)
{
 
  m_spatialUnknownsLocation = CepsLocationFlag::Point;
 
  CEPS_SAYS("Creating Finite Elements...");
  if (ceps::isProfiling())
    getProfiler()->start("buildFE","building finite elements and degrees of freedom");
 
  // Initialize reference finite elements
  CEPS_SAYS_NOENDL("  1/4 Creating reference finite elements of order " << order << "...");
  initializeRefElements();
  CEPS_SAYS_INLINE(" Done.");
 
  ceps::barrier();
 
  // Create all DistributedInfos structures with nodes and elements
  CEPS_SAYS_NOENDL("  2/4 Building finite elements...");
  buildElements();
 
  ceps::barrier();
 
  computeNeighbors();
  // Synchronize all containers
  // The globalId of cells was already set from geometry global id (one geom cell = one finite elem)
  // For nodes, the loop through nodes is not the same as geom, and there may be additional nodes.
  // So we need to recompute the global indices of FE nodes.
  m_volElems->synchronize(false);
  m_bdrElems->synchronize(false);
  m_nodes   ->synchronize(true,&setGlobalIndexFunction<FENode>);
 
  CEPS_SAYS_INLINE(" Done.");
 
  ceps::barrier();
 
  // Build degrees of freedom tree from the problem data, then distributed infos and factory
  CEPS_SAYS_NOENDL("  3/4 Building degrees of freedom on FE nodes...");
  buildDofs();
  CEPS_SAYS_INLINE(" Done.");
 
  CEPS_SAYS_NOENDL("  4/4 Assembling finite elements mass and stiffness matrices...");
  m_mass  = newDofMatrix();
  m_stiff = newDofMatrix();
 
  FEMassAssembler     asbm(this);
  FEDivKGradAssembler asbs(this);
  for (Unknown* u: m_problem->getSpatialUnknowns())
    asbs.setKForUnknown(u,1.);
 
  asbm.setMatrix(m_mass .get());
  asbs.setMatrix(m_stiff.get());
  asbm.assemble();
  asbs.assemble();
  CEPS_SAYS_INLINE(" Done.");
 
  CEPS_SAYS("Finite elements ready");
  if (ceps::isProfiling())
    getProfiler()->stop("buildFE");
 
}
 
FiniteElements::~FiniteElements()
{
  for (auto& it : m_refElems)
    for (auto& ref: it.second)
      ceps::destroyObject(ref);
 
  m_volElems->destroyObjects();
  m_bdrElems->destroyObjects();
  m_nodes   ->destroyObjects();
  ceps::destroyObject(m_volElems);
  ceps::destroyObject(m_bdrElems);
  ceps::destroyObject(m_nodes);
}
 
CepsVector<FEBase*>&
FiniteElements::getFiniteElements()
{
  return m_volElems->getOwned();
}
 
CepsVector<FEBase*>&
FiniteElements::getBoundaryFiniteElements()
{
  return m_bdrElems->getOwned();
}
 
CepsUInt
FiniteElements::getNumberOfHaloFENodes() const
{
  return m_nodes->getNumberOfHalo();
}
 
CepsUInt
FiniteElements::getNumberOfOwnedFENodes() const
{
  return m_nodes->getNumberOfOwned();
}
 
const CepsVector<FENode*>&
FiniteElements::getOwnedFENodes()
{
  return m_nodes->getOwned();
}
 
const CepsVector<FENode*>&
FiniteElements::getHaloFENodes()
{
  return m_nodes->getHalo();
}
 
FENode*
FiniteElements::getFENode(CepsGlobalIndex nID)
{
  if (isOwnedFENode(nID))
    return m_nodes->getOwnedMapping()->atGlobal(nID);
 
  if (isHaloFENode(nID))
    return m_nodes->getHaloMapping()->atGlobal(nID);
 
  return nullptr;
}
 
CepsBool
FiniteElements::isOwnedFENode(CepsGlobalIndex nID) const
{
  return m_nodes->getOwnedMapping()->hasGlobal(nID);
}
 
CepsBool
FiniteElements::isHaloFENode(CepsGlobalIndex nID) const
{
  return m_nodes->getHaloMapping()->hasGlobal(nID);
}
 
CepsVector<DegreeOfFreedom*>
FiniteElements::getDofsOnElement(FEBase* elem)
{
  CepsVector<DegreeOfFreedom*> res;
  CepsUInt iloc = static_cast<CepsUInt> (CepsLocationFlag::Point);
  for (FENode* node: elem->getNodes())
    for (DegreeOfFreedom* dof: m_dofsTree[iloc][node->getGlobalIndex()])
      res.push_back(dof);
 
  return res;
}
 
CepsVector<DegreeOfFreedom*>
FiniteElements::getDofsOnElementForUnknowns(FEBase* elem, const CepsVector<Unknown*>& us)
{
  CepsVector<DegreeOfFreedom*> res = {};
 
  CepsUInt iloc = static_cast<CepsUInt> (CepsLocationFlag::Point);
  for (FENode* node: elem->getNodes())
    for (DegreeOfFreedom* dof: m_dofsTree[iloc][node->getGlobalIndex()])
      if (ceps::contains(us, dof->getUnknown()))
        res.push_back(dof);
 
  iloc = static_cast<CepsUInt> (CepsLocationFlag::Cell);
  for (DegreeOfFreedom* dof: m_dofsTree[iloc][elem->getProperties()->getGlobalIndex()])
      if (ceps::contains(us, dof->getUnknown()))
        res.push_back(dof);
  
  iloc = static_cast<CepsUInt> (CepsLocationFlag::ZeroD);
    for (auto l: m_dofsTree[iloc])
      for (DegreeOfFreedom* dof: l.second)
        if (ceps::contains(us, dof->getUnknown()))
          res.push_back(dof);
 
  return res;
}
 
CepsIndex
FiniteElements::getDofSpatialId(const DegreeOfFreedom* dof) const 
{
  CEPS_ABORT_IF(ceps::isNullPtr(dof),
    "nullptr dof"
  );
  CEPS_ABORT_IF(not dof->getUnknown()->isSpatial(),
    "you are trying to return a spatial index for a non-local "
    "unknown with name '"
    << dof->getUnknown()->getName() << "'"
  );
  return dof->getGeomId();
}
 
// --------------------------------------------------------------------------
 
DMatPtr
FiniteElements::getMassMatrix() const
{
  return m_mass;
}
 
DMatPtr
FiniteElements::getStiffnessMatrix() const
{
  return m_stiff;
}
 
 
// --------------------------------------------------------------------------
 
CepsReal
FiniteElements::dotProduct(
  DHVecPtr                      u,
  DHVecPtr                      v,
  CepsBool                      boundary,
  const CepsSet<CepsAttribute>& attrs,
  const CepsVector<Unknown*>&   unknowns
)
{
 
  CepsBool restrictions = boundary or not attrs.empty() or not unknowns.empty();
 
  if (not restrictions) { // Simply use the precomputed mass matrix
    DVecPtr tmp = newDofVector();
    tmp->mult(*m_mass,*u);
    return v->dot(*tmp);
  }
 
  // Create a special assembler for this purpose
  // FIXME this is higly memory and time inefficient.
  // We should implement integrators in parallel to assemblers, which, instead of
  // of storing the coefficients in the matrix, simply add the element contribution to a sum
  // {
    // FEMassAssembler asb(this,boundary,attrs,unknowns);
    // DMatPtr mat = newDofMatrix();
    // asb.setMatrix(mat.get());
    // asb.assemble();
    // tmp->mult(*mat,*u);
    FEMassIntegrator itg(this,boundary,attrs,unknowns);
    return itg.integrate(v,u);
  // }
 
 
}
 
CepsReal
FiniteElements::stiffProduct(
  DHVecPtr                      u,
  DHVecPtr                      v,
  CepsBool                      boundary,
  const CepsSet<CepsAttribute>& attrs,
  const CepsVector<Unknown*>&   unknowns
)
{
  CepsBool restrictions = boundary or not attrs.empty() or not unknowns.empty();
 
  DVecPtr tmp = newDofVector();
 
  if (not restrictions) { // Simply use the precomputed mass matrix
    tmp->mult(*m_stiff,*u);
  }
  else // Create a special assembler for this purpose
  // this is higly memory and time inefficient.
  {
    FEDivKGradAssembler asb(this,boundary,attrs,unknowns);
    for (Unknown* un: unknowns.empty() ? m_problem->getSpatialUnknowns() : unknowns)
      asb.setKForUnknown(un,1.);
    DMatPtr mat = newDofMatrix();
    asb.setMatrix(mat.get());
    asb.assemble();
    tmp->mult(*mat,*u);
  }
  return v->dot(*tmp);
}
 
CepsReal
FiniteElements::h1Norm(
  DHVecPtr                      v,
  CepsBool                      boundary,
  const CepsSet<CepsAttribute>& attrs,
  const CepsVector<Unknown*>&   unknowns
)
{
  return(std::sqrt(stiffProduct(v,v,boundary,attrs,unknowns)));
}
 
 
// --------------------------------------------------------------------------
// Protected methods
// --------------------------------------------------------------------------
 
void
FiniteElements::initializeRefElements()
{
  for (CepsCellType type : {CepsCellType::Simplex})
  {
    CepsFeLagrangeType fetype = {type,m_order};
    for (CepsUInt dim=0; dim<=3; dim++)
      m_refElems[type][dim] = ceps::getNew<ReferenceFE>(dim,fetype);
  }
  return;
}
 
ReferenceFE*
FiniteElements::getReferenceElementOfDim(CepsCellType type, CepsUInt dim) const
{
  return m_refElems.at(type).at(dim);
}
 
void
FiniteElements::buildElements()
{
  m_nodes    = new DistributedInfos<FENode*,CepsHash3>();
  m_volElems = new DistributedInfos<FEBase*>();
  m_bdrElems = new DistributedInfos<FEBase*>();
 
  CepsUInt nVolCells = 0u, nBdrCells = 0u;
  for (CepsUInt dim=1;dim<=3;dim++)
    if (m_geometry->hasMeshOfDim(dim))
    {
      nVolCells += m_geometry->getMeshOfDim(dim)->getLocalNbCells();
      nBdrCells += m_geometry->getMeshOfDim(dim)->getLocalNbBoundaryCells();
    }
  m_volElems->getOwned().reserve(nVolCells);
  m_bdrElems->getOwned().reserve(nBdrCells);
 
  // Builds volumic and boundary elements
  for (CepsUInt dim=1; dim<=3; dim++)
  {
    if (m_geometry->hasMeshOfDim(dim))
    {
      buildElements(m_geometry->getMeshOfDim(dim),false);
      buildElements(m_geometry->getMeshOfDim(dim),true );
    }
  }
 
}
 
void
FiniteElements::buildElements(Mesh* mesh, CepsBool onBoundary)
{
 
  // Some shortcuts
  ReferenceFE*                 refElem = nullptr;
  LocalGlobalMapping<FEBase*>* toBuild = nullptr;
 
  if (onBoundary)
    toBuild = m_bdrElems->getOwnedMapping();
  else
    toBuild = m_volElems->getOwnedMapping();
 
  CepsCellGlobalIndex cellId;
  CepsBool            foundSameNode;
  // Loop on geometrical cells
  for (GeomCell *cell : onBoundary ? mesh->getBoundaryCells() : mesh->getCells())
  {
    cellId = cell->getGlobalIndex();
    // Get the correct reference element for this geometrical element
    refElem = getReferenceElementOfDim(cell->getCellType(),cell->getDimension());
 
    // Create a FiniteElement with the same global index than geometrical cell.
    // Order is set from reference element.
    FEBase* fElem = ceps::getNew<FEBase>(cell,false,refElem);
 
    // Create the missing FENodes for that element
    const CepsUInt &n = refElem->getNumberOfBasisVertices();
 
    for (CepsUInt i=0U; i<n; ++i)
    {
      // Check if the node was not already added in a neighbouring element
      foundSameNode = false;
      for (CepsCellGlobalIndex neighId : cell->getNodeSharingCellsIndices())
      {
        FEBase* nElem = nullptr;
        if (m_volElems->hasGlobal(neighId,ceps::getRank()))
          nElem = m_volElems->atGlobal(neighId,ceps::getRank());
        else if (m_bdrElems->hasGlobal(neighId,ceps::getRank()))
          nElem = m_bdrElems->atGlobal(neighId,ceps::getRank());
 
        if (ceps::isValidPtr(nElem))
          checkNeighborsBeforeAssign(fElem,nElem,i,onBoundary,foundSameNode);
 
        if (foundSameNode)
          break;
      }
 
      if (not foundSameNode)
        createNewNode(fElem, i, onBoundary, mesh);
    }
 
    toBuild->append(fElem,cellId,cellId);
  }
  return;
}
 
void
FiniteElements::checkNeighborsBeforeAssign(
  FEBase*   fElem,
  FEBase*   nElem,
  CepsUInt  vertexId,
  CepsBool  boundary,
  CepsBool& foundSameNode
)
{
 
  ReferenceFE*        refFElem = fElem->getReferenceElement();
  CepsReal3D          p        = refFElem->getTransformedBasisVertex(vertexId,fElem->getGeomObject());
 
  // Loop through all FE nodes of neighbor element that were already created (!=nullptr)
  for (FENode* node : nElem->getValidNodes())
  {
 
    CepsReal3D q = node->getGeomObject()->getCoordinates();
 
    // p is obtained by multiplication by jacobian, so it is prone to floating point errors
    // We check equality to reach the same level of precision than CepsHash3
    foundSameNode = areSameArray(p,q,CEPS_HASH3_PRECISION/std::numeric_limits<CepsReal>::epsilon());
 
    if (foundSameNode)
    {
      node->getProperties()->setOnBoundary(boundary);
      node->addCell(fElem);
      fElem->setNodeAt(vertexId,node);
      return;
    }
  }
 
  return;
}
 
void
FiniteElements::createNewNode(FEBase* el, CepsUInt vertexId, CepsBool isBoundary, Mesh* mesh)
{
  const ReferenceFE*     refElement      = el->getReferenceElement();
  GeomCell*              gCell           = el->getGeomObject();
  CepsInt                feNodeCellIndex = refElement->getBasisVertexGeomCellIndex(vertexId);
  CepsSet<CepsAttribute> attributes      = {};
  CepsUInt               owner           = 0U;
  FENode*                fenode          = nullptr;
 
  // if the node is built on a geom node, we can directly use this geom node to get some information
  if (feNodeCellIndex!=-1)
  {
    const GeomNode* gNode = gCell->getNodeAt(feNodeCellIndex);
    attributes            = gNode->getAttributes();
    owner                 = gNode->getOwner();
    fenode                = ceps::getNew<FENode>(const_cast<GeomNode*>(gNode),false);
 
  }
  // else, we need to create these informations
  else
  {
 
    // Compute coordinates of new node
    CepsReal3D p = refElement->getTransformedBasisVertex(vertexId,gCell);
    // We need to determine the owner and attributes of this new node. 
    // The rule is given by the reference element.
    // If the node is -on a face : get data from node with min geom index on the face
    //                -inside the cell : use cell data (normally owner is of node with min geomid)
 
    GeomNode* ownerNode = refElement->getNodeOfOwnerForBasisVertex(vertexId,gCell);
 
    // We give the global index -1 to say that it is not part of the initial geometry
    GeomNode* gNode = ceps::getNew<GeomNode>(p[0],p[1],p[2],-1);
 
    // Set owner
    owner = ownerNode->getOwner();
    gNode->setOwner(owner);
 
    // Set attributes
    attributes = ownerNode->getAttributes();
    gNode->setAttributes(attributes.begin(),attributes.end());
 
    fenode = ceps::getNew<FENode>(gNode,true);
  }
 
  fenode->addCell(el);
  fenode->getProperties()->setOnBoundary(isBoundary);
 
  el->setNodeAt(vertexId,fenode);
 
 
  m_nodes->add(fenode,getNodeHash(fenode),owner);
 
  return;
}
 
void
FiniteElements::computeNeighbors()
{
  // FE Nodes
  for (auto &all : {m_volElems, m_bdrElems})
    for (FEBase* el : all->getOwned())
    {      
      CepsVector<FENode*> nodes = el->getNodes();
      CepsUInt n = nodes.size();
      for (CepsUInt i=0U; i<n; ++i)
        for (CepsUInt j=0U; j<n; ++j)
          if (i!=j)
          {
            if (ceps::argVector(nodes[i]->getNodes(),nodes[j])==-1)
            {
              nodes[i]->addNode(nodes[j]);
            }
          }
    }
 
  // Elements
  for (auto &list : {m_nodes->getOwned(), m_nodes->getHalo()})
    for (FENode* nd : list)
    {
      CepsVector<FEBase*> elems = nd->getCells();
      CepsUInt n = elems.size();
      for (CepsUInt i=0U; i<n; ++i)
        for (CepsUInt j=i+1; j<n; ++j)
          if (ceps::argVector(elems[i]->getCells(),elems[j])==-1)
            elems[i]->addCell(elems[j]);
    }
 
  return;
}
 
CepsHash3
FiniteElements::getNodeHash(const FENode* n) const
{
  return ceps::hash::getHash3(n->getGeomObject()->getCoordinates());
}
 
void 
FiniteElements::buildDofsTree()
{
 
  // Set unknowns to be computed on points (FENodes)
  for (Unknown* u: m_problem->getUnknowns())
    if (u->isSpatial())
      registerSpatialUnkown(u);
 
  // Loops on FE Nodes, owned then halo. 
  // A dof for each unknown defined in at least one of a neighboring cell will be created. 
  // Owned/halo distribution for dofs follows the distribution of FE Nodes.
  // Unknowns are defined on points, but the discrete variational formulation involves
  // integrals on cells
  CepsUInt iloc = static_cast<CepsUInt> (CepsLocationFlag::Point);
  for (const CepsVector<FENode*>& toParse : {m_nodes->getOwned(), m_nodes->getHalo()})
  {
    for (FENode* feNode : toParse)
    {
 
      // Build the set of attributes in cells around this node
      // We also add the cell attributes to the dofs afterwards
      // CepsSet<CepsAttribute> attrs; 
      // for (FEBase* nElem : feNode->getCells())
      //   for(CepsAttribute a: nElem->getGeomObject()->getAttributes())
      //     attrs.insert(a);
      auto attrs = feNode->getGeomObject()->getAttributes();
 
      CepsNodeGlobalIndex feId = feNode->getGlobalIndex();
      CepsNodeGlobalIndex gId  = feNode->getGeomObject()->getGlobalIndex();
      // Now check matching unknowns and create dof when needed
 
      for (Unknown* u: m_problem->getUnknowns())
      {
        if  (u->isSpatial() and (u->hasOneOfAttributes(attrs) or u->hasAttribute(CepsUniversal)))
        {
          auto dof = ceps::getNew<DegreeOfFreedom>(u,attrs,gId,0x0,feNode->getGeomObject()->getOwner());
          dof->setDiscrId(feId);
          dof->setVertex(feNode->getGeomObject());
          dof->setOnBoundary(feNode->getGeomObject()->isBoundary()); 
          m_dofsTree[iloc][feId].push_back(dof);
          // The global id of the dof is now 0, 
          // but it will be synchronized when calling buildDofsDistributedInfo
        }
      }
    }
  }
}
 
void 
FiniteElements::setSpatialUnknownsDofsInteractions()
{
 
  for (DistributedInfos<FEBase*>* toParse : {m_volElems, m_bdrElems})
  for (FEBase* elem: toParse->getOwned())
  {
    auto dofs = getDofsOnElement(elem);
    for (DegreeOfFreedom* dof1 : dofs)
    for (DegreeOfFreedom* dof2 : dofs)
    if (dof1 != dof2)
    {
 
      Unknown* u1 = dof1->getUnknown();
      Unknown* u2 = dof2->getUnknown();
 
      // Couple the dofs if they are for the same unknown
      CepsBool coupled = u1==u2;
      // Otherwise check in interactions defined in the problem
      CepsUInt i            = 0;
      auto     interactions = m_problem->getUnknownsInteractions();
      while (not coupled and i<interactions.size())
      {
        CepsBool okAttribute = interactions[i]->hasAttribute(CepsUniversal) 
                            or interactions[i]->getNumberOfAttributes()==0
                            or elem->getGeomObject()->hasOneOfAttributes(interactions[i]->getAttributes());
 
        coupled = okAttribute and 
                  ((u1==interactions[i]->u1 and u2==interactions[i]->u2) or 
                   (u1==interactions[i]->u2 and u2==interactions[i]->u1));
        i++;
      }
 
      if (coupled)
      {
        dof1->addNeighbor(dof2);
        dof2->addNeighbor(dof1);
 
        // If both dofs belong to another processor, there is a case in 3D where 
        // this other processor does not know the existence of the "edge" between the two dofs
        // (all cells around this edge can belong to another processor than the owner of the dofs)
        // We store this relationship, so that can give a safe non-zero stucture for allocation to petsc
        CepsProcId owner1 = dof1->getOwner();
        CepsProcId owner2 = dof2->getOwner();
        if (owner1==owner2 and owner1!=ceps::getRank())
        {
          m_extraAdjacencyToSend[owner1][dof1->getGlobalIndex()].insert(dof2->getGlobalIndex());
          m_extraAdjacencyToSend[owner1][dof2->getGlobalIndex()].insert(dof1->getGlobalIndex());
        }
      }
    }
  }
 
  // Commmunicate the extra adjacency 
  MPI_Status  status;
  MPI_Request request;
  MPI_Comm    comm     = ceps::getCommunicator();
  CepsUInt    gridSize = ceps::getGridSize();
  CepsProcId  rank     = ceps::getRank();
 
  // For now, slow communication matrix, needed barrier at each comm (???)
  for (CepsProcId sndr=0; sndr<gridSize; sndr++)
  for (CepsProcId rcvr=0; rcvr<gridSize; rcvr++)
  if (sndr!=rcvr)
  {
 
    CepsVector<CepsDofGlobalIndex> adj;
    CepsVector<CepsDofGlobalIndex> adjRows;
    CepsVector<CepsUInt>           adjPtr;
 
    CepsUInt adjSize     = 0u;
    CepsUInt adjRowsSize = 0u;
 
    CepsInt tag = 5*(gridSize*sndr+rcvr);
 
    if (rank==sndr)
    {
      for (auto& adjm: m_extraAdjacencyToSend[rcvr])
        adjSize += adjm.second.size();
      adjRowsSize = m_extraAdjacencyToSend[rcvr].size();
      MPI_Isend(&adjSize    ,1,CEPS_MPI_UINT,rcvr,0+tag,comm,&request);
      MPI_Isend(&adjRowsSize,1,CEPS_MPI_UINT,rcvr,1+tag,comm,&request);
    }
    else if(rank==rcvr)
    {
      MPI_Recv(&adjSize    ,1,CEPS_MPI_UINT,sndr,0+tag,comm,&status);
      MPI_Recv(&adjRowsSize,1,CEPS_MPI_UINT,sndr,1+tag,comm,&status);
    }
 
    adj    .resize(adjSize      );
    adjPtr .resize(adjRowsSize+1);
    adjRows.resize(adjRowsSize  );
 
    if (adjSize!= 0 and rank==sndr)
    {
      CepsUInt i     = 0;
      CepsUInt count = 0;
      for (auto& adjm: m_extraAdjacencyToSend[rcvr])
      {
        for (CepsDofGlobalIndex dofId: adjm.second)
        {
          adj[count] = dofId;
          count++;
        }
        adjRows[i  ] = adjm.first;
        adjPtr [i+1] = count;
        i++; 
      }
      MPI_Isend(adj    .data(),adjSize      ,CEPS_MPI_GLOBAL_INDEX,rcvr,2+tag,comm,&request);
      MPI_Isend(adjRows.data(),adjRowsSize  ,CEPS_MPI_GLOBAL_INDEX,rcvr,3+tag,comm,&request);
      MPI_Isend(adjPtr .data(),adjRowsSize+1,CEPS_MPI_UINT        ,rcvr,4+tag,comm,&request);
    }
    else if (adjSize!=0 and rank==rcvr)
    {
      MPI_Recv(adj    .data(),adjSize      ,CEPS_MPI_GLOBAL_INDEX,sndr,2+tag,comm,&status);
      MPI_Recv(adjRows.data(),adjRowsSize  ,CEPS_MPI_GLOBAL_INDEX,sndr,3+tag,comm,&status);
      MPI_Recv(adjPtr .data(),adjRowsSize+1,CEPS_MPI_UINT        ,sndr,4+tag,comm,&status);
 
      for (CepsUInt i=0;i<adjRowsSize;i++)
        for (CepsUInt ptr=adjPtr[i];ptr<adjPtr[i+1];ptr++)
          m_extraAdjacencyToRecv[adjRows[i]].insert(adj[ptr]);
    }
    ceps::barrier();
 
  }
 
}
 
void
FiniteElements::extractValuesForUnknown(
  const CepsUnknownIndex& uId,
  DVecPtr                 vec,
  CepsVector<CepsIndex>&  indices,
  CepsVector<CepsReal>&   values,
  CepsBool                sendToMaster,
  CepsBool                returnGeomIndices
)
{
 
  CepsProcId rank = ceps::getRank();
 
  // The selection criteria for dofs, based on arguments
  std::function<CepsBool(DegreeOfFreedom*)> selector = [&] (DegreeOfFreedom* dof) ->CepsBool {
    
    GeomNode* node = nullptr;
    if (returnGeomIndices)
      node = this->getFENode(dof->getDiscrId())->getGeomObject();
 
    // If extraction is for outputs, only give dofs defined on geometry
    return ((not returnGeomIndices) or node->getGlobalIndex()!=-1)
    // If we regroup on master, use only owned dofs, not halos
       and ((not sendToMaster) or dof->getOwner()==rank);
  };
 
  extractValuesForUnknownInternal(
    uId,vec,indices,values,sendToMaster,returnGeomIndices,selector
  );
  return;
}
 
CepsReal
getDomainMeasure(
  const CepsVector<FEBase*>&    cells,
  const CepsSet<CepsAttribute>& attributes,
  CepsBool                      onlyOfThisProc
)
{
  CEPS_ABORT_IF(attributes.empty(), "missing attributes");
 
  using selector_t = AttributesSelector<CepsVector<FEBase*>::const_iterator>;
  selector_t selec;
  selec.setAttributes(attributes.begin(), attributes.end());
  selec.onlyOnThisProc(true);
  selec.selectBetween(cells.begin(), cells.end());
  return getDomainMeasure(selec.getSelected(), onlyOfThisProc);
}
 
CepsReal
getDomainMeasure(const CepsVector<FEBase*>& cells, CepsBool onlyOfThisProc)
{
  CepsReal locMeas  = 0.0;
  CepsReal globMeas = 0.0;
 
  for (auto cell : cells)
    locMeas += cell->getGeomObject()->getMeasure();
 
  if (onlyOfThisProc)
    return locMeas;
 
  MPI_Allreduce(&locMeas, &globMeas, 1, CEPS_MPI_REAL, MPI_SUM, ceps::getCommunicator());
  return globMeas;
}