#ifndef _Mesh_h
  #define _Mesh_h

  #include <stdio.h>
  #include "Types.h"
  #include "Triangle.h"


  typedef struct _tMesh_Point{
    double2 Point;
    struct _tMesh_VolumeElement *Ball;
    struct _tMesh_SurfaceElement *SurfaceBall;
    BOOL Hanging;  //Bisection: point is a hanging node
    /*    float RiemannMetric[6];  //valid if Anisotropic_Init was called */
    /*    float w[4];  //rho, v, p  */
    struct _tMesh_Point *PeriodicCopy;  //0 if not periodic
    struct _tMesh_Point *Prev,*Next;
    struct _tAdj_Mesh_Point *AdjPoint;
    }tMesh_Point;

  typedef struct _tMesh_VolumeElement{
    tMesh_Point *P[3];  //corner points
    statevec w;
    double q[3]; /* for visualization */
    /*    double rho;
    double2 v;
    double p; */
    struct _tMesh_VolumeElement *Neighbour[3];  //neighbours
    struct _tMesh_SurfaceElement *Surface[3];
    struct _tMesh_VolumeElement *BallPrev[3],*BallNext[3];
    int splitcounter;  //how many times to split
    struct _tMesh_VolumeElement *Parent;  //for bisection
    char r,m[3],f;
    BOOL Last;  //was moved by Mesh_VolumeElementLast?
    struct _tMesh_VolumeElement *Prev,*Next;
    }tMesh_VolumeElement;

  typedef struct _tMesh_SurfaceElement{
    tMesh_Point *P[2];  //corner points
    uint surfnr;  //surface number
    struct _tMesh_SurfaceElement *SurfaceBallPrev[2],*SurfaceBallNext[2];
    struct _tMesh_SurfaceElement *Prev,*Next;
    uint curveelement;  /* the boundary element lies on curved part of bound */
    double2 xb;  /* only if  curveelemnt =1 */
    }tMesh_SurfaceElement;

/*  typedef struct _tMesh_Red{
    tMesh_VolumeElement *Element;  //red element
    struct _tMesh_Red *Prev,*Next;
    }tMesh_Red; */

  typedef struct{
    tMesh_Point *Point;
    tMesh_VolumeElement *VolumeElement,*VolumeElementLast;
    tMesh_SurfaceElement *SurfaceElement;
    double kappa;  //Poisson's constant
/*    tMesh_Red *Red; */
    uint npoin,nelem,nbound,nred,ncurvedelements;
    int nper, iper[4]; /* parametres for periodicity */
    double rper[4];    /* maximum two pairs of boundary parts */
    BOOL NeighboursFound;  //TRUE if Neighbour and Surface items are correct
    tSurface_Surface *Surface;
    }tMesh_Mesh;


  tMesh_Mesh *Mesh_Create(void);
  void Mesh_Destroy(tMesh_Mesh *);

  tMesh_Point *Mesh_PointAdd(tMesh_Mesh *,double,double);
  void Mesh_PointDelete(tMesh_Mesh *,tMesh_Point *);

  tMesh_VolumeElement *Mesh_VolumeElementAdd(tMesh_Mesh *,tMesh_Point *,
    tMesh_Point *,tMesh_Point *,double,double,double,double,double,double,
  double,double,double,double,double,double);
  void Mesh_VolumeElementUnqueue(tMesh_Mesh *,tMesh_VolumeElement *);
  void Mesh_VolumeElementDelete(tMesh_Mesh *,tMesh_VolumeElement *);
  void Mesh_VolumeElementLast(tMesh_Mesh *,tMesh_VolumeElement *);

  tMesh_SurfaceElement *Mesh_SurfaceElementAdd(tMesh_Mesh *,tMesh_Point *,
    tMesh_Point *,uint);
  void Mesh_SurfaceElementDelete(tMesh_Mesh *,tMesh_SurfaceElement *);

/*  tMesh_Red *Mesh_RedAdd(tMesh_Mesh *,tMesh_VolumeElement *);
  void Mesh_RedDelete(tMesh_Mesh *,tMesh_Red *); */

  tMesh_Mesh *Mesh_Read(FILE *,FILE *);
  void Mesh_Write(tMesh_Mesh *,FILE *,FILE *);

  void Mesh_FindNeighbours1(tMesh_VolumeElement *);
  void Mesh_FindNeighbours(tMesh_Mesh *);


  #define Mesh_ContainsPoint(ve,p) ((ve)->P[0]==(p)||(ve)->P[1]==(p)||\
    (ve)->P[2]==(p))
  #define Mesh_ContainsEdge(ve,p1,p2) (Mesh_ContainsPoint(ve,p1)&&\
    Mesh_ContainsPoint(ve,p2))
  #define Mesh_EdgeOpposite(ve,p1,p2)\
    ((tMesh_Point *)((pointer)(ve)->P[0]+(pointer)(ve)->P[1]+\
    (pointer)(ve)->P[2]-(pointer)(p1)-(pointer)(p2)))

  #define Mesh_SurfaceContainsPoint(se,p) ((se)->P[0]==(p)||(se)->P[1]==(p)||\
    (se)->P[2]==(p))
  #define Mesh_SurfaceContainsEdge(se,p1,p2)\
    (Mesh_SurfaceContainsPoint(se,p1)&&Mesh_SurfaceContainsPoint(se,p2))
  #define Mesh_SurfaceIsFace(se,p1,p2,p3)\
    (Mesh_SurfaceContainsPoint(se,p1)&&Mesh_SurfaceContainsPoint(se,p2)&&\
    Mesh_SurfaceContainsPoint(se,p3))
  #define Mesh_EdgeOpposite1(se,p1,p2) ((tMesh_Point *)((pointer)(se)->P[0]+\
    (pointer)(se)->P[1]+(pointer)(se)->P[2]-(pointer)(p1)-(pointer)(p2)))

#endif

