ElementObjectDatabase.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 ElementObjectDatabase.h */

#ifndef AMDIS_ELEMENT_OBJECT_DATABASE_H
#define AMDIS_ELEMENT_OBJECT_DATABASE_H

#include <map>
#include <vector>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>

#include "AMDiS_fwd.h"
#include "Containers.h"
#include "Global.h"
#include "Boundary.h"
#include "Serializer.h"
#include "FiniteElemSpace.h"

namespace AMDiS {

  using namespace std;

  /// Just to templatize the typedef.
  template<typename T>
  struct PerBoundMap {
    typedef map<pair<T, T>, BoundaryType> type;
    typedef typename type::iterator iterator;
  };


  /// Defines one element object. This may be either a vertex, edge or face.
  struct ElementObjectData {
    ElementObjectData(int a = -1, int b = 0)
      : elIndex(a),
	ithObject(b)
    {}

    /// Index of the element this object is part of.
    int elIndex;
    
    /// Index of the object within the element.
    int ithObject;
    
    /// Write this element object to disk.
    void serialize(ostream &out) const
    {
      SerUtil::serialize(out, elIndex);
      SerUtil::serialize(out, ithObject);
    }

    /// Read this element object from disk.
    void deserialize(istream &in)
    {
      SerUtil::deserialize(in, elIndex);
      SerUtil::deserialize(in, ithObject);
    }

    /// Compare this element object with another one.
    bool operator==(ElementObjectData& cmp) const
    {
      return (elIndex == cmp.elIndex && ithObject == cmp.ithObject);
    }

    /// Define a strict order on element objects.
    bool operator<(const ElementObjectData& rhs) const
    {
      return (elIndex < rhs.elIndex || 
	      (elIndex == rhs.elIndex && ithObject < rhs.ithObject));
    }
  };



  /** \brief
   * This class is a database of element objects. An element object is either a
   * vertex, edge or the face of a specific element. This database is used to
   * store all objects of all elements of a mesh. The information is stored in a
   * way that makes it possible to identify all elements, which have a given
   * vertex, edge or face in common. If is is known which element is owned by 
   * which rank in parallel computations, it is thus possible to get all interior
   * boundaries on object level. This is required, because two elements may share
   * a common vertex without beging neighbours in the definition of AMDiS.
   */
  class ElementObjectDatabase {
  public:
    ElementObjectDatabase()
      : feSpace(NULL),
	mesh(NULL),
	iterGeoPos(CENTER)
    {}

    void setFeSpace(const FiniteElemSpace *fe)
    {
      feSpace = fe;
      mesh = feSpace->getMesh();
    }
  
    Mesh* getMesh()
    {
      return mesh;
    }

    void create();

    void createMacroElementInfo(vector<MacroElement*> &mel);

    /** \brief
     * Create for a filled object database the membership information for all
     * element objects. An object is owned by a rank, if the rank has the
     * heighest rank number of all ranks where the object is part of.
     *
     * \param[in]  macroElementRankMap   Maps to each macro element of the mesh
     *                                   the rank that owns this macro element.
     */
    void createRankData(map<int, int>& macroElementRankMap,
			MeshLevelData& levelData);


    /** \brief
     * Iterates over all elements for one geometrical index, i.e., over all
     * vertices, edges or faces in the mesh. The function returns true, if the
     * result is valid. Otherwise the iterator is at the end position.
     *
     * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the
     *                   elements that should be traversed.
     */
    bool iterate(GeoIndex pos)
    {
      // CENTER marks the variable "iterGeoPos" to be in an undefined state. I.e.,
      // there is no iteration that is actually running.

      if (iterGeoPos == CENTER) {
	iterGeoPos = pos;
	switch (iterGeoPos) {
	case VERTEX:
	  vertexIter = vertexInRank.begin();
	  break;
	case EDGE:
	  edgeIter = edgeInRank.begin();
	  break;
	case FACE:
	  faceIter = faceInRank.begin();
	  break;
	default:
	  ERROR_EXIT("Not GeoIndex %d!\n", iterGeoPos);
	}
      } else {
	switch (iterGeoPos) {
	case VERTEX:
	  ++vertexIter;
	  break;
	case EDGE:
	  ++edgeIter;
	  break;
	case FACE:
	  ++faceIter;
	  break;
	default:
	  ERROR_EXIT("Not GeoIndex %d!\n", iterGeoPos);
	}
      }

      switch (iterGeoPos) {
      case VERTEX:
	if (vertexIter == vertexInRank.end()) {
	  iterGeoPos = CENTER;
	  return false;
	}
	break;
      case EDGE:
	if (edgeIter == edgeInRank.end()) {
	  iterGeoPos = CENTER;
	  return false;
	}
	break;
      case FACE:
	if (faceIter == faceInRank.end()) {
	  iterGeoPos = CENTER;
	  return false;
	}
	break;
      default:
	ERROR_EXIT("Should not happen!\n");	
      }

      return true;
    }


    /// Returns the data of the current iterator position.
    map<int, ElementObjectData>& getIterateData()
    {
      switch (iterGeoPos) {
      case VERTEX:
	return vertexIter->second;
	break;
      case EDGE:
	return edgeIter->second;
	break;
      case FACE:
	return faceIter->second;
	break;
      default:
	ERROR_EXIT("Should not happen!\n");

	// Will never be reached, just to avoid compiler warnings.
	return faceIter->second;
      }
    }


    /// Returns the rank owner of the current iterator position.
    int getIterateOwner()
    {
      switch (iterGeoPos) {
      case VERTEX:
	return vertexOwner[vertexIter->first];
	break;
      case EDGE:
	return edgeOwner[edgeIter->first];
	break;
      case FACE:
	return faceOwner[faceIter->first];
	break;
      default:
	ERROR_EXIT("Should not happen!\n");

	// Will never be reached, just to avoid compiler warnings.
	return -1;
      }
    }

    /// Returns the rank owner of the current iterator position.
    int getIterateMaxLevel()
    {
      switch (iterGeoPos) {
      case VERTEX:
	return vertexMaxLevel[vertexIter->first];
	break;
      case EDGE:
	return edgeMaxLevel[edgeIter->first];
	break;
      case FACE:
	return faceMaxLevel[faceIter->first];
	break;
      default:
	ERROR_EXIT("Should not happen!\n");

	// Will never be reached, just to avoid compiler warnings.
	return -1;
      }
    }

    /// Returns the rank owner of a vertex DOF.
    int getOwner(DegreeOfFreedom vertex)
    {
      return vertexOwner[vertex];
    }

    /// Returns the rank owner of an edge.
    int getOwner(DofEdge edge)
    {
      return edgeOwner[edge];
    }

    /// Returns the rank owner of an face.
    int getOwner(DofFace face, int level)
    {
      return faceOwner[face];
    }


    /// Checks if a given vertex DOF is in a given rank.
    int isInRank(DegreeOfFreedom vertex, int rank)
    {
      return (vertexInRank[vertex].count(rank));
    }

    /// Checks if a given edge is in a given rank.
    int isInRank(DofEdge edge, int rank)
    {
      return (edgeInRank[edge].count(rank));
    }

    /// Checks if a given face is in a given rank.
    int isInRank(DofFace face, int rank)
    {
      return (faceInRank[face].count(rank));
    }


    /// Returns a vector with all macro elements that have a given vertex DOF 
    /// in common.
    vector<ElementObjectData>& getElements(DegreeOfFreedom vertex)
    {
      return vertexElements[vertex];
    }

    /// Returns a vector with all macro elements that have a given edge in common.
    vector<ElementObjectData>& getElements(DofEdge edge)
    {
      return edgeElements[edge];
    }

    /// Returns a vector with all macro elements that have a given face in common.
    vector<ElementObjectData>& getElements(DofFace face)
    {
      return faceElements[face];
    }


    /// Returns a vector with all macro elements that have a given vertex DOF 
    /// in common.
    vector<ElementObjectData>& getElementsVertex(int elIndex, int ithVertex)
    {
      ElementObjectData elObj(elIndex, ithVertex);
      DegreeOfFreedom vertex = vertexLocalMap[elObj];
      return vertexElements[vertex];
    }
    
    /// Returns a vector with all macro elements that have a given edge in common.
    vector<ElementObjectData>& getElementsEdge(int elIndex, int ithEdge)
    {
      ElementObjectData elObj(elIndex, ithEdge);
      DofEdge edge = edgeLocalMap[elObj];
      return edgeElements[edge];
    }

    /// Returns a vector with all macro elements that have a given face in common.
    vector<ElementObjectData>& getElementsFace(int elIndex, int ithFace)
    {
      ElementObjectData elObj(elIndex, ithFace);
      DofFace face = faceLocalMap[elObj];
      return faceElements[face];
    }


    
    /// Returns a map that maps to each rank all macro elements in this rank that
    /// have a given vertex DOF in common.
    map<int, ElementObjectData>& getElementsInRank(DegreeOfFreedom vertex)
    {
      return vertexInRank[vertex];
    }

    /// Returns a map that maps to each rank all macro elements in this rank that
    /// have a given edge in common.
    map<int, ElementObjectData>& getElementsInRank(DofEdge edge)
    {
      return edgeInRank[edge];
    }

    /// Returns a map that maps to each rank all macro elements in this rank that
    /// have a given face in common.
    map<int, ElementObjectData>& getElementsInRank(DofFace face)
    {
      return faceInRank[face];
    }

    /// Returns to an element object data the appropriate vertex DOF.
    DegreeOfFreedom getVertexLocalMap(ElementObjectData &data)
    {
      TEST_EXIT_DBG(vertexLocalMap.count(data))("Should not happen!\n");

      return vertexLocalMap[data];
    }

    /// Returns to an element object data the appropriate edge.
    DofEdge getEdgeLocalMap(ElementObjectData &data)
    {
      TEST_EXIT_DBG(edgeLocalMap.count(data))("Should not happen!\n");

      return edgeLocalMap[data];
    }

    /// Returns to an element object data the appropriate face.
    DofFace getFaceLocalMap(ElementObjectData &data)
    {
      TEST_EXIT_DBG(faceLocalMap.count(data))("Should not happen!\n");

      return faceLocalMap[data];
    }

    PerBoundMap<DegreeOfFreedom>::type& getPeriodicVertices()
    {
      return periodicVertices;
    }

    PerBoundMap<DofEdge>::type& getPeriodicEdges()
    {
      return periodicEdges;
    }

    PerBoundMap<DofFace>::type& getPeriodicFaces()
    {
      return periodicFaces;
    }

    inline bool getEdgeReverseMode(ElementObjectData &obj0, 
				   ElementObjectData &obj1)
    {
      if (mesh->getDim() == 2)
	return true;

      TEST_EXIT_DBG(edgeReverseMode.count(make_pair(obj0, obj1)))
	("Should not happen!\n");

      return edgeReverseMode[make_pair(obj0, obj1)];
    }

    inline bool getFaceReverseMode(ElementObjectData &obj0, 
				   ElementObjectData &obj1)
    {
      TEST_EXIT_DBG(faceReverseMode.count(make_pair(obj0, obj1)))
	("Should not happen!\n");

      return faceReverseMode[make_pair(obj0, obj1)];
    }

    /// Returns true if there is periodic data.
    bool hasPeriodicData()
    {
      return (periodicVertices.size() != 0);
    }

    /// Returns true if the given boundary type is larger or equal to the smallest
    /// periodic boundary ID in mesh. See \ref smallestPeriodicBcType for more
    /// information.
    bool isValidPeriodicType(BoundaryType t) const
    {
      return (t >= smallestPeriodicBcType);
    }

    inline Element* getElementPtr(int index)
    {
      return macroElIndexMap[index];
    }

    inline int getElementType(int index)
    {
      return macroElIndexTypeMap[index];
    }

    /// Write the element database to disk.
    void serialize(ostream &out);
    
    /// Read the element database from disk.
    void deserialize(istream &in);

  protected:
    /** \brief
     * Adds an element to the object database. If the element is part of a
     * periodic boundary, all information about subobjects of the element on
     * this boundary are collected.
     *
     * \param[in]  elInfo    ElInfo object of the element. 
     */
    void addElement(ElInfo *elInfo);

    /// Adds the i-th DOF vertex of an element to the object database.
    void addVertex(Element *el, int ith)
    {
      DegreeOfFreedom vertex = el->getDof(ith, 0);
      int elIndex = el->getIndex();
      ElementObjectData elObj(elIndex, ith);

      vertexElements[vertex].push_back(elObj);
      vertexLocalMap[elObj] = vertex;
    }

    /// Adds the i-th edge of an element to the object database.
    void addEdge(Element *el, int ith)
    {
      DofEdge edge = el->getEdge(ith);
      int elIndex = el->getIndex();
      ElementObjectData elObj(elIndex, ith);

      edgeElements[edge].push_back(elObj);
      edgeLocalMap[elObj] = edge;
    }

    /// Adds the i-th face of an element to the object database.
    void addFace(Element *el, int ith)
    {
      DofFace face = el->getFace(ith);
      int elIndex = el->getIndex();
      ElementObjectData elObj(elIndex, ith);

      faceElements[face].push_back(elObj);
      faceLocalMap[elObj] = face;
    }

    /** \brief
     * Creates final data of the periodic boundaries. Must be called after all
     * elements of the mesh are added to the object database. Then this functions
     * search for indirectly connected vertices in periodic boundaries. This is
     * only the case, if there are more than one boundary conditions. Then, e.g., 
     * in 2D, all edges of a square are iterectly connected. In 3D, if the macro 
     * mesh is a box, all eight vertex nodes and always four of the 12 edges are 
     * indirectly connected.
     */
    void createPeriodicData();


    /// Creates on all boundaries the reverse mode flag.
    void createReverseModeData();

    BoundaryType getNewBoundaryType();

    BoundaryType provideConnectedPeriodicBoundary(BoundaryType b0, 
						  BoundaryType b1);

    /// Some auxiliary function to write the element object database to disk.
    void serialize(ostream &out, vector<ElementObjectData>& elVec);

    /// Some auxiliary function to read the element object database from disk.
    void deserialize(istream &in, vector<ElementObjectData>& elVec);

    /// Some auxiliary function to write the element object database to disk.
    void serialize(ostream &out, map<int, ElementObjectData>& data);

    /// Some auxiliary function to read the element object database from disk.
    void deserialize(istream &in, map<int, ElementObjectData>& data);

  private:
    const FiniteElemSpace* feSpace;

    /// The mesh that is used to store all its element information in 
    /// the database.
    Mesh *mesh;
    
    /// Maps to each vertex DOF all element objects that represent this vertex.
    map<DegreeOfFreedom, vector<ElementObjectData> > vertexElements;

    /// Maps to each edge all element objects that represent this edge.
    map<DofEdge, vector<ElementObjectData> > edgeElements;

    /// Maps to each face all element objects that represent this edge.
    map<DofFace, vector<ElementObjectData> > faceElements;

    
    /// Maps to an element object the corresponding vertex DOF.
    map<ElementObjectData, DegreeOfFreedom> vertexLocalMap;

    /// Maps to an element object the corresponding edge.
    map<ElementObjectData, DofEdge> edgeLocalMap;

    /// Maps to an element object the corresponding face.
    map<ElementObjectData, DofFace> faceLocalMap;

    /// Defines for all vertex DOFs the rank that ownes this vertex DOF.
    map<DegreeOfFreedom, int> vertexOwner;

    /// Defines for all edges the rank that ownes this edge.
    map<DofEdge, int> edgeOwner;

    /// Defines for all faces the rank that ownes this face.
    map<DofFace, int> faceOwner;

    
    map<DegreeOfFreedom, int> vertexMaxLevel;

    map<DofEdge, int> edgeMaxLevel;

    map<DofFace, int> faceMaxLevel;


    /// Defines to each vertex DOF a map that maps to each rank number the element
    /// objects that have this vertex DOF in common.
    map<DegreeOfFreedom, map<int, ElementObjectData> > vertexInRank;

    /// Defines to each edge a map that maps to each rank number the element 
    /// objects that have this edge in common.
    map<DofEdge, map<int, ElementObjectData> > edgeInRank;

    /// Defines to each face a map that maps to each rank number the element 
    /// objects that have this face in common.
    map<DofFace, map<int, ElementObjectData> > faceInRank;


    /// Vertex iterator to iterate over \ref vertexInRank
    map<DegreeOfFreedom, map<int, ElementObjectData> >::iterator vertexIter;

    /// Edge iterator to iterate over \ref edgeInRank
    map<DofEdge, map<int, ElementObjectData> >::iterator edgeIter;

    /// Face iterator to iterate over \ref faceInRank
    map<DofFace, map<int, ElementObjectData> >::iterator faceIter;


    /// Defines the geometrical iteration index of the iterators. I.e., the value
    /// is either VERTEX, EDGE or FACE and the corresponding element objects are
    /// traversed. The value CENTER is used to define a not defined states of the
    /// iterators, i.e., if no iteration is running.
    GeoIndex iterGeoPos;

    map<pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;

    /// The following three data structures store periodic DOFs, edges and faces.
    PerBoundMap<DegreeOfFreedom>::type periodicVertices;
    PerBoundMap<DofEdge>::type periodicEdges;
    PerBoundMap<DofFace>::type periodicFaces;

    /// Defines the smallest boudary ID for periodic boundary conditions. This is
    /// required to distinguish between "real" periodic boundaries and periodic
    /// boundary IDs that are set by the parallel algorithm for indirectly 
    /// connected boundaries.
    BoundaryType smallestPeriodicBcType;

    /// Stores to each vertex all its periodic associations.
    map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;

    /// Stores to each edge all its periodic associations.
    map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;

    map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode;

    map<pair<ElementObjectData, ElementObjectData>, bool> faceReverseMode;

    /// Maps to each macro element index a pointer to the corresponding element.
    map<int, Element*> macroElIndexMap;
    
    /// Maps to each macro element index the type of this element.
    map<int, int> macroElIndexTypeMap;
  };

}

#endif