// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file Mesh.h */

/** \defgroup Triangulation Triangulation module
 * @{ <img src="triangulation.png"> @}
 *
 * Example:
 *
 * @{ <img src="hierarchicalMesh.png"> @}
 *
 * \brief
 * Contains all triangulation classes.
 */

#ifndef AMDIS_MESH_H
#define AMDIS_MESH_H

#include <deque>
#include <set>
#include <stdio.h>
#include "AMDiS_fwd.h"
#include "DOFAdmin.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "Element.h"
#include "ElInfo.h"
#include "FixVec.h"
#include "Serializable.h"
#include "BoundaryCondition.h"

namespace AMDiS {

  /** \ingroup Triangulation 
   * \brief
   * A Mesh holds all information about a triangulation. 
   */
  class Mesh : public Serializable
  {
  public:
    /// Creates a mesh with the given name of dimension dim
    Mesh(std::string name, int dim);

    /// Destructor
    virtual ~Mesh();

    /// Reads macro triangulation.
    void initialize();

    /// Assignment operator
    Mesh& operator=(const Mesh&);

    /** \name static methods used while mesh traversal 
     * \{
     */

    /// Used while dof compress
    static int newDOFFct1(ElInfo* e);

    /// Used while dof compress
    static int newDOFFct2(ElInfo* e);
    /** \} */

    /** \name getting methods
     * \{
     */

    /** \brief
     * Returns geometric information about this mesh. With GeoIndex p it is 
     * specifiedm which information is requested.
     */
    inline int getGeo(GeoIndex p) const 
    { 
      return Global::getGeo(p, dim); 
    }

    /// Returns \ref name of the mesh
    inline std::string getName() const 
    { 
      return name; 
    }

    /// Returns \ref dim of the mesh
    inline int getDim() const \
    { 
      return dim; 
    }

    /// Returns \ref nDOFEl of the mesh
    inline const int getNumberOfAllDOFs() const 
    { 
      return nDOFEl; 
    }

    /// Returns \ref nNodeEl of the mesh
    inline const int getNumberOfNodes() const 
    { 
      return nNodeEl; 
    }

    /// Returns \ref nVertices of the mesh
    inline const int getNumberOfVertices() const 
    { 
      return nVertices; 
    }

    /// Returns \ref nEdges of the mesh 
    inline const int getNumberOfEdges() const 
    { 
      return nEdges; 
    }

    /// Returns \ref nFaces of the mesh 
    inline const int getNumberOfFaces() const 
    { 
      return nFaces; 
    }

    /// Returns \ref nLeaves of the mesh 
    inline const int getNumberOfLeaves() const 
    { 
      return nLeaves; 
    }

    /// Returns \ref nElements of the mesh
    inline const int getNumberOfElements() const 
    { 
      return nElements; 
    }

    /// Returns \ref maxEdgeNeigh of the mesh
    inline const int getMaxEdgeNeigh() const 
    { 
      return maxEdgeNeigh; 
    }

    /// Returns \ref parametric of the mesh
    inline Parametric *getParametric() const 
    { 
      return parametric; 
    }

    /// Returns \ref diam of the mesh
    inline const WorldVector<double>& getDiameter() const 
    { 
      return diam; 
    }

    /// Returns nDOF[i] of the mesh
    inline const int getNumberOfDOFs(int i) const 
    { 
      return nDOF[i]; 
    }

    /// Returns \ref elementPrototype of the mesh
    inline Element* getElementPrototype() 
    { 
      return elementPrototype; 
    }

    /// Returns \ref leafDataPrototype of the mesh
    inline ElementData* getElementDataPrototype() 
    { 
      return elementDataPrototype; 
    }

    /// Returns node[i] of the mesh 
    inline int getNode(int i) const 
    { 
      return node[i]; 
    }

    /** \brief
     * Allocates the number of DOFs needed at position and registers the DOFs
     * at the DOFAdmins. The number of needed DOFs is the sum over the needed
     * DOFs of all DOFAdmin objects belonging to this mesh. 
     * The return value is a pointer to the first allocated DOF. 
     */
    DegreeOfFreedom *getDOF(GeoIndex position);

    /// Returns *(\ref admin[i]) of the mesh
    inline const DOFAdmin& getDOFAdmin(int i) const 
    {
      return *(admin[i]);
    }

    /** \brief
     * Creates a DOFAdmin with name lname. nDOF specifies how many DOFs 
     * are needed at the different positions (see \ref DOFAdmin::nrDOF).
     * A pointer to the created DOFAdmin is returned.
     */
    const DOFAdmin* createDOFAdmin(std::string lname, DimVec<int> nDOF);

    /** \brief
     * Returns the size of \ref admin which is the number of the DOFAdmins
     * belonging to this mesh
     */
    const int getNumberOfDOFAdmin() const 
    {
      return admin.size();
    }

    /** \brief
     * Returns the size of \ref macroElements which is the number of
     * of macro elements of this mesh
     */
    const int getNumberOfMacros() const 
    {
      return macroElements.size();
    }

    /// Returns a DOFAdmin which at least manages vertex DOFs
    const DOFAdmin* getVertexAdmin() const;

    /// Allocates a array of DOF pointers. The array holds one pointer for each node.
    DegreeOfFreedom **createDOFPtrs();

    /// Returns \ref preserveCoarseDOFs of the mesh
    inline bool queryCoarseDOFs() const 
    { 
      return preserveCoarseDOFs;
    }

    /// Returns an iterator to the begin of \ref macroElements
    inline std::deque<MacroElement*>::iterator firstMacroElement() 
    {
      return macroElements.begin();
    }

    /// Returns macroElements[i].
    inline MacroElement *getMacroElement(int i) 
    { 
      return macroElements[i]; 
    }

    /// Returns an iterator to the end of \ref macroElements
    inline std::deque<MacroElement*>::iterator endOfMacroElements() 
    {
      return macroElements.end();
    }

    /** \} */

    /** \name setting methods
     * \{
     */

    /// Sets \ref name of the mesh
    inline void setName(std::string aName) 
    { 
      name = aName;
    }

    /// Sets \ref nVertices of the mesh
    inline void setNumberOfVertices(int n) 
    { 
      nVertices = n; 
    }

    /// Sets \ref nFaces of the mesh
    inline void setNumberOfFaces(int n) 
    { 
      nFaces = n; 
    }

    /// Increments \ref nVertices by inc
    inline void incrementNumberOfVertices(int inc) 
    { 
      nVertices += inc; 
    }
 
    /// Sets \ref nEdges of the mesh
    inline void setNumberOfEdges(int n) 
    { 
      nEdges = n; 
    }

    /// Increments \ref nEdges by inc
    inline void incrementNumberOfEdges(int inc) 
    { 
      nEdges += inc; 
    }

    /// Increments \ref nFaces by inc
    inline void incrementNumberOfFaces(int inc) 
    { 
      nFaces += inc; 
    }

    /// Sets \ref nLeaves of the mesh
    inline void setNumberOfLeaves(int n) 
    { 
      nLeaves = n; 
    }

    /// Increments \ref nLeaves by inc
    inline void incrementNumberOfLeaves(int inc) 
    { 
      nLeaves += inc; 
    }

    /// Sets \ref nElements of the mesh
    inline void setNumberOfElements(int n) 
    { 
      nElements = n; 
    }

    /// Increments \ref nElements by inc
    inline void incrementNumberOfElements(int inc) 
    { 
      nElements += inc; 
    }

    /// Sets *\ref diam to w
    void setDiameter(const WorldVector<double>& w);

    /// Sets (*\ref diam)[i] to d
    void setDiameter(int i, double d);

    /// Sets \ref preserveCoarseDOFs = true
    inline void retainCoarseDOFs() 
    {
      preserveCoarseDOFs = true;
    }

    /// Sets \ref preserveCoarseDOFs = b
    inline void setPreserveCoarseDOFs(bool b) 
    {
      preserveCoarseDOFs = b;
    }

    /// Sets \ref preserveCoarseDOFs = false
    inline void noCoarseDOFs() 
    {
      preserveCoarseDOFs = false;
    }

    /// Sets \ref elementPrototype of the mesh
    inline void setElementPrototype(Element* prototype) 
    {
      elementPrototype = prototype;
    }
    
    /// Sets \ref elementDataPrototype of the mesh
    inline void setElementDataPrototype(ElementData* prototype) 
    {
      elementDataPrototype = prototype;
    }

    ///
    inline void setParametric(Parametric *param) 
    {
      parametric = param;
    }

    ///
    inline void setMaxEdgeNeigh(int m) 
    { 
      maxEdgeNeigh = m; 
    }
  
    /** \} */

    /// Creates a new Element by cloning \ref elementPrototype
    Element* createNewElement(Element *parent = NULL);

    /// Creates a new ElInfo dependent of \ref dim of the mesh
    ElInfo* createNewElInfo();

    /// Frees DOFs at the given position pointed by dof 
    void freeDOF(DegreeOfFreedom* dof, GeoIndex position);

    /// Frees memory for the given element el
    void freeElement(Element* el);

    /// Performs DOF compression for all DOFAdmins (see \ref DOFAdmin::compress)
    void dofCompress();

    /// Adds a DOFAdmin to the mesh
    virtual void addDOFAdmin(DOFAdmin *admin);

    /** \brief
     * Traverses the mesh. The argument level specifies the element level if 
     * CALL_EL_LEVEL or CALL_LEAF_EL_LEVEL, or the multigrid level if 
     * CALL_MG_LEVEL is set. Otherwise this variable is ignored. By the argument
     * fillFlag the elements to be traversed and data to be filled into ElInfo is
     * selected, using bitwise or of one CALL_... flag and several FILL_... 
     * flags. The argument elFct is a pointer to a function which is called on 
     * every element selected by the CALL_... part of fillFlag.
     * It is possible to use the recursive mesh traversal recursively, by calling
     * traverse() from elFct.
     */
    int traverse(int level, 
		 const Flag fillFlag,
		 int (*elFct)(ElInfo*));

    /// Recalculates the number of leave elements.
    void updateNumberOfLeaves();

    /// Clears \ref macroElements
    inline void clearMacroElements() 
    { 
      macroElements.clear();
    }
  
    /// Adds a macro element to the mesh
    void addMacroElement(MacroElement* me);

    /* \brief
     * Removes a set of macro elements from the mesh. This works only for the case, 
     * that there are no global or local refinements, i.e., all macro elements have 
     * no children.
     */
    void removeMacroElements(std::vector<MacroElement*>& macros,
			     const FiniteElemSpace* feSpace);

    /// Frees the array of DOF pointers (see \ref createDOFPtrs)
    void freeDOFPtrs(DegreeOfFreedom **ptrs);

    /// Used by \ref findElementAtPoint. 
    bool findElInfoAtPoint(const WorldVector<double>& xy,
			   ElInfo *el_info,
			   DimVec<double>& bary,
			   const MacroElement *start_mel,
			   const WorldVector<double> *xy0,
			   double *sp);

    /** \brief
     * Access to an element at world coordinates xy. Some applications need the 
     * access to elements at a special location in world coordinates. Examples 
     * are characteristic methods for convection problems, or the implementation
     * of a special right hand side like point evaluations or curve integrals.
     * For such purposes, a routine is available which returns an element pointer
     * and corresponding barycentric coordinates.
     *
     * \param xy world coordinates of point
     * \param elp return address for a pointer to the element at xy
     * \param pary returns barycentric coordinates of xy
     * \param start_mel initial guess for the macro element containing xy or NULL
     * \param xy0 start point from a characteristic method, see below, or NULL
     * \param sp return address for relative distance to domain boundary in a 
     *        characteristic method, see below, or NULL
     * \return true is xy is inside the domain , false otherwise
     * 
     * For a characteristic method, where \f$ xy = xy_0 - V\tau \f$, it may be 
     * convenient to know the point on the domain's boundary which lies on the 
     * line segment between the old point xy0 and the new point xy, in case that 
     * xy is outside the domain. Such information is returned when xy0 and a 
     * pointer sp!=NULL are supplied: *sp is set to the value s such that 
     * \f$ xy_0 +s (xy -xy_0) \in \partial Domain \f$, and the element and local 
     * coordinates corresponding to that boundary point will be returned via elp 
     * and bary.
     *
     * The implementation of findElementAtPoint() is based on the transformation 
     * from world to local coordinates, available via the routine worldToCoord(),
     * At the moment, findElementAtPoint() works correctly only for domains with 
     * non-curved boundary. This is due to the fact that the implementation first
     * looks for the macro-element containing xy and then finds its path through 
     * the corresponding element tree based on the macro barycentric coordinates.
     * For non-convex domains, it is possible that in some cases a point inside
     * the domain is considered as external.
     */
    bool findElementAtPoint(const WorldVector<double>& xy,
			    Element **elp, 
			    DimVec<double>& bary,
			    const MacroElement *start_mel,
			    const WorldVector<double> *xy0,
			    double *sp);

    /** \brief
     * Returns for a given dof its world coordinates in this mesh. Because we do
     * not have any direct connection between dofs and coordinates, this function
     * has to search for the element in this mesh, that contains the dof. Than the
     * coordinates can be computed. Therefore, this function is very costly and
     * should be used for debugging purpose only.
     *
     * @param[in]    dof       A pointer to the dof we have to search for.
     * @param[in]    feSpace   The fe space to be used for the search.
     * @param[out]   coords    World vector that stores the coordinates of the dof.
     *
     * The function returns true, if the dof was found, otherwise false.
     */
    bool getDofIndexCoords(const DegreeOfFreedom* dof, 
			   const FiniteElemSpace* feSpace,
			   WorldVector<double>& coords);

    bool getDofIndexCoords(DegreeOfFreedom dof, 
			   const FiniteElemSpace* feSpace,
			   WorldVector<double>& coords);

    /// Returns FILL_ANY_?D
    inline static const Flag& getFillAnyFlag(int dim) 
    {
      switch (dim) {
      case 1:
	return FILL_ANY_1D;
	break;
      case 2:
	return FILL_ANY_2D;
	break;
      case 3:
	return FILL_ANY_3D;
	break;
      default:
	ERROR_EXIT("invalid dim\n");
	return FILL_ANY_1D;
      }
    }

    /// Serialize the mesh to a file.
    void serialize(std::ostream &out);

    /// Deserialize a mesh from a file.
    void deserialize(std::istream &in);

    /// Returns \ref elementIndex and increments it by 1.
    inline int getNextElementIndex() 
    { 
      return elementIndex++; 
    }

    /// Returns \ref initialized.
    inline bool isInitialized() 
    {
      return initialized; 
    }
  
    ///
    inline std::map<BoundaryType, VertexVector*>& getPeriodicAssociations() 
    {
      return periodicAssociations;
    }

    ///
    bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

    ///
    bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

    /// Returns \macroFileInfo
    inline MacroInfo* getMacroFileInfo() 
    { 
      return macroFileInfo;
    }

    ///
    void clearMacroFileInfo();

    /** \brief
     * Traverse this mesh to compute the number of non zero elements the assembled 
     * matrix will have in each row.
     */
    void computeMatrixFillin(const FiniteElemSpace *feSpace,
			     std::vector<int> &nnzInRow, int &overall, int &average);

    ///
    int calcMemoryUsage();

  public:
    ///
    static const Flag FILL_NOTHING;

    ///
    static const Flag FILL_COORDS; 

    ///
    static const Flag FILL_BOUND; 

    ///
    static const Flag FILL_NEIGH; 

    ///
    static const Flag FILL_OPP_COORDS; 

    ///
    static const Flag FILL_ORIENTATION; 

    ///
    static const Flag FILL_ADD_ALL; 
  
    ///
    static const Flag FILL_ANY_1D; 

    ///
    static const Flag FILL_ANY_2D; 

    ///
    static const Flag FILL_ANY_3D; 

    ///
    static const Flag FILL_DET;

    ///
    static const Flag FILL_GRD_LAMBDA;

    //**************************************************************************
    //  flags for Mesh traversal                                                
    //**************************************************************************

    ///
    static const Flag CALL_EVERY_EL_PREORDER;      

    ///
    static const Flag CALL_EVERY_EL_INORDER;     

    ///
    static const Flag CALL_EVERY_EL_POSTORDER;    

    ///
    static const Flag CALL_LEAF_EL;   

    ///
    static const Flag CALL_LEAF_EL_LEVEL;  

    ///
    static const Flag CALL_EL_LEVEL; 

    ///
    static const Flag CALL_MG_LEVEL;

  protected:
    ///
    bool findElementAtPointRecursive(ElInfo *elinfo,
				     const DimVec<double>& lambda,
				     int outside,
				     ElInfo *final_el_info);

  protected:
    /// maximal number of DOFs at one position
    static const int MAX_DOF;

    /// Name of this Mesh
    std::string name;

    /// Dimension of this Mesh. Doesn't have to be equal to dimension of world.
    int dim;

    /// Number of vertices in this Mesh
    int nVertices;

    /// Number of Edges in this Mesh
    int nEdges;

    /// Number of leaf elements in this Mesh
    int nLeaves;

    /// Total number of elements in this Mesh
    int nElements;

    /// Number of faces in this Mesh
    int nFaces;

    /** \brief
     * Maximal number of elements that share one edge; used to allocate memory 
     * to store pointers to the neighbour at the refinement/coarsening edge 
     * (only 3d);
     */
    int maxEdgeNeigh;

    /// Diameter of the mesh in the DIM_OF_WORLD directions
    WorldVector<double> diam;

    /** \brief
     * Is pointer to NULL if mesh contains no parametric elements else pointer 
     * to a Parametric object containing coefficients of the parameterization 
     * and related information
     */
    Parametric *parametric;

    /** \brief
     * If the value is non zero then preserve all DOFs on all levels (can
     * be used for multigrid, e.g.); otherwise all DOFs on the parent that are 
     * not handed over to a child are removed during refinement and added again 
     * on the parent during coarsening.
     */
    bool preserveCoarseDOFs;

    /// Number of all DOFs on a single element
    int nDOFEl;

    /** \brief
     * Number of DOFs at the different positions VERTEX, EDGE, (FACE,) CENTER on
     * an element:
     *
     * - nDOF[VERTEX]: number of DOFs at a vertex (>= 1)
     *
     * - nDOF[EDGE]: number of DOFs at an edge; if no DOFs are associated to
     *   edges, then this value is 0
     *
     * - nDOF[FACE]: number of DOFs at a face; if no DOFs are associated to
     *   faces, then this value is 0 (only 3d)
     *
     * - nDOF[CENTER]: number of DOFs at the barycenter; if no DOFs are 
     *   associated to the barycenter, then this value is 0
     */
    DimVec<int> nDOF;

    /** \brief
     * Number of nodes on a single element where DOFs are located; needed for 
     * the (de-) allocation of the dof-vector on the element (\ref Element::dof);
     */
    int nNodeEl;

    /** \brief
     * Gives the index of the first node at vertex, edge, face (only 3d), and 
     * barycenter:
     *
     * - node[VERTEX]: has always value 0; dof[0],...,dof[N_VERTICES-1] are 
     *   always DOFs at the vertices;
     *
     * - node[EDGE]: dof[node[EDGE]],..., dof[node[EDGE]+N_EDGES-1] are the DOFs
     *   at the N_EDGES edges, if DOFs are located at edges;
     *
     * - node[FACE]: dof[node[FACE]],..., dof[node[FACE]+N_FACES-1] are the DOFs
     *   at the N_FACES faces, if DOFs are located at faces (only 3d);
     *
     * - node[CENTER]: dof[node[CENTER]] are the DOFs at the barycenter, if DOFs
     *   are located at the barycenter;
     */
    DimVec<int> node;

    /// List of all DOFAdmins
    std::vector<DOFAdmin*> admin;

    /// List of all MacroElements of this Mesh
    std::deque<MacroElement*> macroElements;

    /// Needed during DOF compression (\ref DOFAdmin::compress).
    std::vector<DegreeOfFreedom> newDOF;

    /// Needed during DOF compression (\ref DOFAdmin::compress).
    static DOFAdmin *compressAdmin;

    /** \brief
     * Used for recursive mesh traversal. Static pointer to the mesh
     * that should be traversed. This allows access to the mesh even
     * from the static traverse routines
     */
    static Mesh* traversePtr;

    /// Used by check functions
    static std::vector<DegreeOfFreedom> dof_used;

    /** \brief
     * This map is used for serialization and deserialization of mesh elements. 
     * During the serialization process, all elements are visited and their dof indices
     * are written to the file. If a dof index at a position, i.e. vertex, line or face,
     * was written to file, the combination of dof index and position is inserted to
     * this map. That ensures that the same dof at the same position, but being part of
     * another element, is not written twice to the file. 
     * When a state should be deserialized, the information can be used to construct
     * exactly the same dof structure.
     */
    static std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> serializedDOFs;

    /** \brief
     * Used while mesh refinement. To create new elements 
     * elementPrototype->clone() is called, which returns a Element of the
     * same type as elementPrototype. So e.g. Elements of the different
     * dimensions can be created in a uniform way. 
     */
    Element* elementPrototype;

    /// Prototype for leaf data. Used for creation of new leaf data while refinement.
    ElementData* elementDataPrototype;

    /// Used for enumeration of all mesh elements
    int elementIndex;

    /// True if the mesh is already initialized, false otherwise.
    bool initialized;

    /// Map of managed periodic vertex associations.
    std::map<BoundaryType, VertexVector*> periodicAssociations;

    /** \brief
     * If the mesh has been created by reading a macro file, here 
     * the information are stored about the content of the file.
     */
    MacroInfo *macroFileInfo;

  protected:
    /// for findElement-Fcts
    DimVec<double> final_lambda;

    /** \brief
     * Temporary variables that are used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
     */
    const WorldVector<double> *g_xy0, *g_xy;

    /** \brief
     * Temporary variable that is used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
     */    
    double *g_sp;
   
    friend class MacroInfo;
    friend class MacroReader;
    friend class MacroWriter;
    friend class MacroElement;
    friend class Element;
    friend void Element::newDOFFct1(const DOFAdmin*);
    friend void Element::newDOFFct2(const DOFAdmin*);
  };

}

#endif  // AMDIS_MESH_H