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

#ifndef _MPRISM_H_
#define _MPRISM_H_

#include "MElement.h"

/*
 * MPrism
 *
 *               w
 *               ^
 *               |
 *               3
 *             ,/|`\
 *           ,/  |  `\
 *         ,/    |    `\
 *        4------+------5
 *        |      |      |
 *        |    ,/|`\    |
 *        |  ,/  |  `\  |
 *        |,/    |    `\|
 *       ,|      |      `\
 *     ,/ |      0      | `\
 *    u   |    ,/ `\    |    v
 *        |  ,/     `\  |
 *        |,/         `\|
 *        1-------------2
 *
 */
class MPrism : public MElement {
 protected:
  MVertex *_v[6];
  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[edges_prism(num, 0)];
    v[1] = _v[edges_prism(num, 1)];
  }
  void _getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[faces_prism(num, 0)];
    v[1] = _v[faces_prism(num, 1)];
    v[2] = _v[faces_prism(num, 2)];
    if(num >= 2) v[3] = _v[faces_prism(num, 3)];
  }
 public :
  MPrism(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
         MVertex *v5, int num=0, int part=0)
    : MElement(num, part)
  {
    _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
    _v[4] = v4; _v[5] = v5;
  }
  MPrism(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MElement(num, part)
  {
    for(int i = 0; i < 6; i++) _v[i] = v[i];
  }
  ~MPrism(){}
  virtual int getDim() const { return 3; }
  virtual int getNumVertices() const { return 6; }
  virtual double getInnerRadius();
  virtual MVertex *getVertex(int num){ return _v[num]; }
  virtual const MVertex *getVertex(int num)const{ return _v[num]; }
  virtual void setVertex(int num,   MVertex *v){ _v[num] = v; }
  virtual int getNumEdges(){ return 9; }
  virtual MEdge getEdge(int num) const
  {
    return MEdge(_v[edges_prism(num, 0)], _v[edges_prism(num, 1)]);
  }
  virtual int getNumEdgesRep(){ return 9; }
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
  {
    static const int f[9] = {0, 1, 2, 0, 2, 3, 1, 1, 1};
    MEdge e(getEdge(num));
    _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, f[num]);
  }
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(2);
    _getEdgeVertices(num, v);
  }
  virtual int getNumFaces(){ return 5; }
  virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign, int &rot) const;
  virtual MFace getFace(int num)
  {
    if(num < 2)
      return MFace(_v[faces_prism(num, 0)],
                   _v[faces_prism(num, 1)],
                   _v[faces_prism(num, 2)]);
    else
      return MFace(_v[faces_prism(num, 0)],
                   _v[faces_prism(num, 1)],
                   _v[faces_prism(num, 2)],
                   _v[faces_prism(num, 3)]);
  }
  virtual int getNumFacesRep(){ return 8; }
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
  {
    static const int f[8][3] = {
      {0, 2, 1},
      {3, 4, 5},
      {0, 1, 4}, {0, 4, 3},
      {0, 3, 5}, {0, 5, 2},
      {1, 2, 5}, {1, 5, 4}
    };
    _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
                x, y, z, n);
  }
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize((num < 2) ? 3 : 4);
    _getFaceVertices(num, v);
  }
  virtual int getType() const { return TYPE_PRI; }
  virtual int getTypeForMSH() const { return MSH_PRI_6; }
  virtual int getTypeForUNV() const { return 112; } // solid linear wedge
  virtual int getTypeForVTK() const { return 13; }
  virtual const char *getStringForPOS() const { return "SI"; }
  virtual const char *getStringForBDF() const { return "CPENTA"; }
  virtual const char *getStringForINP() const { return "C3D6"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
    tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
  }
  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
  virtual int getVolumeSign();
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    switch(num) {
    case 0 : u = 0.; v = 0.; w = -1.; break;
    case 1 : u = 1.; v = 0.; w = -1.; break;
    case 2 : u = 0.; v = 1.; w = -1.; break;
    case 3 : u = 0.; v = 0.; w =  1.; break;
    case 4 : u = 1.; v = 0.; w =  1.; break;
    case 5 : u = 0.; v = 1.; w =  1.; break;
    default: u = 0.; v = 0.; w =  0.; break;
    }
  }
  virtual SPoint3 barycenterUVW() const
  {
    return SPoint3(1/3., 1/3., 0.);
  }
  virtual bool isInside(double u, double v, double w) const
  {
    double tol = _isInsideTolerance;
    if(w > (1. + tol) || w < -(1. + tol) || u < (-tol)
       || v < (-tol) || u > ((1. + tol) - v))
      return false;
    return true;
  }
  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
  static int edges_prism(const int edge, const int vert)
  {
    static const int e[9][2] = {
      {0, 1},
      {0, 2},
      {0, 3},
      {1, 2},
      {1, 4},
      {2, 5},
      {3, 4},
      {3, 5},
      {4, 5}
    };
    return e[edge][vert];
  }
  static int faces_prism(const int face, const int vert)
  {
    static const int f[5][4] = {
      {0, 2, 1, -1},
      {3, 4, 5, -1},
      {0, 1, 4,  3},
      {0, 3, 5,  2},
      {1, 2, 5,  4}
    };
    return f[face][vert];
  }
};

/*
 * MPrism15
 *
 *               3
 *             ,/|`\
 *           12  |  13
 *         ,/    |    `\
 *        4------14-----5
 *        |      8      |
 *        |      |      |
 *        |      |      |
 *        |      |      |
 *        10     |      11
 *        |      0      |
 *        |    ,/ `\    |
 *        |  ,6     `7  |
 *        |,/         `\|
 *        1------9------2
 *
 */
class MPrism15 : public MPrism {
 protected:
  MVertex *_vs[9];
 public :
  MPrism15(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
           MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
           MVertex *v10, MVertex *v11, MVertex *v12, MVertex *v13, MVertex *v14,
           int num=0, int part=0)
    : MPrism(v0, v1, v2, v3, v4, v5, num, part)
  {
    _vs[0] = v6; _vs[1] = v7; _vs[2] = v8; _vs[3] = v9; _vs[4] = v10;
    _vs[5] = v11; _vs[6] = v12; _vs[7] = v13; _vs[8] = v14;
    for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
  }
  MPrism15(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MPrism(v, num, part)
  {
    for(int i = 0; i < 9; i++) _vs[i] = v[6 + i];
    for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MPrism15(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 15; }
  virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
  virtual const MVertex *getVertex(int num) const { return num < 6 ? _v[num] : _vs[num - 6]; }
  virtual MVertex *getVertexUNV(int num)
  {
    static const int map[15] = {0, 6, 1, 9, 2, 7, 8, 10, 11, 3, 12, 4, 14, 5, 13};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexBDF(int num)
  {
    static const int map[15] = {0, 1, 2, 3, 4, 5, 6, 9, 7, 8, 10, 11, 12, 14, 13};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexINP(int num){ return getVertexBDF(num); }
  virtual int getNumEdgeVertices() const { return 9; }
  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);
    MPrism::_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((num < 2) ? 6 : 8);
    MPrism::_getFaceVertices(num, v);
    static const int f[5][4] = {
      {1, 3, 0, -1},
      {6, 8, 7, -1},
      {0, 4, 6,  2},
      {2, 7, 5,  1},
      {3, 5, 8,  4}
    };
    const int i = (num < 2) ? 3 : 4;
    v[i  ] = _vs[f[num][0]];
    v[i+1] = _vs[f[num][1]];
    v[i+2] = _vs[f[num][2]];
    if (num >= 2) v[7] = _vs[f[num][3]];
  }
  virtual int getTypeForMSH() const { return MSH_PRI_15; }
  virtual int getTypeForUNV() const { return 113; } // solid parabolic wedge
  virtual const char *getStringForBDF() const { return "CPENTA"; }
  virtual const char *getStringForINP() const { return "C3D15"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
    tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
    tmp = _vs[1]; _vs[1] = _vs[3]; _vs[3] = tmp;
    tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp;
    tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp;
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

/*
 * MPrism18
 *
 *               3
 *             ,/|`\
 *           12  |  13
 *         ,/    |    `\
 *        4------14-----5
 *        |      8      |
 *        |    ,/|`\    |
 *        |  15  |  16  |
 *        |,/    |    `\|
 *        10-----17-----11
 *        |      0      |
 *        |    ,/ `\    |
 *        |  ,6     `7  |
 *        |,/         `\|
 *        1------9------2
 *
 */
class MPrism18 : public MPrism {
 protected:
  MVertex *_vs[12];
 public :
  MPrism18(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
           MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
           MVertex *v10, MVertex *v11, MVertex *v12, MVertex *v13, MVertex *v14,
           MVertex *v15, MVertex *v16, MVertex *v17, int num=0, int part=0)
    : MPrism(v0, v1, v2, v3, v4, v5, num, part)
  {
    _vs[0] = v6; _vs[1] = v7; _vs[2] = v8; _vs[3] = v9; _vs[4] = v10;
    _vs[5] = v11; _vs[6] = v12; _vs[7] = v13; _vs[8] = v14; _vs[9] = v15;
    _vs[10] = v16; _vs[11] = v17;
    for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
  }
  MPrism18(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MPrism(v, num, part)
  {
    for(int i = 0; i < 12; i++) _vs[i] = v[6 + i];
    for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MPrism18(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 18; }
  virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
  virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num - 6]; }
  virtual int getNumEdgeVertices() const { return 9; }
  virtual int getNumFaceVertices() 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);
    MPrism::_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((num < 2) ? 6 : 9);
    MPrism::_getFaceVertices(num, v);
    static const int f[5][4] = {
      {1, 3, 0, -1},
      {6, 8, 7, -1},
      {0, 4, 6,  2},
      {2, 7, 5,  1},
      {3, 5, 8,  4}
    };
    const int i = (num < 2) ? 3 : 4;
    v[i  ] = _vs[f[num][0]];
    v[i+1] = _vs[f[num][1]];
    v[i+2] = _vs[f[num][2]];
    if (num >= 2) {
      v[7] = _vs[f[num][3]];
      v[8] = _vs[7+num];
    }
  }
  virtual int getTypeForMSH() const { return MSH_PRI_18; }
  virtual const char *getStringForPOS() const { return "SI2"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
    tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
    // edge vertices
    tmp = _vs[1]; _vs[1] = _vs[3]; _vs[3] = tmp;
    tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp;
    tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp;
    // quad face vertices
    tmp = _vs[10]; _vs[10] = _vs[11]; _vs[11] = tmp;
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

#endif
