MueLu_IndexManager_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-Vergiat (lberge@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_INDEXMANAGER_DEF_HPP
#define MUELU_INDEXMANAGER_DEF_HPP

#include "Teuchos_OrdinalTraits.hpp"

#include "MueLu_ConfigDefs.hpp"
#include <MueLu_IndexManager_decl.hpp>

/*****************************************************************************

****************************************************************************/

namespace MueLu {

template <class LocalOrdinal, class GlobalOrdinal, class Node>
IndexManager<LocalOrdinal, GlobalOrdinal, Node>::
    IndexManager(const RCP<const Teuchos::Comm<int> > comm,
                 const bool coupled,
                 const bool singleCoarsePoint,
                 const int NumDimensions,
                 const int interpolationOrder,
                 const Array<GO> GFineNodesPerDir,
                 const Array<LO> LFineNodesPerDir)
  : comm_(comm)
  , coupled_(coupled)
  , singleCoarsePoint_(singleCoarsePoint)
  , numDimensions(NumDimensions)
  , interpolationOrder_(interpolationOrder)
  , gFineNodesPerDir(GFineNodesPerDir)
  , lFineNodesPerDir(LFineNodesPerDir) {
  coarseRate.resize(3);
  endRate.resize(3);
  gCoarseNodesPerDir.resize(3);
  lCoarseNodesPerDir.resize(3);
  ghostedNodesPerDir.resize(3);

  offsets.resize(3);
  coarseNodeOffsets.resize(3);
  startIndices.resize(6);
  startGhostedCoarseNode.resize(3);

}  // Constructor

template <class LocalOrdinal, class GlobalOrdinal, class Node>
void IndexManager<LocalOrdinal, GlobalOrdinal, Node>::
    computeMeshParameters() {
  RCP<Teuchos::FancyOStream> out;
  if (const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
    out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setShowAllFrontMatter(false).setShowProcRank(true);
  } else {
    out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
  }

  if (coupled_) {
    gNumFineNodes10 = gFineNodesPerDir[1] * gFineNodesPerDir[0];
    gNumFineNodes   = gFineNodesPerDir[2] * gNumFineNodes10;
  } else {
    gNumFineNodes10 = Teuchos::OrdinalTraits<GO>::invalid();
    gNumFineNodes   = Teuchos::OrdinalTraits<GO>::invalid();
  }
  lNumFineNodes10 = lFineNodesPerDir[1] * lFineNodesPerDir[0];
  lNumFineNodes   = lFineNodesPerDir[2] * lNumFineNodes10;
  for (int dim = 0; dim < 3; ++dim) {
    if (dim < numDimensions) {
      if (coupled_) {
        if (startIndices[dim] == 0) {
          meshEdge[2 * dim] = true;
        }
        if (startIndices[dim + 3] + 1 == gFineNodesPerDir[dim]) {
          meshEdge[2 * dim + 1] = true;
          endRate[dim]          = startIndices[dim + 3] % coarseRate[dim];
        }
      } else {  // With uncoupled problem each rank might require a different endRate
        meshEdge[2 * dim]     = true;
        meshEdge[2 * dim + 1] = true;
        endRate[dim]          = (lFineNodesPerDir[dim] - 1) % coarseRate[dim];
      }
      if (endRate[dim] == 0) {
        endRate[dim] = coarseRate[dim];
      }

      // If uncoupled aggregation is used, offsets[dim] = 0, so nothing to do.
      if (coupled_) {
        offsets[dim] = Teuchos::as<LO>(startIndices[dim]) % coarseRate[dim];
        if (offsets[dim] == 0) {
          coarseNodeOffsets[dim] = 0;
        } else if (startIndices[dim] + endRate[dim] == lFineNodesPerDir[dim]) {
          coarseNodeOffsets[dim] = endRate[dim] - offsets[dim];
        } else {
          coarseNodeOffsets[dim] = coarseRate[dim] - offsets[dim];
        }

        if (interpolationOrder_ == 0) {
          int rem = startIndices[dim] % coarseRate[dim];
          if ((rem != 0) && (rem <= Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
            ghostInterface[2 * dim] = true;
          }
          rem = startIndices[dim + 3] % coarseRate[dim];
          // uncoupled by nature does not require ghosts nodes
          if (coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
              (rem > Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
            ghostInterface[2 * dim + 1] = true;
          }

        } else if (interpolationOrder_ == 1) {
          if (coupled_ && (startIndices[dim] % coarseRate[dim] != 0 ||
                           startIndices[dim] == gFineNodesPerDir[dim] - 1)) {
            ghostInterface[2 * dim] = true;
          }
          if (coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
              ((lFineNodesPerDir[dim] == 1) || (startIndices[dim + 3] % coarseRate[dim] != 0))) {
            ghostInterface[2 * dim + 1] = true;
          }
        }
      }
    } else {  // Default value for dim >= numDimensions
      endRate[dim] = 1;
    }
  }

  *out << "singleCoarsePoint? " << singleCoarsePoint_ << std::endl;
  *out << "gFineNodesPerDir: " << gFineNodesPerDir << std::endl;
  *out << "lFineNodesPerDir: " << lFineNodesPerDir << std::endl;
  *out << "endRate: " << endRate << std::endl;
  *out << "ghostInterface: {" << ghostInterface[0] << ", " << ghostInterface[1] << ", "
       << ghostInterface[2] << ", " << ghostInterface[3] << ", " << ghostInterface[4] << ", "
       << ghostInterface[5] << "}" << std::endl;
  *out << "meshEdge: {" << meshEdge[0] << ", " << meshEdge[1] << ", "
       << meshEdge[2] << ", " << meshEdge[3] << ", " << meshEdge[4] << ", "
       << meshEdge[5] << "}" << std::endl;
  *out << "startIndices: " << startIndices << std::endl;
  *out << "offsets: " << offsets << std::endl;
  *out << "coarseNodeOffsets: " << coarseNodeOffsets << std::endl;

  // Here one element can represent either the degenerate case of one node or the more general
  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
  // one node. This helps generating a 3D space from tensorial products...
  // A good way to handle this would be to generalize the algorithm to take into account the
  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
  // discretization will have a unique node per element. This way 1D discretization can be
  // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
  // element in the z direction.
  // !!! Operations below are aftecting both local and global values that have two         !!!
  // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
  // coarseRate, endRate and offsets are in the global basis, as well as all the variables
  // starting with a g.
  // !!! while the variables starting with an l are in the local basis.                    !!!
  for (int dim = 0; dim < 3; ++dim) {
    if (dim < numDimensions) {
      // Check whether the partition includes the "end" of the mesh which means that endRate
      // will apply. Also make sure that endRate is not 0 which means that the mesh does not
      // require a particular treatment at the boundaries.
      if (meshEdge[2 * dim + 1]) {
        lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] - endRate[dim] + offsets[dim] - 1) / coarseRate[dim] + 1;
        if (offsets[dim] == 0) {
          ++lCoarseNodesPerDir[dim];
        }
        // We might want to coarsening the direction
        // into a single layer if there are not enough
        // points left to form two aggregates
        if (singleCoarsePoint_ && lFineNodesPerDir[dim] - 1 < coarseRate[dim]) {
          lCoarseNodesPerDir[dim] = 1;
        }
      } else {
        lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] + offsets[dim] - 1) / coarseRate[dim];
        if (offsets[dim] == 0) {
          ++lCoarseNodesPerDir[dim];
        }
      }

      // The first branch of this if-statement will be used if the rank contains only one layer
      // of nodes in direction i, that layer must also coincide with the boundary of the mesh
      // and coarseRate[i] == endRate[i]...
      if (interpolationOrder_ == 0) {
        startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
        int rem                     = startIndices[dim] % coarseRate[dim];
        if (rem > (Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
          ++startGhostedCoarseNode[dim];
        }
      } else {
        if ((startIndices[dim] == gFineNodesPerDir[dim] - 1) &&
            (startIndices[dim] % coarseRate[dim] == 0)) {
          startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim] - 1;
        } else {
          startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
        }
      }

      // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next
      // level.
      gCoarseNodesPerDir[dim] = (gFineNodesPerDir[dim] - 1) / coarseRate[dim];
      if ((gFineNodesPerDir[dim] - 1) % coarseRate[dim] == 0) {
        ++gCoarseNodesPerDir[dim];
      } else {
        gCoarseNodesPerDir[dim] += 2;
      }
    } else {  // Default value for dim >= numDimensions
      // endRate[dim] = 1;
      gCoarseNodesPerDir[dim] = 1;
      lCoarseNodesPerDir[dim] = 1;
    }  // if (dim < numDimensions)

    // This would happen if the rank does not own any nodes but in that case a subcommunicator
    // should be used so this should really not be a concern.
    if (lFineNodesPerDir[dim] < 1) {
      lCoarseNodesPerDir[dim] = 0;
    }
    ghostedNodesPerDir[dim] = lCoarseNodesPerDir[dim];
    // Check whether face *low needs ghost nodes
    if (ghostInterface[2 * dim]) {
      ghostedNodesPerDir[dim] += 1;
    }
    // Check whether face *hi needs ghost nodes
    if (ghostInterface[2 * dim + 1]) {
      ghostedNodesPerDir[dim] += 1;
    }
  }  // Loop for dim=0:3

  // With uncoupled aggregation we need to communicate to compute the global number of coarse points
  if (!coupled_) {
    for (int dim = 0; dim < 3; ++dim) {
      gCoarseNodesPerDir[dim] = -1;
    }
  }

  // Compute cummulative values
  lNumCoarseNodes10 = lCoarseNodesPerDir[0] * lCoarseNodesPerDir[1];
  lNumCoarseNodes   = lNumCoarseNodes10 * lCoarseNodesPerDir[2];
  numGhostedNodes10 = ghostedNodesPerDir[1] * ghostedNodesPerDir[0];
  numGhostedNodes   = numGhostedNodes10 * ghostedNodesPerDir[2];
  numGhostNodes     = numGhostedNodes - lNumCoarseNodes;

  *out << "lCoarseNodesPerDir: " << lCoarseNodesPerDir << std::endl;
  *out << "gCoarseNodesPerDir: " << gCoarseNodesPerDir << std::endl;
  *out << "ghostedNodesPerDir: " << ghostedNodesPerDir << std::endl;
  *out << "lNumCoarseNodes=" << lNumCoarseNodes << std::endl;
  *out << "numGhostedNodes=" << numGhostedNodes << std::endl;
}

}  // namespace MueLu

#define MUELU_INDEXMANAGER_SHORT
#endif  // MUELU_INDEXMANAGER_DEF_HPP