MueLu_LocalLexicographicIndexManager_def.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu        (jhu@sandia.gov)
//                    Ray Tuminaro       (rstumin@sandia.gov)
//                    Luc Berger-Vergoat (lberge@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_LOCALLEXICOGRAPHICINDEXMANAGER_DEF_HPP_
#define MUELU_LOCALLEXICOGRAPHICINDEXMANAGER_DEF_HPP_

#include <MueLu_LocalLexicographicIndexManager_decl.hpp>
#include <Xpetra_MapFactory.hpp>

namespace MueLu {

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  LocalLexicographicIndexManager(const RCP<const Teuchos::Comm<int> > comm, const bool coupled,
                                 const int NumDimensions, const int interpolationOrder,
                                 const int MyRank, const int NumRanks,
                                 const Array<GO> GFineNodesPerDir, const Array<LO> LFineNodesPerDir,
                                 const Array<LO> CoarseRate, const Array<GO> MeshData) :
  IndexManager(comm, coupled, NumDimensions, interpolationOrder, GFineNodesPerDir, LFineNodesPerDir),
  myRank(MyRank), numRanks(NumRanks) {

    // Allocate data based on user input
    meshData.resize(numRanks);
    rankIndices.resize(numRanks);
    coarseMeshData.resize(numRanks);

    // Load coarse rate, being careful about formating
    for(int dim = 0; dim < 3; ++dim) {
      if(dim < this->numDimensions) {
        if(CoarseRate.size() == 1) {
          this->coarseRate[dim] = CoarseRate[0];
        } else if(CoarseRate.size() == this->numDimensions) {
          this->coarseRate[dim] = CoarseRate[dim];
        }
      } else {
        this->coarseRate[dim] = 1;
      }
    }

    // Load meshData for local lexicographic case
    for(int rank = 0; rank < numRanks; ++rank) {
      meshData[rank].resize(10);
      for(int entry = 0; entry < 10; ++entry) {
        meshData[rank][entry] = MeshData[10*rank + entry];
      }
    }

    if(this->coupled_) {
      myBlock = meshData[myRank][2];
      sortLocalLexicographicData();
    }

    // Start simple parameter calculation
    myRankIndex = rankIndices[myRank];
    for(int dim = 0; dim < 3; ++dim) {
      this->startIndices[dim]     = meshData[myRankIndex][2*dim + 3];
      this->startIndices[dim + 3] = meshData[myRankIndex][2*dim + 4];
    }

    this->computeMeshParameters();
    computeGlobalCoarseParameters();
    computeCoarseLocalLexicographicData();
  } // Constructor

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  computeGlobalCoarseParameters() {
    this->gNumCoarseNodes10 = this->gCoarseNodesPerDir[0]*this->gCoarseNodesPerDir[1];
    this->gNumCoarseNodes   = this->gNumCoarseNodes10*this->gCoarseNodesPerDir[2];
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getGhostedNodesData(const RCP<const Map>/* fineMap */,
                      Array<LO>& ghostedNodeCoarseLIDs,
                      Array<int>& ghostedNodeCoarsePIDs,
                      Array<GO>& ghostedNodeCoarseGIDs) const {

    // First we allocated memory for the outputs
    ghostedNodeCoarseLIDs.resize(this->getNumLocalGhostedNodes());
    ghostedNodeCoarsePIDs.resize(this->getNumLocalGhostedNodes());
    ghostedNodeCoarseGIDs.resize(this->numGhostedNodes);

    // Now the tricky part starts, the coarse nodes / ghosted coarse nodes need to be imported.
    // This requires finding what their GID on the fine mesh is. They need to be ordered
    // lexicographically to allow for fast sweeps through the mesh.

    // We loop over all ghosted coarse nodes by increasing global lexicographic order
    Array<LO>  ghostedCoarseNodeCoarseIndices(3), ghostedCoarseNodeFineIndices(3);
    Array<LO>  lCoarseNodeCoarseIndices(3);
    Array<GO>  lCoarseNodeCoarseGIDs(this->lNumCoarseNodes);
    LO currentIndex = -1, countCoarseNodes = 0;
    for(int k = 0; k < this->ghostedNodesPerDir[2]; ++k) {
      for(int j = 0; j < this->ghostedNodesPerDir[1]; ++j) {
        for(int i = 0; i < this->ghostedNodesPerDir[0]; ++i) {
          currentIndex = k*this->numGhostedNodes10 + j*this->ghostedNodesPerDir[0] + i;
          ghostedCoarseNodeCoarseIndices[0] = this->startGhostedCoarseNode[0] + i;
          ghostedCoarseNodeFineIndices[0] = ghostedCoarseNodeCoarseIndices[0]*this->coarseRate[0];
          if(ghostedCoarseNodeFineIndices[0] > this->gFineNodesPerDir[0] - 1) {
            ghostedCoarseNodeFineIndices[0] = this->gFineNodesPerDir[0] - 1;
          }
          ghostedCoarseNodeCoarseIndices[1] = this->startGhostedCoarseNode[1] + j;
          ghostedCoarseNodeFineIndices[1] = ghostedCoarseNodeCoarseIndices[1]*this->coarseRate[1];
          if(ghostedCoarseNodeFineIndices[1] > this->gFineNodesPerDir[1] - 1) {
            ghostedCoarseNodeFineIndices[1] = this->gFineNodesPerDir[1] - 1;
          }
          ghostedCoarseNodeCoarseIndices[2] = this->startGhostedCoarseNode[2] + k;
          ghostedCoarseNodeFineIndices[2] = ghostedCoarseNodeCoarseIndices[2]*this->coarseRate[2];
          if(ghostedCoarseNodeFineIndices[2] > this->gFineNodesPerDir[2] - 1) {
            ghostedCoarseNodeFineIndices[2] = this->gFineNodesPerDir[2] - 1;
          }

          GO myGID = -1, myCoarseGID = -1;
          LO myLID = -1, myPID = -1, myCoarseLID = -1;
          getGIDLocalLexicographic(i, j, k, ghostedCoarseNodeFineIndices, myGID, myPID, myLID);

          int rankIndex = rankIndices[myPID];
          for(int dim = 0; dim < 3; ++dim) {
            if(dim < this->numDimensions) {
              lCoarseNodeCoarseIndices[dim] = ghostedCoarseNodeCoarseIndices[dim]
                - coarseMeshData[rankIndex][3 + 2*dim];
            }
          }
          LO myRankIndexCoarseNodesInDir0 = coarseMeshData[rankIndex][4]
            - coarseMeshData[rankIndex][3] + 1;
          LO myRankIndexCoarseNodes10 = (coarseMeshData[rankIndex][6]
                                         - coarseMeshData[rankIndex][5] + 1)
            *myRankIndexCoarseNodesInDir0;
          myCoarseLID = lCoarseNodeCoarseIndices[2]*myRankIndexCoarseNodes10
            + lCoarseNodeCoarseIndices[1]*myRankIndexCoarseNodesInDir0
            + lCoarseNodeCoarseIndices[0];
          myCoarseGID = myCoarseLID + coarseMeshData[rankIndex][9];

          ghostedNodeCoarseLIDs[currentIndex] = myCoarseLID;
          ghostedNodeCoarsePIDs[currentIndex] = myPID;
          ghostedNodeCoarseGIDs[currentIndex] = myCoarseGID;

          if(myPID == myRank) {
            lCoarseNodeCoarseGIDs[countCoarseNodes] = myCoarseGID;
            ++countCoarseNodes;
          }
        }
      }
    }
  }

  template<class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodesData(const RCP<const Map> fineCoordinatesMap,
                     Array<GO>& coarseNodeCoarseGIDs,
                     Array<GO>& coarseNodeFineGIDs) const {

    // Allocate sufficient storage space for outputs
    coarseNodeCoarseGIDs.resize(this->getNumLocalCoarseNodes());
    coarseNodeFineGIDs.resize(this->getNumLocalCoarseNodes());

    // Load all the GIDs on the fine mesh
    ArrayView<const GO> fineNodeGIDs = fineCoordinatesMap->getNodeElementList();

    Array<GO> coarseStartIndices(3);
    for(int dim = 0; dim < 3; ++dim) {
      coarseStartIndices[dim] = this->coarseMeshData[myRankIndex][2*dim + 3];
    }

    // Extract the fine LIDs of the coarse nodes and store the corresponding GIDs
    LO fineLID;
    for(LO coarseLID = 0; coarseLID < this->getNumLocalCoarseNodes(); ++coarseLID) {
      Array<LO> coarseIndices(3), fineIndices(3), gCoarseIndices(3);
      this->getCoarseNodeLocalTuple(coarseLID,
                                    coarseIndices[0],
                                    coarseIndices[1],
                                    coarseIndices[2]);
      getCoarseNodeFineLID(coarseIndices[0],coarseIndices[1],coarseIndices[2],fineLID);
      coarseNodeFineGIDs[coarseLID] = fineNodeGIDs[fineLID];

      LO myRankIndexCoarseNodesInDir0 = coarseMeshData[myRankIndex][4]
            - coarseMeshData[myRankIndex][3] + 1;
      LO myRankIndexCoarseNodes10 = (coarseMeshData[myRankIndex][6]
                                  - coarseMeshData[myRankIndex][5] + 1)
                                  *myRankIndexCoarseNodesInDir0;
      LO myCoarseLID = coarseIndices[2]*myRankIndexCoarseNodes10
                     + coarseIndices[1]*myRankIndexCoarseNodesInDir0
                     + coarseIndices[0];
      GO myCoarseGID = myCoarseLID + coarseMeshData[myRankIndex][9];
      coarseNodeCoarseGIDs[coarseLID] = myCoarseGID;
    }

  }

  template<class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getGIDLocalLexicographic(const LO iGhosted, const LO jGhosted, const LO kGhosted,
                           const Array<LO> coarseNodeFineIndices,
                           GO& myGID, LO& myPID, LO& myLID) const {

    LO ni = -1, nj = -1, li = -1, lj = -1, lk = -1;
    LO myRankGuess = myRankIndex;
    // We try to make a logical guess as to which PID owns the current coarse node
    if(iGhosted == 0 && this->ghostInterface[0]) {
      --myRankGuess;
    } else if((iGhosted == this->ghostedNodesPerDir[0] - 1) && this->ghostInterface[1]) {
      ++myRankGuess;
    }
    if(jGhosted == 0 && this->ghostInterface[2]) {
      myRankGuess -= pi;
    } else if((jGhosted == this->ghostedNodesPerDir[1] - 1) && this->ghostInterface[3]) {
      myRankGuess += pi;
    }
    if(kGhosted == 0 && this->ghostInterface[4]) {
      myRankGuess -= pj*pi;
    } else if((kGhosted == this->ghostedNodesPerDir[2] - 1) && this->ghostInterface[5]) {
      myRankGuess += pj*pi;
    }
    if(coarseNodeFineIndices[0] >= meshData[myRankGuess][3]
       && coarseNodeFineIndices[0] <= meshData[myRankGuess][4]
       && coarseNodeFineIndices[1] >= meshData[myRankGuess][5]
       && coarseNodeFineIndices[1] <= meshData[myRankGuess][6]
       && coarseNodeFineIndices[2] >= meshData[myRankGuess][7]
       && coarseNodeFineIndices[2] <= meshData[myRankGuess][8]
       && myRankGuess < numRanks - 1) {
      myPID = meshData[myRankGuess][0];
      ni = meshData[myRankGuess][4] - meshData[myRankGuess][3] + 1;
      nj = meshData[myRankGuess][6] - meshData[myRankGuess][5] + 1;
      li = coarseNodeFineIndices[0] - meshData[myRankGuess][3];
      lj = coarseNodeFineIndices[1] - meshData[myRankGuess][5];
      lk = coarseNodeFineIndices[2] - meshData[myRankGuess][7];
      myLID = lk*nj*ni + lj*ni + li;
      myGID = meshData[myRankGuess][9] + myLID;
    } else { // The guess failed, let us use the heavy artilery: std::find_if()
      // It could be interesting to monitor how many times this branch of the code gets
      // used as it is far more expensive than the above one...
      auto nodeRank = std::find_if(myBlockStart, myBlockEnd,
                                   [coarseNodeFineIndices](const std::vector<GO>& vec){
                                     if(coarseNodeFineIndices[0] >= vec[3]
                                        && coarseNodeFineIndices[0] <= vec[4]
                                        && coarseNodeFineIndices[1] >= vec[5]
                                        && coarseNodeFineIndices[1] <= vec[6]
                                        && coarseNodeFineIndices[2] >= vec[7]
                                        && coarseNodeFineIndices[2] <= vec[8]) {
                                       return true;
                                     } else {
                                       return false;
                                     }
                                   });
      myPID = (*nodeRank)[0];
      ni = (*nodeRank)[4] - (*nodeRank)[3] + 1;
      nj = (*nodeRank)[6] - (*nodeRank)[5] + 1;
      li = coarseNodeFineIndices[0] - (*nodeRank)[3];
      lj = coarseNodeFineIndices[1] - (*nodeRank)[5];
      lk = coarseNodeFineIndices[2] - (*nodeRank)[7];
      myLID = lk*nj*ni + lj*ni + li;
      myGID = (*nodeRank)[9] + myLID;
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  sortLocalLexicographicData() {

    std::sort(meshData.begin(), meshData.end(),
              [](const std::vector<GO>& a, const std::vector<GO>& b)->bool {
                // The below function sorts ranks by blockID, kmin, jmin and imin
                if(a[2] < b[2]) {
                  return true;
                } else if(a[2] == b[2]) {
                  if(a[7] < b[7]) {
                    return true;
                  } else if(a[7] == b[7]) {
                    if(a[5] < b[5]) {
                      return true;
                    } else if(a[5] == b[5]) {
                      if(a[3] < b[3]) {return true;}
                    }
                  }
                }
                return false;
              });

    numBlocks = meshData[numRanks - 1][2] + 1;
    // Find the range of the current block
    myBlockStart = std::lower_bound(meshData.begin(), meshData.end(), myBlock - 1,
                                    [] (const std::vector<GO>& vec, const GO val)->bool {
                                      return (vec[2] < val) ? true : false;
                                    });
    myBlockEnd = std::upper_bound(meshData.begin(), meshData.end(), myBlock,
                                  [] (const GO val, const std::vector<GO>& vec)->bool {
                                    return (val < vec[2]) ? true : false;
                                  });
    // Assuming that i,j,k and ranges are split in pi, pj and pk processors
    // we search for these numbers as they will allow us to find quickly the PID of processors
    // owning ghost nodes.
    auto myKEnd = std::upper_bound(myBlockStart, myBlockEnd, (*myBlockStart)[3],
                                   [] (const GO val, const std::vector<GO>& vec)->bool {
                                     return (val < vec[7]) ? true : false;
                                   });
    auto myJEnd = std::upper_bound(myBlockStart, myKEnd, (*myBlockStart)[3],
                                   [] (const GO val, const std::vector<GO>& vec)->bool {
                                     return (val < vec[5]) ? true : false;
                                   });
    pi = std::distance(myBlockStart, myJEnd);
    pj = std::distance(myBlockStart, myKEnd) / pi;
    pk = std::distance(myBlockStart, myBlockEnd) / (pj*pi);

    // We also look for the index of the local rank in the current block.
    const int MyRank = myRank;
    myRankIndex = std::distance(meshData.begin(),
                                std::find_if(myBlockStart, myBlockEnd,
                                             [MyRank] (const std::vector<GO>& vec)->bool {
                                               return (vec[0] == MyRank) ? true : false;
                                             })
                                );
    // We also construct a mapping of rank to rankIndex in the meshData vector,
    // this will allow us to access data quickly later on.
    for(int rankIndex = 0; rankIndex < numRanks; ++rankIndex) {
      rankIndices[meshData[rankIndex][0]] = rankIndex;
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  computeCoarseLocalLexicographicData() {
    Array<LO> rankOffset(3);
    for(int rank = 0; rank < numRanks; ++rank) {
      coarseMeshData[rank].resize(10);
      coarseMeshData[rank][0] = meshData[rank][0];
      coarseMeshData[rank][1] = meshData[rank][1];
      coarseMeshData[rank][2] = meshData[rank][2];
      for(int dim = 0; dim < 3; ++dim) {
        coarseMeshData[rank][3 + 2*dim] = meshData[rank][3 + 2*dim] / this->coarseRate[dim];
        if(meshData[rank][3 + 2*dim] % this->coarseRate[dim] > 0) {
          ++coarseMeshData[rank][3 + 2*dim];
        }
        coarseMeshData[rank][3 + 2*dim + 1] = meshData[rank][3 + 2*dim + 1] / this->coarseRate[dim];
        if(meshData[rank][3 + 2*dim + 1] == this->gFineNodesPerDir[dim] - 1 &&
           meshData[rank][3 + 2*dim + 1] % this->coarseRate[dim] > 0) {
           //this->endRate[dim] < this->coarseRate[dim]) {
          ++coarseMeshData[rank][3 + 2*dim + 1];
        }
      }
      if(rank > 0) {
        coarseMeshData[rank][9] = coarseMeshData[rank - 1][9]
          + (coarseMeshData[rank - 1][8] - coarseMeshData[rank - 1][7] + 1)
          * (coarseMeshData[rank - 1][6] - coarseMeshData[rank - 1][5] + 1)
          * (coarseMeshData[rank - 1][4] - coarseMeshData[rank - 1][3] + 1);
      }
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  std::vector<std::vector<GlobalOrdinal> > LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseMeshData() const {return coarseMeshData;}

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getFineNodeGlobalTuple(const GO /* myGID */, GO& /* i */, GO& /* j */, GO& /* k */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getFineNodeLocalTuple(const LO myLID, LO& i, LO& j, LO& k) const {
    LO tmp;
    k   = myLID / this->lNumFineNodes10;
    tmp = myLID % this->lNumFineNodes10;
    j   = tmp   / this->lFineNodesPerDir[0];
    i   = tmp   % this->lFineNodesPerDir[0];
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getFineNodeGhostedTuple(const LO myLID, LO& i, LO& j, LO& k) const {
    LO tmp;
    k   = myLID / this->lNumFineNodes10;
    tmp = myLID % this->lNumFineNodes10;
    j   = tmp   / this->lFineNodesPerDir[0];
    i   = tmp   % this->lFineNodesPerDir[0];

    k += this->offsets[2];
    j += this->offsets[1];
    i += this->offsets[0];
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getFineNodeGID(const GO /* i */, const GO /* j */, const GO /* k */, GO& /* myGID */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getFineNodeLID(const LO /* i */, const LO /* j */, const LO /* k */, LO& /* myLID */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeGlobalTuple(const GO /* myGID */, GO& /* i */, GO& /* j */, GO& /* k */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeLocalTuple(const LO myLID, LO& i, LO& j, LO& k) const {
    LO tmp;
    k   = myLID / this->lNumCoarseNodes10;
    tmp = myLID % this->lNumCoarseNodes10;
    j   = tmp / this->lCoarseNodesPerDir[0];
    i   = tmp % this->lCoarseNodesPerDir[0];
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeGID(const GO /* i */, const GO /* j */, const GO /* k */, GO& /* myGID */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeLID(const LO /* i */, const LO /* j */, const LO /* k */, LO& /* myLID */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeGhostedLID(const LO i, const LO j, const LO k, LO& myLID) const {
    myLID = k*this->numGhostedNodes10 + j*this->ghostedNodesPerDir[0] + i;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getCoarseNodeFineLID(const LO i, const LO j, const LO k, LO& myLID) const {
    // Assumptions: (i,j,k) is a tuple on the coarse mesh
    //              myLID is the corresponding local ID on the fine mesh
    const LO multiplier[3] = {1, this->lFineNodesPerDir[0], this->lNumFineNodes10};
    const LO indices[3] = {i, j, k};

    myLID = 0;
    for(int dim = 0; dim < 3; ++dim) {
      if((indices[dim] == this->getLocalCoarseNodesInDir(dim) - 1) && this->meshEdge[2*dim + 1]) {
        // We are dealing with the last node on the mesh in direction dim
        // so we can simply use the number of nodes on the fine mesh in that direction
        myLID += (this->getLocalFineNodesInDir(dim) - 1)*multiplier[dim];
      } else {
        myLID += (indices[dim]*this->getCoarseningRate(dim) + this->getCoarseNodeOffset(dim))
          *multiplier[dim];
      }
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getGhostedNodeFineLID(const LO /* i */, const LO /* j */, const LO /* k */, LO& /* myLID */) const {
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void LocalLexicographicIndexManager<LocalOrdinal, GlobalOrdinal, Node>::
  getGhostedNodeCoarseLID(const LO /* i */, const LO /* j */, const LO /* k */, LO& /* myLID */) const {
  }

} //namespace MueLu

#endif /* MUELU_LOCALLEXICOGRAPHICINDEXMANAGER_DEF_HPP_ */