ParallelDofMapping.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 */

#ifndef AMDIS_FE_SPACE_MAPPING_H
#define AMDIS_FE_SPACE_MAPPING_H

#include <mpi.h>
#include <vector>
#include <map>
#include <set>
#include <petsc.h>
#include <petscis.h>
#include <boost/container/flat_map.hpp>
#include <boost/container/flat_set.hpp>

#include "AMDiS_fwd.h"
#include "parallel/DofComm.h"
#include "parallel/MpiHelper.h"
#include "parallel/ParallelTypes.h"
#include "parallel/StdMpi.h"

namespace AMDiS {

  using namespace std;

  /** \brief
   * Is used to store matrix indices to all DOFs in rank's subdomain. Thus, the
   * class defines a mapping from component number and DOF index to a global
   * matrix index. This class does not calculate the indices by itself!
   */
  class DofToMatIndex
  {
  public:
    DofToMatIndex() {}

    /// Reset the data structure.
    inline void clear()
    {
      data.clear();
    }

    /** Add a new mapping for a given DOF.
     * 
     * \param[in]   component       Component number for which the mapping 
     *                              is defined.
     * \param[in]   dof             DOF index
     * \param[in]   globalMatIndex  Global matrix index.
     */
    inline void add(int component, DegreeOfFreedom dof, int globalMatIndex)
    {
      data[component][dof] = globalMatIndex;
    }

    /// Maps a global DOF index to the global matrix index for a specific 
    /// system component number.
    inline int get(int component, DegreeOfFreedom dof)
    {
      FUNCNAME("DofToMatIndex::get()");

      TEST_EXIT_DBG(data.count(component))
	("No mapping data for component %d available!\n", component);

      TEST_EXIT_DBG(data[component].count(dof))
	("Mapping for DOF %d in component %d does not exists!\n",
	 dof, component);

      return data[component][dof];
    }

    /// Returns for a given matrix index the component and (local or global) DOF
    /// index. As the data structure is not made for this kind of reverse
    /// search, this is very slow and should be only used for debugging.
    void getReverse(int rowIndex, int &component, int &dofIndex);

  private:
    /// The mapping data. For each system component there is a specific map that
    /// maps global DOF indices to global matrix indices.
    map<int, boost::container::flat_map<DegreeOfFreedom, int> > data;
  };


  /**
   * This class defines the parallel mapping of DOFs for one FE space. It is used
   * by the class \ref ParallelDofMapping to specifiy the mapping for a set of 
   * FE spaces.
   */
  class ComponentDofMap
  {
  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.
    ComponentDofMap() 
    {
      ERROR_EXIT("Should not be called!\n");
    }
     
    /// This is the only valid constructur to be used. 
    ComponentDofMap(MeshLevelData* ld);

    /// Clears all data of the mapping.
    void clear();

    /// Maps a DOF index to both, the local and global index of the mapping. The
    /// global index must not be set.
    MultiIndex& operator[](DegreeOfFreedom d)
    {
      TEST_EXIT_DBG(dofMap.count(d))("DOF %d is not in map!\n", d);

      return dofMap[d];
    }

    /** \brief
     * Searchs the map for a given DOF. It does not fail, if the DOF is not 
     * mapped by this mapping. In this case, it returns false. If the DOF is 
     * mapped, the result is stored and the function returns true.
     *
     * \param[in]    dof     DOF index for which a mapping is searched.
     * \param[out]   index   In the case that the DOF is mapped, the result
     *                       is stored here.
     */
    inline bool find(DegreeOfFreedom dof, MultiIndex& index)
    {
      DofMap::iterator it = dofMap.find(dof);
      if (it == dofMap.end())
	return false;

      index = it->second;
      return true;
    }
    
    /// Inserts a new DOF to rank's mapping. The DOF is assumed to be owend by
    /// the rank.
    void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
    {
      FUNCNAME("ComponentDofMap::insertRankDof()");
      
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
      
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      if (dof1 != -1)
	nRankDofs++;
    }
    
    /// Inserts a new DOF to rank's mapping. The DOF exists in rank's subdomain
    /// but is owned by a different rank, thus it is part of an interior boundary.
    void insertNonRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
    {
      FUNCNAME("ComponentDofMap::insertNonRankDof()");
      
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
      
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      nonRankDofs.insert(dof0);
    }
    
    /// Checks if a given DOF is in the DOF mapping.
    bool isSet(DegreeOfFreedom dof)
    {
      return static_cast<bool>(dofMap.count(dof));
    }

    /// Checks if a given DOF is a rank owned DOF of the DOF mapping. The DOF
    /// must be a DOF of the mapping (this is not checked here), otherwise the
    /// result is meaningsless.
    bool isRankDof(DegreeOfFreedom dof)
    {
      return !(static_cast<bool>(nonRankDofs.count(dof)));
    }

    bool isRankGlobalDof(int dof)
    {
      return (dof >= rStartDofs && dof < rStartDofs + nRankDofs);
    }
    
    /// Returns number of DOFs in the mapping.
    unsigned int size()
    {
      return dofMap.size();
    }
    
    /// Returns the raw data of the mapping.
    DofMap& getMap()
    {
      return dofMap;
    }

    DofMap::iterator begin()
    {
      return dofMap.begin();
    }

    DofMap::iterator end()
    {
      return dofMap.end();
    }

    /// Recomputes the mapping.
    void update();

    /// Sets the FE space this mapping corresponds to.
    void setFeSpace(const FiniteElemSpace *fe)
    {
      feSpace = fe;
    }

    /// Informs the mapping whether the mapping will include DOFs that are not
    /// owned by the rank.
    void setNonLocal(bool b)
    {
      isNonLocal = b;
    }

    /// Informs the mapping whether a global index must be computed.
    void setNeedGlobalMapping(bool b)
    {
      needGlobalMapping = b;
    }

    /// Sets the DOF communicator.
    void setDofComm(DofComm &dc)
    {
      dofComm = &dc;
    }

    void setMpiComm(MPI::Intracomm &m, int l)
    {
      mpiComm = m;
      meshLevel = l;
    }

  private:
    /// Computes a global mapping from the local one.
    void computeGlobalMapping();

    /// Computes the global indices of all DOFs in the mapping that are not owned
    /// by the rank.
    void computeNonLocalIndices();

  private:
    MeshLevelData *levelData;

    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;

    MPI::Intracomm mpiComm;

    int meshLevel;

    /// The FE space this mapping belongs to. This is used only the get the
    /// correct DOF communicator in \ref dofComm.
    const FiniteElemSpace *feSpace;

    /// Mapping data from DOF indices to local and global indices.
    DofMap dofMap;

    /// Set of all DOFs that are in mapping but are not owned by the rank.
    boost::container::flat_set<DegreeOfFreedom> nonRankDofs;

    /// If true, a global index mapping will be computed for all DOFs.
    bool needGlobalMapping;

    /// Is true if there are DOFs in at least one subdomain that are not owned
    /// by the rank. If the value is false, each rank contains only DOFs that
    /// are also owned by this rank.
    bool isNonLocal;

  public:
    /// 
    int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;
  };


  class ComponentIterator {
  public:
    virtual ComponentDofMap& operator*() = 0;

    virtual ComponentDofMap* operator->() = 0;

    virtual bool end() = 0;
    
    virtual void next() = 0;

    virtual void reset() = 0;
  };


  class ComponentDataInterface
  {
  public:
    virtual ComponentDofMap& operator[](int compNumber) = 0;

    virtual ComponentDofMap& operator[](const FiniteElemSpace *feSpace) = 0;

    virtual bool isDefinedFor(int compNumber) const = 0;

    virtual ComponentIterator& getIteratorData() = 0;
    
    virtual ComponentIterator& getIteratorComponent() = 0;

    virtual void init(vector<const FiniteElemSpace*> &f0,
		      vector<const FiniteElemSpace*> &f1,
		      bool isNonLocal,
		      MeshLevelData &levelData) = 0;

    vector<const FiniteElemSpace*>& getFeSpaces()
    {
      return feSpaces;
    }

  protected:
    /// The FE spaces for all components.
    vector<const FiniteElemSpace*> feSpaces;

    /// The set of all FE spaces. It uniquly contains all different FE spaces
    /// from \ref feSpaces.
    vector<const FiniteElemSpace*> feSpacesUnique;
  };


  class ComponentDataEqFeSpace : ComponentDataInterface
  {
  public:
    ComponentDataEqFeSpace()
      : iterData(this),
	iterComponent(this)
    {}

    ComponentDofMap& operator[](int compNumber)
    {
      const FiniteElemSpace *feSpace = feSpaces[compNumber];
      return componentData.find(feSpace)->second;
    }

    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
    {
      TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");;
      
      return componentData.find(feSpace)->second;
    }

    bool isDefinedFor(int compNumber) const
    {
      const FiniteElemSpace *feSpace = feSpaces[compNumber];
      return static_cast<bool>(componentData.count(feSpace));
    }

    ComponentIterator& getIteratorData()
    {
      iterData.reset();
      return iterData;
    }
    
    ComponentIterator& getIteratorComponent()
    {
      iterComponent.reset();
      return iterComponent;
    }

    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
	      bool isNonLocal,
	      MeshLevelData &levelData);


  protected:

    void addFeSpace(const FiniteElemSpace* feSpace,
		    MeshLevelData &levelData);

    class IteratorData : public ComponentIterator {
    public:
      IteratorData(ComponentDataEqFeSpace *d)
	: data(d)
      {}

      ComponentDofMap& operator*()
      {
	(*data)[*it];
      }

      ComponentDofMap* operator->()
      {
	&((*data)[*it]);
      }

      bool end()
      {
	return (it != data->feSpacesUnique.end());
      }

      void next()
      {
	++it;
      }

      void reset()
      {
	it = data->feSpacesUnique.begin();
      }

    protected:
      ComponentDataEqFeSpace *data;

      vector<const FiniteElemSpace*>::iterator it;
    };

    class IteratorComponent : public ComponentIterator {
    public:
      IteratorComponent(ComponentDataEqFeSpace *d)
	: data(d)
      {}

      ComponentDofMap& operator*()
      {
	(*data)[*it];
      }

      ComponentDofMap* operator->()
      {
	&((*data)[*it]);
      }

      bool end()
      {
	return (it != data->feSpaces.end());
      }

      void next()
      {
	++it;
      }

      void reset()
      {
	it = data->feSpaces.begin();
      }

    protected:
      ComponentDataEqFeSpace *data;

      vector<const FiniteElemSpace*>::iterator it;
    };
    

    map<const FiniteElemSpace*, ComponentDofMap> componentData;

    IteratorData iterData;

    IteratorComponent iterComponent;

    friend class IteratorData;
    
    friend class IteratorComponent;
  };


  class ComponentDataDiffFeSpace : ComponentDataInterface
  {
  public:
    ComponentDataDiffFeSpace()
      : iter(this)
    {}

    ComponentDofMap& operator[](int compNumber)
    {
      TEST_EXIT_DBG(componentData.count(compNumber))("No data for component %d!\n", compNumber);

      return componentData.find(compNumber)->second;
    }

    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
    {
      ERROR_EXIT("BLUB\n");
    }
    
    ComponentIterator& getIteratorData()
    {
      iter.reset();
      return iter;
    }
    
    ComponentIterator& getIteratorComponent()
    {
      iter.reset();
      return iter;
    }

    bool isDefinedFor(int compNumber) const
    {
      return (static_cast<unsigned int>(compNumber) < componentData.size());
    }

    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
	      bool isNonLocal,
	      MeshLevelData &levelData);

  protected:
    void addComponent(unsigned int component,
		      const FiniteElemSpace* feSpace,
		      MeshLevelData &levelData);

    class Iterator : public ComponentIterator {
    public:
      Iterator(ComponentDataDiffFeSpace *d)
	: data(d),
	  componentCounter(-1)
      {}

      ComponentDofMap& operator*()
      {
	(*data)[componentCounter];
      }

      ComponentDofMap* operator->()
      {
	&((*data)[componentCounter]);
      }


      bool end()
      {
	return (it != data->feSpaces.end());
      }

      void next()
      {
	++it;
	++componentCounter;

	if (it == data->feSpaces.end())
	  componentCounter = -1;
      }

      void reset()
      {
	it = data->feSpaces.begin();
	componentCounter = 0;
      }

    protected:
      ComponentDataDiffFeSpace *data;

      vector<const FiniteElemSpace*>::iterator it;

      int componentCounter;
    };

    map<unsigned int, ComponentDofMap> componentData;

    Iterator iter;

    friend class Iterator;
  };


  /**
   * Implements the mapping from sets of distributed DOF indices to local and
   * global indices. The mapping works for a given set of FE spaces. Furthermore,
   * this class may compute the matrix indices of the set of DOF indices.
   */
  class ParallelDofMapping
  {
  public:
    ParallelDofMapping() 
      : levelData(NULL),
	dofComm(NULL),
	isNonLocal(true),
	needMatIndexFromGlobal(false),
	nRankDofs(1),
	nLocalDofs(1),
	nOverallDofs(1),
	rStartDofs(1)
    {
      nRankDofs = -1;
      nLocalDofs = -1;
      nOverallDofs = -1;
      rStartDofs = -1;
    } 

    /** \brief Initialize the parallel DOF mapping.
     *
     * \param[in]  m                  MPI communicator.
     * \param[in]  fe                 The FE spaces of all components of the 
     *                                PDE to be solved.
     * \param[in]  uniqueFe           Unique list of FE spaces. Thus, two
     *                                arbitrary elements of this list are always
     *                                different.
     * \param[in]  isNonLocal         If true, at least one rank's mapping con-
     *                                taines DOFs that are not owend by the rank.
     */
    void init(MeshLevelData& mld,
	      vector<const FiniteElemSpace*> &fe,
	      vector<const FiniteElemSpace*> &uniqueFe,
	      bool isNonLocal = true);

    /// In the case of having only one FE space, this init function can be used.
    void init(MeshLevelData& mld,
	      const FiniteElemSpace *feSpace,
	      bool isNonLocal = true)
    {
      vector<const FiniteElemSpace*> feSpaces;
      feSpaces.push_back(feSpace);
      init(mld, feSpaces, feSpaces, isNonLocal);
    }

    void setMpiComm(MPI::Intracomm &m, int l);
    
    /// Clear all data.
    void clear();

    /// Set the DOF communicator objects that are required to exchange information
    /// about DOFs that are on interior boundaries.
    void setDofComm(DofComm &dofComm);

    /// Returns the DOF communicator.
    DofComm& getDofComm()
    {
      FUNCNAME("ParallelDofMapping::getDofComm()");

      TEST_EXIT_DBG(dofComm)("No DOF communicator object defined!\n");

      return *dofComm;
    }

    /// Changes the computation of matrix indices based of either local or
    /// global DOF indices, see \ref needMatIndexFromGlobal
    void setComputeMatIndex(bool global)
    {
      needMatIndexFromGlobal = global;
    }

    inline bool isMatIndexFromGlobal()
    {
      return needMatIndexFromGlobal;
    }


    /// Access the DOF mapping for a given component number.
    inline ComponentDofMap& operator[](int compNumber)
    {
      return (*data)[compNumber];
    }

    /// Access the DOF mapping for a given FE space    
    inline ComponentDofMap& operator[](const FiniteElemSpace *feSpace) 
    {
      return (*data)[feSpace];
    }

    inline bool isDefinedFor(int compNumber) const
    {
      return data->isDefinedFor(compNumber);
    }

    /// Returns the number of solution components the mapping is defined on.
    inline int getNumberOfComponents() const
    {      
      return static_cast<int>(data->getFeSpaces().size());
    }

    /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank.
    inline int getRankDofs()
    {
      TEST_EXIT_DBG(nRankDofs >= 0)("Should not happen!\n");

      return nRankDofs;
    }

    /// Returns \ref nLocalDofs, thus the number of DOFs in ranks subdomain.
    inline int getLocalDofs()
    {
      TEST_EXIT_DBG(nLocalDofs >= 0)("Should not happen!\n");

      return nLocalDofs;
    }

    /// Returns \ref nOverallDofs, thus the number of all DOFs in this mapping.
    inline int getOverallDofs()
    {
      TEST_EXIT_DBG(nOverallDofs >= 0)("Should not happen!\n");

      return nOverallDofs;
    }

    /// Returns \ref rStartDofs, thus the smallest global index of a DOF that is
    /// owned by the rank.
    inline int getStartDofs()
    {
      TEST_EXIT_DBG(rStartDofs >= 0)("Should not happen!\n");

      return rStartDofs;
    }

    /// Update the mapping.
    void update();

    /// Update the mapping.
    void update(vector<const FiniteElemSpace*>& feSpaces);

    /// Returns the global matrix index of a given DOF for a given 
    /// component number.
    inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
    {
      return dofToMatIndex.get(ithComponent, d);
    }

    /// Returns the component number and local/global DOF index for a given
    /// matrix row index. Should be used for debugging only!
    inline void getReverseMatIndex(int index, int &component, int &dofIndex)
    {
      dofToMatIndex.getReverse(index, component, dofIndex);
    }

    /// Returns the local matrix index of a given DOF for a given 
    /// component number.
    inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d)
    {
      return dofToMatIndex.get(ithComponent, d) - rStartDofs;
    }

    // inline const FiniteElemSpace* getFeSpace(int i = 0)
    // {      
    //   TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n");
    //   return feSpacesUnique[i];
    // }

    // Writes all data of this object to an output stream.
    void serialize(ostream &out)
    {
      ERROR_EXIT("MUST BE IMPLEMENTED!\n");
    }

    // Reads the object data from an input stream.
    void deserialize(istream &in)
    {
      ERROR_EXIT("MUST BE IMPLEMENTED!\n");
    }

    /// Compute local and global matrix indices.
    void computeMatIndex(bool globalIndex);

    void createIndexSet(IS &is, 
			int firstComponent, 
			int nComponents);

    /// Create a parallel distributed PETSc vector based on this mapping.
    inline void createVec(Vec &vec)
    {
      VecCreateMPI(mpiComm, getRankDofs(), getOverallDofs(), &vec);
    }

    /// Create a parallel distributed PETsc vector based on this mapping but
    /// with a different (larger) global size. This is used in multi-level 
    /// method to embed a local vector into a subdomain spaned by several
    /// ranks.
    inline void createVec(Vec &vec, int nGlobalRows)
    {
      VecCreateMPI(mpiComm, getRankDofs(), nGlobalRows, &vec);
    }

    inline void createLocalVec(Vec &vec)
    {
      VecCreateSeq(PETSC_COMM_SELF, getRankDofs(), &vec);
    }

  protected:
    /// Compute \ref nRankDofs.
    int computeRankDofs();

    /// Compute \ref nLocalDofs.
    int computeLocalDofs();

    /// Compute \ref nOverallDofs.
    int computeOverallDofs();

    /// Compute \ref rStartDofs.
    int computeStartDofs();

  private:
    MPI::Intracomm mpiComm;

    int meshLevel;

    MeshLevelData *levelData;

    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;

    /// Is true if there are DOFs in at least one subdomain that are not owned
    /// by the rank. If the value is false, each rank contains only DOFs that
    /// are also owned by this rank.
    bool isNonLocal;

    /// If matrix indices should be computed, this variable defines if the
    /// mapping from DOF indices to matrix row indices is defined on local
    /// or global DOF indices. If true, the mapping is to specify and to use
    /// on global ones, otherwise on local DOF indices.
    /// In most scenarios the mapping stored on local DOF indices is what we
    /// want to have. Only when periodic boundary conditions are used together
    /// with some global matrix approache, the matrix indices must be stored
    /// also for DOFs that are not part of the local subdomain. Thus, the
    /// mapping will be stored on global DOF indices.
    bool needMatIndexFromGlobal;

    /// Maps from components to DOF mappings.
    ComponentDataInterface *data;

    /// Number of DOFs owned by rank.
    int nRankDofs;

    /// Number of DOFs in rank's subdomain.
    int nLocalDofs;

    /// Number of global DOFs (this value is thus the same on all ranks).
    int nOverallDofs;

    /// Smallest global index of a DOF owned by the rank.
    int rStartDofs;

    /// Mapping from global DOF indices to global matrix indices under 
    /// consideration of possibly multiple components.
    DofToMatIndex dofToMatIndex;
  };
}

#endif