FeSpaceMapping.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// 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.



/** \file FeSpaceMapping.h */

#include <vector>
#include <map>
#include <set>
#include "parallel/DofComm.h"
#include "parallel/MpiHelper.h"
#include "parallel/ParallelTypes.h"

#ifndef AMDIS_FE_SPACE_MAPPING_H
#define AMDIS_FE_SPACE_MAPPING_H

namespace AMDiS {
  using namespace std;

  struct MultiIndex
  {
    int local, global;
  };

  class GlobalDofMap
  {
  public:
    /// This constructor exists only to create std::map of this class and make
    /// use of the operator [] for read access. Should never be called.
    GlobalDofMap() 
    {
      ERROR_EXIT("Should not be called!\n");
    }
      
    GlobalDofMap(MPI::Intracomm* m)
      : mpiComm(m),
	feSpace(NULL),
	sendDofs(NULL),
	recvDofs(NULL),
	needGlobalMapping(false),
	nRankDofs(0),
	nLocalDofs(0),
	nOverallDofs(0),
	rStartDofs(0),
	overlap(false)
    {}
    
    void clear();
    
    MultiIndex& operator[](DegreeOfFreedom d)
    {
      TEST_EXIT_DBG(dofMap.count(d))("Should not happen!\n");

      return dofMap[d];
    }
    
    void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
    {
      FUNCNAME("GlobalDofMap::insertRankDof()");
      
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
      
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      if (dof1 != -1)
	nRankDofs++;
    }
    
    void insert(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
    {
      FUNCNAME("GlobalDofMap::insert()");
      
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
      
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      nonRankDofs.insert(dof0);
    }
    
    bool isSet(DegreeOfFreedom dof)
    {
      return static_cast<bool>(dofMap.count(dof));
    }
    
    unsigned int size()
    {
      return dofMap.size();
    }
    
    map<DegreeOfFreedom, MultiIndex>& getMap()
    {
      return dofMap;
    }
    
    void update();

    void addOffset(int offset);

    void computeNonLocalIndices();

    void computeMatIndex(int offset);

    void print();

    void setFeSpace(const FiniteElemSpace *fe)
    {
      feSpace = fe;
    }

    void setOverlap(bool b)
    {
      overlap = b;
    }

    void setDofComm(DofComm &pSend, DofComm &pRecv)
    {
      sendDofs = &pSend;
      recvDofs = &pRecv;
    }

    void setNeedGlobalMapping(bool b)
    {
      needGlobalMapping = b;
    }

  private:
    MPI::Intracomm* mpiComm;

    const FiniteElemSpace *feSpace;

    /// 
    map<DegreeOfFreedom, MultiIndex> dofMap;

    std::set<DegreeOfFreedom> nonRankDofs;

    DofComm *sendDofs;
    
    DofComm *recvDofs;

    bool needGlobalMapping;
  public:
    /// 
    int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;

    bool overlap;
  };
  

  template<typename T>
  class FeSpaceData
  {
  public:
    FeSpaceData() 
      : mpiComm(NULL),
	feSpaces(0),
	nRankDofs(-1),
	nOverallDofs(-1),
	rStartDofs(-1)
    {} 

    T& operator[](const FiniteElemSpace* feSpace)
    {
      FUNCNAME("FeSpaceData::operator[]()");

      TEST_EXIT_DBG(data.count(feSpace))("Should not happen!\n");

      return data.find(feSpace)->second;
    }

    void addFeSpace(const FiniteElemSpace* feSpace)
    {
      FUNCNAME("FeSpaceData::addFeSpace()");
      
      if (data.count(feSpace))
	data.find(feSpace)->second.clear();
      else
	data.insert(make_pair(feSpace, T(mpiComm)));

      data.find(feSpace)->second.setFeSpace(feSpace);
    }    

    int getRankDofs(vector<const FiniteElemSpace*> &fe)
    {
      FUNCNAME("FeSpaceData::getRankDofs()");

      int result = 0;
      for (unsigned int i = 0; i < fe.size(); i++) {
	TEST_EXIT_DBG(data.count(fe[i]))("Cannot find FE space: %p\n", fe[i]);
	result += data[fe[i]].nRankDofs;
      }

      return result;
    }

    inline int getRankDofs()
    {
      TEST_EXIT_DBG(nRankDofs >= 0)("Should not happen!\n");

      return nRankDofs;
    }

    int getLocalDofs(vector<const FiniteElemSpace*> &fe)
    {
      FUNCNAME("FeSpaceData::getLocalDofs()");

      int result = 0;
      for (unsigned int i = 0; i < fe.size(); i++) {
	TEST_EXIT_DBG(data.count(fe[i]))("Cannot find FE space: %p\n", fe[i]);
	result += data[fe[i]].nLocalDofs;
      }

      return result;
    }

    inline int getLocalDofs()
    {
      TEST_EXIT_DBG(nLocalDofs >= 0)("Should not happen!\n");

      return nLocalDofs;
    }

    int getOverallDofs(vector<const FiniteElemSpace*> &feSpaces)
    {
      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;
      }

      return result;
    }

    inline int getOverallDofs()
    {
      TEST_EXIT_DBG(nOverallDofs >= 0)("Should not happen!\n");

      return nOverallDofs;
    }

    int getStartDofs(vector<const FiniteElemSpace*> &feSpaces)
    {
      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;
      }

      return result;
    }

    int getStartDofs()
    {
      TEST_EXIT_DBG(rStartDofs >= 0)("Should not happen!\n");

      return rStartDofs;
    }

    void init(MPI::Intracomm *m,
	      vector<const FiniteElemSpace*> &fe,
	      bool needGlobalMapping)
    {
      mpiComm = m;
      feSpaces = fe;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	feSpacesUnique.insert(feSpaces[i]);
	
	addFeSpace(feSpaces[i]);
	data[feSpaces[i]].setNeedGlobalMapping(needGlobalMapping);
      }
    }

    void update()
    {
      for (std::set<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	   it != feSpacesUnique.end(); ++it)
	data[*it].update();

      nRankDofs = getRankDofs(feSpaces);
      nLocalDofs = getLocalDofs(feSpaces);
      nOverallDofs = getOverallDofs(feSpaces);
      rStartDofs = getStartDofs(feSpaces);
    }

    inline int mapLocal(int index, int ithFeSpace)
    {
      int result = 0;
      for (int i = 0; i < ithFeSpace; i++)
	result += data[feSpaces[i]].nRankDofs;
      result += data[feSpaces[ithFeSpace]][index].local;

      return result;
    }

    inline int mapLocal(int index, const FiniteElemSpace *feSpace)
    {
      for (unsigned int i = 0; i < feSpaces.size(); i++)
	if (feSpaces[i] == feSpace)
	  return mapLocal(index, feSpace, i);

      return -1;
    }

    inline int mapGlobal(int index, int ithFeSpace)
    {
      int result = rStartDofs;
      for (int i = 0; i < ithFeSpace; i++)
	result += data[feSpaces[i]].nRankDofs;
      result += data[feSpaces[ithFeSpace]][index].local;

      return result;
    }

    inline int mapGlobal(int index, const FiniteElemSpace *feSpace)
    {
      for (unsigned int i = 0; i < feSpaces.size(); i++)
	if (feSpaces[i] == feSpace)
	  return mapGlobal(index, feSpace, i);

      return -1;
    }

  private:
    MPI::Intracomm* mpiComm;
    
    map<const FiniteElemSpace*, T> data;

    vector<const FiniteElemSpace*> feSpaces;

    std::set<const FiniteElemSpace*> feSpacesUnique;

    int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;
  };

}

#endif