// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#ifndef _MTRIANGLE_H_
#define _MTRIANGLE_H_

#include "MElement.h"

/* 
 * MTriangle
 *
 *   v
 *   ^
 *   |
 *   2
 *   |`\
 *   |  `\
 *   |    `\
 *   |      `\
 *   |        `\
 *   0----------1 --> u
 *
 */
class MTriangle : public MElement {
 protected:
  MVertex *_v[3];
  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[edges_tri(num, 0)];
    v[1] = _v[edges_tri(num, 1)];
  }
  void _getFaceVertices(std::vector<MVertex*> &v) const
  {
    v[0] = _v[0];
    v[1] = _v[1];
    v[2] = _v[2];
  }
 public :
  MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0) 
    : MElement(num, part)
  {
    _v[0] = v0; _v[1] = v1; _v[2] = v2;
  }
  MTriangle(std::vector<MVertex*> &v, int num=0, int part=0) 
    : MElement(num, part)
  {
    for(int i = 0; i < 3; i++) _v[i] = v[i];
  }
  ~MTriangle(){}
  virtual int getDim(){ return 2; }
  virtual double gammaShapeMeasure();
  virtual double distoShapeMeasure();
  virtual int getNumVertices() const { return 3; }
  virtual MVertex *getVertex(int num){ return _v[num]; }
  virtual MVertex *getVertexMED(int num)
  {
    static const int map[3] = {0, 2, 1};
    return getVertex(map[num]); 
  }
  virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
  { 
    if(_v[0] != v1 && _v[0] != v2) return _v[0];
    if(_v[1] != v1 && _v[1] != v2) return _v[1];
    if(_v[2] != v1 && _v[2] != v2) return _v[2];
    return 0;
  }
  virtual int getNumEdges(){ return 3; }
  virtual MEdge getEdge(int num)
  {
    return MEdge(_v[edges_tri(num, 0)], _v[edges_tri(num, 1)]);
  }
  virtual int getNumEdgesRep(){ return 3; }
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
  { 
    MEdge e(getEdge(num));
    _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, 0);
  }
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(2);
    _getEdgeVertices(num, v);
  }
  virtual int getNumFaces(){ return 1; }
  virtual MFace getFace(int num)
  { 
    return MFace(_v[0], _v[1], _v[2]); 
  }
  virtual int getNumFacesRep(){ return 1; }
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
  { 
    _getFaceRep(_v[0], _v[1], _v[2], x, y, z, n);
  }
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    _getFaceVertices(v);
  }
  virtual int getType() const { return TYPE_TRI; }
  virtual int getTypeForMSH() const { return MSH_TRI_3; }
  virtual int getTypeForUNV() const { return 91; } // thin shell linear triangle
  virtual int getTypeForVTK() const { return 5; }
  virtual const char *getStringForPOS() const { return "ST"; }
  virtual const char *getStringForBDF() const { return "CTRIA3"; }
  virtual const char *getStringForDIFF() const { return "ElmT3n2D"; }
  virtual void revert() 
  {
    MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
  }
  virtual const functionSpace* getFunctionSpace(int o=-1) const;
  virtual bool isInside(double u, double v, double w)
  {
    double tol = _isInsideTolerance;
    if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
      return false; 
    return true;
  }
  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
  virtual SPoint3 circumcenter();
 private:
  int edges_tri(const int edge, const int vert) const
  {
    static const int e[3][2] = {
      {0, 1},
      {1, 2},
      {2, 0}
    };
    return e[edge][vert];
  }
};

/*
 * MTriangle6
 *
 *   2
 *   |`\
 *   |  `\
 *   5    `4
 *   |      `\
 *   |        `\
 *   0-----3----1
 *
 */
class MTriangle6 : public MTriangle {
 protected:
  MVertex *_vs[3];
 public :
  MTriangle6(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
             MVertex *v5, int num=0, int part=0) 
    : MTriangle(v0, v1, v2, num, part)
  {
    _vs[0] = v3; _vs[1] = v4; _vs[2] = v5;
    for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
  }
  MTriangle6(std::vector<MVertex*> &v, int num=0, int part=0) 
    : MTriangle(v, num, part)
  {
    for(int i = 0; i < 3; i++) _vs[i] = v[3 + i];
    for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MTriangle6(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 6; }
  virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
  virtual MVertex *getVertexUNV(int num)
  {
    static const int map[6] = {0, 3, 1, 4, 2, 5};
    return getVertex(map[num]); 
  }
  virtual MVertex *getVertexMED(int num)
  {
    static const int map[6] = {0, 2, 1, 5, 4, 3};
    return getVertex(map[num]); 
  }
  virtual int getNumEdgeVertices() const { return 3; }
  virtual int getNumEdgesRep();
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    MTriangle::_getEdgeVertices(num, v);
    v[2] = _vs[num];
  }
  virtual int getNumFacesRep();
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(6);
    MTriangle::_getFaceVertices(v);
    v[3] = _vs[0];
    v[4] = _vs[1];
    v[5] = _vs[2];
  }
  virtual int getTypeForMSH() const { return MSH_TRI_6; }
  virtual int getTypeForUNV() const { return 92; } // thin shell parabolic triangle
  //virtual int getTypeForVTK() const { return 22; }
  virtual const char *getStringForPOS() const { return "ST2"; }
  virtual const char *getStringForBDF() const { return "CTRIA6"; }
  virtual const char *getStringForDIFF() const { return "ElmT6n2D"; }
  virtual void revert() 
  {
    MVertex *tmp;
    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
    tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
  }
};

/*
 * MTriangleN  FIXME: check the plot
 *
 *   2
 *   |`\                E = order - 1;
 *   |  `\              N = total number of vertices
 * 3+2E   2+2E
 *   |      `\          Interior vertex numbers
 *  ...       ...         for edge 0 <= i <= 2: 3+i*E to 2+(i+1)*E
 *   |          `\        in volume           : 3+3*E to N-1
 * 2+3E           3+E
 *   |  3+3E to N-1 `\
 *   |                `\
 *   0---3--...---2+E---1
 *
 */
class MTriangleN : public MTriangle {
 protected:
  std::vector<MVertex *> _vs;
  const char _order;
 public:
  MTriangleN(MVertex *v0, MVertex *v1, MVertex *v2, 
             std::vector<MVertex*> &v, char order, int num=0, int part=0) 
    : MTriangle(v0, v1, v2, num, part), _vs(v), _order(order)
  {
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
  MTriangleN(std::vector<MVertex*> &v, char order, int num=0, int part=0) 
    : MTriangle(v[0], v[1], v[2], num, part), _order(order)
  {
    for(unsigned int i = 3; i < v.size(); i++) _vs.push_back(v[i]);
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
  ~MTriangleN(){}
  virtual int getPolynomialOrder() const { return _order; }
  virtual int getNumVertices() const { return 3 + _vs.size(); }
  virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
  virtual int getNumFaceVertices() const 
  {
    if(_order == 3 && _vs.size() == 6) return 0;
    if(_order == 3 && _vs.size() == 7) return 1;
    if(_order == 4 && _vs.size() == 9) return 0;
    if(_order == 4 && _vs.size() == 12) return 3;
    if(_order == 5 && _vs.size() == 12) return 0;
    if(_order == 5 && _vs.size() == 18) return 6;
    return 0;
  }
  virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
  virtual int getNumEdgesRep();
  virtual int getNumFacesRep();
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(_order + 1);
    MTriangle::_getEdgeVertices(num, v);
    int j = 2;
    const int ie = (num + 1) * (_order - 1);
    for(int i = num * (_order-1); i != ie; ++i) v[j++] = _vs[i];
  }
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3 + _vs.size());
    MTriangle::_getFaceVertices(v);
    for(unsigned int i = 0; i != _vs.size(); ++i) v[i + 3] = _vs[i];
  }
  virtual int getTypeForMSH() const
  {
    if(_order == 2 && _vs.size() == 3) return MSH_TRI_6; 
    if(_order == 3 && _vs.size() == 6) return MSH_TRI_9; 
    if(_order == 3 && _vs.size() == 7) return MSH_TRI_10; 
    if(_order == 4 && _vs.size() == 9) return MSH_TRI_12; 
    if(_order == 4 && _vs.size() == 12) return MSH_TRI_15; 
    if(_order == 5 && _vs.size() == 12) return MSH_TRI_15I; 
    if(_order == 5 && _vs.size() == 18) return MSH_TRI_21;
    return 0;
  }
  virtual void revert() 
  {
    MVertex *tmp;
    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
    std::vector<MVertex*> inv;
    inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
    _vs = inv;
  }
};

template <class T> 
void sort3(T *t[3])
{
  T *temp;
  if(t[0] > t[1]){
    temp = t[1];
    t[1] = t[0];
    t[0] = temp;
  }
  if(t[1] > t[2]){
    temp = t[2];
    t[2] = t[1];
    t[1] = temp;
  }
  if(t[0] > t[1]){
    temp = t[1];
    t[1] = t[0];
    t[0] = temp;
  }
}

struct compareMTriangleLexicographic
{
  bool operator () (MTriangle *t1, MTriangle *t2) const
  {
    MVertex *_v1[3] = {t1->getVertex(0), t1->getVertex(1), t1->getVertex(2)};
    MVertex *_v2[3] = {t2->getVertex(0), t2->getVertex(1), t2->getVertex(2)};
    sort3(_v1);
    sort3(_v2);
    if(_v1[0] < _v2[0]) return true;
    if(_v1[0] > _v2[0]) return false;
    if(_v1[1] < _v2[1]) return true;
    if(_v1[1] > _v2[1]) return false;
    if(_v1[2] < _v2[2]) return true;
    return false;
  }
};

#endif
