ParallelDofMapping.cc

//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "parallel/ParallelDofMapping.h"
#include "parallel/StdMpi.h"

namespace AMDiS {

  using namespace std;

  void FeSpaceDofMap::clear()
  {
    int nLevel = levelData->getLevelNumber();
    dofMap.clear();
    dofMap.resize(nLevel);

    nonRankDofs.clear();
    nonRankDofs.resize(nLevel);

    for (int i = 0; i < nLevel; i++) {
      nRankDofs[i] = 0;
      nLocalDofs[i] = 0;
      nOverallDofs[i] = 0;
      rStartDofs[i] = 0;
    }
  }


  void FeSpaceDofMap::update()
  {
    FUNCNAME("FeSpaceDofMap::update()");

    int nLevel = levelData->getLevelNumber();
    for (int i = 0; i < nLevel; i++) {

      // === Compute local indices for all rank owned DOFs. ===
      
      for (DofMap::iterator it = dofMap[i].begin(); it != dofMap[i].end(); ++it)
	if (it->second.local == -1 && nonRankDofs[i].count(it->first) == 0)
	  it->second.local = nRankDofs[i]++;
      
      // === Compute number of local and global DOFs in the mapping. ===
      
      nOverallDofs[i] = 0;
      rStartDofs[i] = 0;
      mpi::getDofNumbering(levelData->getMpiComm(i), 
			   nRankDofs[i], rStartDofs[i], nOverallDofs[i]);
      
      // === If required, compute also the global indices. ===
      
      if (needGlobalMapping) {
	computeGlobalMapping(i);
	
	if (hasNonLocalDofs)
	  computeNonLocalIndices(i);
      }
    }
  }


  void FeSpaceDofMap::computeGlobalMapping(int level)
  {
    FUNCNAME("FeSpaceDofMap::computeGlobalMapping()");

    for (DofMap::iterator it = dofMap[level].begin(); 
	 it != dofMap[level].end(); ++it)
      it->second.global = it->second.local + rStartDofs[level];
  }


  void FeSpaceDofMap::computeNonLocalIndices(int level)
  {
    FUNCNAME("FeSpaceDofMap::computeNonLocalIndices()");

    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

    // === Send all global indices of DOFs that are owned by the rank to all ===
    // === other ranks that also include this DOF.                           ===

    StdMpi<vector<int> > stdMpi(levelData->getMpiComm(level));

    for (DofComm::Iterator it(*sendDofs, level, feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	if (dofMap[level].count(it.getDofIndex()) && !nonRankDofs[level].count(it.getDofIndex()))
	  stdMpi.getSendData(it.getRank()).push_back(dofMap[level][it.getDofIndex()].global);

    stdMpi.updateSendDataSize();


    // === Check from which ranks this rank must receive some data. ===

    for (DofComm::Iterator it(*recvDofs, level, feSpace); 
	 !it.end(); it.nextRank()) {
      bool recvFromRank = false;
      for (; !it.endDofIter(); it.nextDof()) {
	if (nonRankDofs[level].count(it.getDofIndex())) {
	  recvFromRank = true;
	  break;
	}
      }

      if (recvFromRank)
	stdMpi.recv(it.getRank());
    }


    // === Start communication to exchange global indices. ===

    stdMpi.startCommunication();


    // === And set the global indices for all DOFs that are not owned by rank. ===
    
    for (DofComm::Iterator it(*recvDofs, level, feSpace);
	 !it.end(); it.nextRank()) {
      int i = 0;
      for (; !it.endDofIter(); it.nextDof())
	if (nonRankDofs[level].count(it.getDofIndex()))
	  dofMap[level][it.getDofIndex()].global = stdMpi.getRecvData(it.getRank())[i++];
    }
  }


  void ParallelDofMapping::init(MeshLevelData &ldata,
				vector<const FiniteElemSpace*> &fe,
				vector<const FiniteElemSpace*> &uniqueFe,
				bool needGlobalMapping,
				bool bNonLocalDofs)
  {
    FUNCNAME("ParallelDofMapping::init()");

    levelData = &ldata;
    feSpaces = fe;
    feSpacesUnique = uniqueFe;
    hasNonLocalDofs = bNonLocalDofs;
    

    // === Init the mapping for all different FE spaces. ===

    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	 it != feSpacesUnique.end(); ++it) {
      addFeSpace(*it);
      data[*it].setNeedGlobalMapping(needGlobalMapping);
      data[*it].setNonLocalDofs(hasNonLocalDofs);
    }
  }


  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	 it != feSpacesUnique.end(); ++it)
      data[*it].clear();

    int nLevel = levelData->getLevelNumber();
    nRankDofs.resize(nLevel);
    nLocalDofs.resize(nLevel);
    nOverallDofs.resize(nLevel);
    rStartDofs.resize(nLevel);

    for (int i = 0; i < nLevel; i++) {
      nRankDofs[i] = -1;
      nLocalDofs[i] = -1;
      nOverallDofs[i] = -1;
      rStartDofs[i] = -1;
    }

    dofToMatIndex.clear();
  }


  void ParallelDofMapping::setDofComm(DofComm &pSend, DofComm &pRecv)
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

    sendDofs = &pSend;
    recvDofs = &pRecv;

    // Add the DOF communicator also to all FE space DOF mappings.
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	 it != feSpacesUnique.end(); ++it)
      data[*it].setDofComm(pSend, pRecv);
  }


  void ParallelDofMapping::addFeSpace(const FiniteElemSpace* feSpace)
  {
    FUNCNAME("ParallelDofMapping::addFeSpace()");
    
    if (data.count(feSpace))
      data.find(feSpace)->second.clear();
    else
      data.insert(make_pair(feSpace, FeSpaceDofMap(levelData)));
    
    data.find(feSpace)->second.setFeSpace(feSpace);
  }    


  int ParallelDofMapping::computeRankDofs(int level)
  {
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
    
    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
      result += data[feSpaces[i]].nRankDofs[level];
    }
    
    return result;
  }


  int ParallelDofMapping::computeLocalDofs(int level)
  {
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
    
    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
      result += data[feSpaces[i]].nLocalDofs[level];
    }
    
    return result;
  }


  int ParallelDofMapping::computeOverallDofs(int level)
  {
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");

    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
      result += data.find(feSpaces[i])->second.nOverallDofs[level];
    }
    
    return result;
  }


  int ParallelDofMapping::computeStartDofs(int level)
  {
    FUNCNAME("ParallelDofMapping::computeStartDofs()");

    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
      result += data.find(feSpaces[i])->second.rStartDofs[level];
    }
    
    return result;
  }


  void ParallelDofMapping::update()
  {
    FUNCNAME("ParallelDofMapping::update()");

    // First, update all FE space DOF mappings.
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	 it != feSpacesUnique.end(); ++it)
      data[*it].update();
    
    int nLevel = levelData->getLevelNumber();
    for (int i = 0; i < nLevel; i++) {
      // Compute all numbers from this mappings.

      nRankDofs[i] = computeRankDofs(i);
      nLocalDofs[i] = computeLocalDofs(i);
      nOverallDofs[i] = computeOverallDofs(i);
      rStartDofs[i] = computeStartDofs(i);
      
      // And finally, compute the matrix indices.
      if (needMatIndex)
	computeMatIndex(needMatIndexFromGlobal, i);
    }
  }
  

  void ParallelDofMapping::update(vector<const FiniteElemSpace*>& fe)
  {
    FUNCNAME("ParallelDofMapping::update()");

    for (vector<const FiniteElemSpace*>::iterator it = fe.begin();
	 it != fe.end(); ++it)
      if (find(feSpacesUnique.begin(), feSpacesUnique.end(), *it) == 
	  feSpacesUnique.end())
	ERROR_EXIT("Wrong FE space!\n");

    feSpaces = fe;
    update();
  }

  
  void ParallelDofMapping::computeMatIndex(bool globalIndex, int level)
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");

    TEST_EXIT(level == 0)("Not yet implemented for non-zero level!\n");

    dofToMatIndex.clear();
    
    // The offset is always added to the local matrix index. The offset for the
    // DOFs in the first FE spaces is the smalled global index of a DOF that is
    // owned by the rank.
    int offset = rStartDofs[level];
    

    // === Create the matrix indices for all component FE spaces. ===

    for (unsigned int i = 0; i < feSpaces.size(); i++) {

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
      DofMap& dofMap = data[feSpaces[i]].getMap(level);
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
	if (data[feSpaces[i]].isRankDof(it->first)) {
	  int globalMatIndex = it->second.local + offset;
	  if (globalIndex)
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
	  else
	    dofToMatIndex.add(i, it->first, globalMatIndex);
	}
      }
      
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
      offset += data[feSpaces[i]].nRankDofs[level];
	
      // If there are no non local DOFs, continue with the next FE space.
      if (!hasNonLocalDofs)
	continue;
      
      TEST_EXIT_DBG(sendDofs != NULL && recvDofs != NULL)
	("No communicator given!\n");
      
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

      StdMpi<vector<DegreeOfFreedom> > stdMpi(levelData->getMpiComm(level));
      for (DofComm::Iterator it(*sendDofs, level, feSpaces[i]); 
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
	    if (globalIndex)
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
	    else
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
	
	stdMpi.send(it.getRank(), sendGlobalDofs);
      }
      
      for (DofComm::Iterator it(*recvDofs, level, feSpaces[i]); 
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());
      
      stdMpi.startCommunication();
      
      {
	for (DofComm::Iterator it(*recvDofs, level, feSpaces[i]); 
	     !it.end(); it.nextRank()) {
	  int counter = 0;
	  for (; !it.endDofIter(); it.nextDof()) {
	    if (dofMap.count(it.getDofIndex())) {
	      DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
	      if (globalIndex)
		dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d);
	      else
		dofToMatIndex.add(i, it.getDofIndex(), d);
	    }
	  }
	}
      }
    }
  }
  
}