Mesh nodes numbering as implemented by FreeFem++

ConstructDataFElement_NodesNumbering.cpp

// This is a simplified version of the DoF numbering logic for a triangle mesh in FreeFem++
// * it's not dealing with periodic boundary conditions;
// * it's not dealing with the mortar method;
// It first numbers geometric features (verts, edges, tris) and build the DoF numbering from this intermediate numbering

ConstructDataFElement::ConstructDataFElement(const Mesh &Th, const  TypeOfFE &TFE) 
: // Number of tris
  m_NbOfElements(0),
  // Number of node per triangle: 3 if only verts nodes, 6 if verts + edge nodes ...
  m_NbNodePerElement(0),
  // Number of DoF per triangle: a node can have several DoF, for example P2 elems got 2 DoF per edge
  m_NbDFPerElement(0),
  // Array of node ids: for example if elems got 6 nodes 
  // m_NodesOfElement[6 * tid + nodeIndex] is the global id of the nodeIndexth node of triangle tid
  m_NodesOfElement(0),
  // m_FirstDfOfNode[nodeid] is the smaller DoF id for that node
  m_FirstDfOfNode(0)
{ 
  Make(Th,TFEs, tm);
}

ConstructDataFElement::~ConstructDataFElement()
{
  delete [] m_NodesOfElement;
  delete [] m_FirstDfOfNode;   
}

// This function will number nodes of mesh Th if they're used by finite element type TFE
void ConstructDataFElement::Make(const Mesh &Th, const  TypeOfFE &TFE)
{  
  assert(TFE.NbDoF == 3 * TFE.NbDfOnVertex + 3 * TFE.NbDfOnEdge + TFE.NbDfOnElement);
  
  // First thing we need to do is count how many DoF per node. For example:
  // P2 got 1 DoF per verts, 2 DoF per edge, 0 Dof per tri
  // We keep separate count for different verts, edges though.
  KN<int> NbDFonNode(TFE.NbNode);
  NbDFonNode = 0;
  for (int df = 0; df < TFE.NbDoF; ++df)
  {
    int node = TFE.NodeOfDF[df];    
    int ndfonnode = TFE.DFIdOnNode[df];
    // The largest local (to a node) DoF id tells us how many DoF there is for a node
    // There is no gab and local DoF ids start from 0, hence the +1
    NbDFonNode[node] = Max(NbDFonNode[node], ndfonnode + 1);
  }

  assert(NbDFonNode.sum() == TFE.NbDoF);
     
  m_NbOfElements = Th.nt;
  m_NbDFPerElement = TFE.NbDoF;  
  m_NbNodePerElement = TFE.NbNode;
  
  // Allocate and initialize the node index we're about to compute
  int totalNbNode =  Th.nt * TFE.NbNode;
  m_NodesOfElement = new int[totalNbNode];   
  for (int nodei = 0; nodei < totalNbNode; ++nodei)
  {
    m_NodesOfElement[nodei] = -1;
  }
 
  // We're going to number *all* verts first,
  // then we number T0 edge+elem nodes, T1 edge+elem nodes ... Tn-1 edge+elem nodes
  // but we need to cope with the case where we don't have verts nodes, for example
  // P0, where we only got a single elem node/DoF.
  int nbVertNodes = 0, edgeAndElemId = 0;
  if(TFE.HasVertexNode) 
  { 
    // So if we got verts nodes, we know we must pack 3 verts node ids before the other edges+elem nodes
    // Hence offset is 3
    nbVertNodes = 3;
    // and we weill number edges+elem nodes starting from the last verts node id
    edgeAndElemId =  Th.nv;
  }
        
  // We walk all the triangles, 
  for (int elemi = 0, nodesOfelemIdx = 0; elemi < Th.nt; ++elemi)
  {
    // verts nodes first
    if(TFE.HasVertexNode) 
    { 
      // Simply reuse the verts ids from the mesh
      for (int verti = 0; verti < 3; ++verti, ++nodesOfelemIdx)
      {
        m_NodesOfElement[nodesOfelemIdx] = Th(elemi, verti);
      }
    }

    // edge nodes next: we must not number a same edge twice if it's linked to two triangles (non-boundary edge)
    if(TFE.HasEdgeNode) 
    {
      for (int edgei = 0; edgei < 3; ++edgei, ++nodesOfelemIdx)
      {
        // If we haven't assigned an id to that edge yet
        if (m_NodesOfElement[nodesOfelemIdx] < 0) 
        {       
          // We lookup the other triangle connected to that edge (if any)
          int edgej = edgei;
          int elemj = Th.ElementAdj(elemi, edgej);
          assert(elemj >= elemi);
          // Compute the node offset within the index for the other triangle
          int nodeOfOtherElemIdx = elemj * TFE.NbNode + nbVertNodes + edgej;
          // And set the node id for both of these triangles in one shot
          m_NodesOfElement[nodesOfelemIdx] = m_NodesOfElement[nodeOfOtherElemIdx] = edgeAndElemId;
          // Id for the next edge or elem
          ++edgeAndElemId;
        }          
      }
    }

    // elem: one elem node per triangle, no sharing to worry about
    if (TFE.HasElementNode)    
    {
      m_NodesOfElement[nodesOfelemIdx] = edgeAndElemId;      
      ++edgeAndElemId;
      ++nodesOfelemIdx;
    }
  }                          

  // At that stage we numbered all the nodes: we must now number the DoFs.
  // In order to do so, we simply prefix sum a "number of DoFs per node" array, using the node mumbering we computed. 
  assert(edgeAndElemId == totalNbNode);

  // Allocate and initialize the array we will use for the prefix sum (hence the +1)
  m_FirstDfOfNode = new int [totalNbNode + 1];  
  for (int nodei = 0; nodei <= totalNbNode; ++nodei) 
  {
    m_FirstDfOfNode[nodei] = -1;
  }
 
  // m_FirstDfOfNode[nodei + 1] is the number of DoF for node nodei
  for (int nodesOfelemIdx = 0, elemi = 0; elemi < Th.nt; ++elemi)
  {
    for (int nodei = 0; nodei < TFE.NbNode; ++nodei, ++nodesOfelemIdx)
    {
      m_FirstDfOfNode[m_NodesOfElement[nodesOfelemIdx] + 1] = NbDFonNode[nodei];
    }
  }
   
  // Prefix sum: m_FirstDfOfNode[nodei] will contain the smallest DoF id assigned to node nodei
  m_FirstDfOfNode[0] = 0;
  for (int nodei = 0; nodei < totalNbNode; ++nodei) 
  {
    m_FirstDfOfNode[nodei + 1] += m_FirstDfOfNode[nodei];
  }    
}