/usr/include/gmsh/MHexahedron.h

// Gmsh - Copyright (C) 1997-2015 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 _MHEXAHEDRON_H_
#define _MHEXAHEDRON_H_

#include "MElement.h"

/*
 * MHexahedron
 *
 *          v
 *   3----------2
 *   |\     ^   |\
 *   | \    |   | \
 *   |  \   |   |  \
 *   |   7------+---6
 *   |   |  +-- |-- | -> u
 *   0---+---\--1   |
 *    \  |    \  \  |
 *     \ |     \  \ |
 *      \|      w  \|
 *       4----------5
 *
 */
class MHexahedron : public MElement {
 protected:
  MVertex *_v[8];
  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[edges_hexa(num, 0)];
    v[1] = _v[edges_hexa(num, 1)];
  }
  void _getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[faces_hexa(num, 0)];
    v[1] = _v[faces_hexa(num, 1)];
    v[2] = _v[faces_hexa(num, 2)];
    v[3] = _v[faces_hexa(num, 3)];
  }
 public :
  MHexahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
              MVertex *v5, MVertex *v6, MVertex *v7, 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; _v[6] = v6; _v[7] = v7;
  }
  MHexahedron(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MElement(num, part)
  {
    for(int i = 0; i < 8; i++) _v[i] = v[i];
  }
  ~MHexahedron(){}
  virtual int getDim() const { return 3; }
  virtual int getNumVertices() const { return 8; }
  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 const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
  virtual MVertex *getVertexDIFF(int num)
  {
    static const int map[8] = {2, 3, 7, 6, 0, 1, 5, 4};
    return getVertex(map[num]);
  }
  virtual int getNumEdges()const{ return 12; }
  virtual MEdge getEdge(int num) const
  {
    return MEdge(_v[edges_hexa(num, 0)], _v[edges_hexa(num, 1)]);
  }
  virtual int getNumEdgesRep(bool curved){ return 12; }
  virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
  {
    static const int f[12] = {0, 0, 1, 0, 1, 0, 3, 2, 1, 2, 3, 4};
    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 6; }
  virtual MFace getFace(int num)
  {
    return MFace(_v[faces_hexa(num, 0)],
                 _v[faces_hexa(num, 1)],
                 _v[faces_hexa(num, 2)],
                 _v[faces_hexa(num, 3)]);
  }
  virtual double getInnerRadius();
  virtual double angleShapeMeasure();
  virtual void getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const;
  virtual int getNumFacesRep(bool curved){ return 12; }
  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
  {
    static const int f[12][3] = {
      {0, 3, 2}, {0, 2, 1},
      {0, 1, 5}, {0, 5, 4},
      {0, 4, 7}, {0, 7, 3},
      {1, 2, 6}, {1, 6, 5},
      {2, 3, 7}, {2, 7, 6},
      {4, 5, 6}, {4, 6, 7}
    };
    _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(4);
    _getFaceVertices(num, v);
  }
  virtual int getType() const { return TYPE_HEX; }
  virtual int getTypeForMSH() const { return MSH_HEX_8; }
  virtual int getTypeForUNV() const { return 115; } // solid linear brick
  virtual int getTypeForVTK() const { return 12; }
  virtual const char *getStringForPOS() const { return "SH"; }
  virtual const char *getStringForBDF() const { return "CHEXA"; }
  virtual const char *getStringForDIFF() const { return "ElmB8n3D"; }
  virtual const char *getStringForINP() const { return "C3D8"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
    tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
  }
  virtual int getVolumeSign();
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    switch(num) {
    case 0 : u = -1.; v = -1.; w = -1.; break;
    case 1 : u =  1.; v = -1.; w = -1.; break;
    case 2 : u =  1.; v =  1.; w = -1.; break;
    case 3 : u = -1.; v =  1.; w = -1.; break;
    case 4 : u = -1.; v = -1.; w =  1.; break;
    case 5 : u =  1.; v = -1.; w =  1.; break;
    case 6 : u =  1.; v =  1.; w =  1.; break;
    case 7 : u = -1.; v =  1.; w =  1.; break;
    default: u =  0.; v =  0.; w =  0.; break;
    }
  }
  virtual SPoint3 barycenterUVW() const
  {
    return SPoint3(0., 0., 0.);
  }
  virtual bool isInside(double u, double v, double w) const
  {
    double tol = _isInsideTolerance;
    if(u < -(1. + tol) || v < -(1. + tol) || w < -(1. + tol) ||
       u > (1. + tol) || v > (1. + tol) || w > (1. + tol))
      return false;
    return true;
  }
  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
  static int edges_hexa(const int edge, const int vert)
  {
    static const int e[12][2] = {
      {0, 1},
      {0, 3},
      {0, 4},
      {1, 2},
      {1, 5},
      {2, 3},
      {2, 6},
      {3, 7},
      {4, 5},
      {4, 7},
      {5, 6},
      {6, 7}
    };
    return e[edge][vert];
  }
  static int faces_hexa(const int face, const int vert)
  {
    static const int f[6][4] = {
      {0, 3, 2, 1},
      {0, 1, 5, 4},
      {0, 4, 7, 3},
      {1, 2, 6, 5},
      {2, 3, 7, 6},
      {4, 5, 6, 7}
    };
    return f[face][vert];
  }
  static int faces2edge_hexa(const int face, const int edge)
  {
    // return -iedge - 1 if edge is inverted
    //         iedge + 1 otherwise
    static const int e[6][4] = {
      {2, -6,  -4,  -1},
      {1,  5,  -9,  -3},
      {3, 10,  -8,  -2},
      {4,  7, -11,  -5},
      {6,  8, -12,  -7},
      {9, 11,  12, -10}
    };
    return e[face][edge];
  }
};

/*
 * MHexahedron20
 *
 *   3----13----2
 *   |\         |\
 *   | 15       | 14
 *   9  \       11 \
 *   |   7----19+---6
 *   |   |      |   |
 *   0---+-8----1   |
 *    \  17      \  18
 *    10 |        12|
 *      \|         \|
 *       4----16----5
 *
 */
class MHexahedron20 : public MHexahedron {
 protected:
  MVertex *_vs[12];
 public :
  MHexahedron20(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, MVertex *v18, MVertex *v19,
                int num=0, int part=0)
    : MHexahedron(v0, v1, v2, v3, v4, v5, v6, v7, num, part)
  {
    _vs[0] = v8; _vs[1] = v9; _vs[2] = v10; _vs[3] = v11; _vs[4] = v12;
    _vs[5] = v13; _vs[6] = v14; _vs[7] = v15; _vs[8] = v16; _vs[9] = v17;
    _vs[10] = v18; _vs[11] = v19;
    for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
  }
  MHexahedron20(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MHexahedron(v, num, part)
  {
    for(int i = 0; i < 12; i++) _vs[i] = v[8 + i];
    for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MHexahedron20(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 20; }
  virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
  virtual MVertex *getVertexUNV(int num)
  {
    static const int map[20] = {0, 8, 1, 11, 2, 13, 3, 9, 10, 12,
                                14, 15, 4, 16, 5, 18, 6, 19, 7, 17};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexBDF(int num)
  {
    static const int map[20] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
                                9, 10, 12, 14, 15, 16, 18, 19, 17};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexINP(int num)
  {
    static const int map[20]={0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
                              9, 16, 18, 19, 17, 10, 12, 14, 15};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexDIFF(int num)
  {
    static const int map[20] = {2, 3, 7, 6, 0, 1, 5, 4, 9, 18, 12,
                                19, 14, 11, 15, 13, 8, 16, 17, 10};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexVTK(int num)
  {
    static const int map[20] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
                                9, 16, 18, 19, 17, 10, 12, 14, 15};
    return getVertex(map[num]);
  }
  virtual int getNumEdgeVertices() const { return 12; }
  virtual int getNumEdgesRep(bool curved);
  virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    MHexahedron::_getEdgeVertices(num, v);
    v[2] = _vs[num];
  }
  virtual int getNumFacesRep(bool curved);
  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(8);
    MHexahedron::_getFaceVertices(num, v);
    static const int f[6][4] = {
      {1,  5,  3, 0},
      {0,  4,  8, 2},
      {2,  9,  7, 1},
      {3,  6, 10, 4},
      {5,  7, 11, 6},
      {8, 10, 11, 9}
    };
    v[4] = _vs[f[num][0]];
    v[5] = _vs[f[num][1]];
    v[6] = _vs[f[num][2]];
    v[7] = _vs[f[num][3]];
  }
  virtual int getTypeForMSH() const { return MSH_HEX_20; }
  virtual int getTypeForUNV() const { return 116; } // solid parabolic brick
  virtual int getTypeForVTK() const { return 25; }
  virtual const char *getStringForBDF() const { return "CHEXA"; }
  virtual const char *getStringForINP() const { return "C3D20"; }
  virtual const char *getStringForDIFF() const { return "ElmB20n3D"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
    tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
    MVertex *old[12];
    for(int i = 0; i < 12; i++) old[i] = _vs[i];
    _vs[0] = old[3]; _vs[3] = old[0];
    _vs[1] = old[5]; _vs[5] = old[1];
    _vs[2] = old[6]; _vs[6] = old[2];
    _vs[8] = old[10]; _vs[10] = old[8];
    _vs[9] = old[11]; _vs[11] = old[9];
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

/*
 * MHexahedron27
 *
 *   3----13----2
 *   |\         |\
 *   |15    24  | 14
 *   9  \ 20    11 \
 *   |   7----19+---6
 *   |22 |  26  | 23|
 *   0---+-8----1   |
 *    \ 17    25 \  18
 *    10 |  21    12|
 *      \|         \|
 *       4----16----5
 *
 */
class MHexahedron27 : public MHexahedron {
 protected:
  MVertex *_vs[19];
 public :
  MHexahedron27(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, MVertex *v18, MVertex *v19,
                MVertex *v20, MVertex *v21, MVertex *v22, MVertex *v23, MVertex *v24,
                MVertex *v25, MVertex *v26, int num=0, int part=0)
    : MHexahedron(v0, v1, v2, v3, v4, v5, v6, v7, num, part)
  {
    _vs[0] = v8; _vs[1] = v9; _vs[2] = v10; _vs[3] = v11; _vs[4] = v12;
    _vs[5] = v13; _vs[6] = v14; _vs[7] = v15; _vs[8] = v16; _vs[9] = v17;
    _vs[10] = v18; _vs[11] = v19; _vs[12] = v20; _vs[13] = v21; _vs[14] = v22;
    _vs[15] = v23; _vs[16] = v24; _vs[17] = v25; _vs[18] = v26;
    for(int i = 0; i < 19; i++) _vs[i]->setPolynomialOrder(2);
  }
  MHexahedron27(const std::vector<MVertex*> &v, int num=0, int part=0)
    : MHexahedron(v, num, part)
  {
    for(int i = 0; i < 19; i++) _vs[i] = v[8 + i];
    for(int i = 0; i < 19; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MHexahedron27(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 27; }
  virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }

  virtual MVertex *getVertexDIFF(int num)
  {
    static const int map[27] = {6, 8, 26, 24, 0, 2, 20, 18, 7, 15, 3, 17, 5, 25,
                                23, 21, 1, 9, 11, 19, 16, 4, 12, 14, 22, 10, 13};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexVTK(int num)
  {
    static const int map[27] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13, 9, 16, 18, 19,
                                17, 10, 12, 14,15, 22, 23, 21, 24, 20, 25, 26};
    return getVertex(map[num]);
  }
  virtual int getNumEdgeVertices() const { return 12; }
  virtual int getNumFaceVertices() const { return 6; }
  virtual int getNumVolumeVertices() const { return 1; }
  virtual int getNumEdgesRep(bool curved);
  virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    MHexahedron::_getEdgeVertices(num, v);
    v[2] = _vs[num];
  }
  virtual int getNumFacesRep(bool curved);
  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(9);
    MHexahedron::_getFaceVertices(num, v);
    static const int f[6][4] = {
      {1,  5,  3, 0},
      {0,  4,  8, 2},
      {2,  9,  7, 1},
      {3,  6, 10, 4},
      {5,  7, 11, 6},
      {8, 10, 11, 9}
    };
    v[4] = _vs[f[num][0]];
    v[5] = _vs[f[num][1]];
    v[6] = _vs[f[num][2]];
    v[7] = _vs[f[num][3]];
    v[8] = _vs[12+num];
  }
  virtual int getTypeForMSH() const { return MSH_HEX_27; }
  virtual int getTypeForVTK() const { return 29; }
  virtual const char *getStringForPOS() const { return "SH2"; }
  virtual const char *getStringForDIFF() const { return "ElmB27n3D"; }
  virtual void reverse()
  {
    MVertex *tmp;
    tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
    tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
    MVertex *old[19];
    for(int i = 0; i < 19; i++) old[i] = _vs[i];
    // edge vertices
    _vs[0] = old[3]; _vs[3] = old[0];
    _vs[1] = old[5]; _vs[5] = old[1];
    _vs[2] = old[6]; _vs[6] = old[2];
    _vs[8] = old[10]; _vs[10] = old[8];
    _vs[9] = old[11]; _vs[11] = old[9];
    // face vertices
    _vs[13] = old[15]; _vs[15] = old[13];
    _vs[14] = old[16]; _vs[16] = old[14];
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

/*
 * MHexahedronN
 *
 *   3----13----2
 *   |\         |\
 *   |15    24  | 14
 *   9  \ 20    11 \
 *   |   7----19+---6
 *   |22 |  26  | 23|
 *   0---+-8----1   |
 *    \ 17    25 \  18
 *    10 |  21    12|
 *      \|         \|
 *       4----16----5
 *
 *    N = 0RDER - 1
 *    8 := 8 --> 8 + N
 *    9 := 8 + N + 1 --> 8 +  2 * N
 *    :
 *   19 := 8 + 11 * N + 1 --> 8 +  12 * N
 *   20 := 8 + 12 * N + 1 --> 8 +  12 * N + N^2
 *   21 := 8 + 12 * N + N^2 --> 8 +  12 * N + 2 * N^2
 *    :
 *   25 := 8 + 12 * N + 5 * N^2 + 1 --> 8 +  12 * N + 6 * N^2
 *   26 := 8 + 12 * N + 6 * N^2 + 1 --> 8 +  12 * N + 6 * N^2 + N^3
 *
 */


class MHexahedronN : public MHexahedron {
 protected:
  const char _order;
  std::vector<MVertex*> _vs;
 public :
 MHexahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
	      MVertex *v4, MVertex *v5, MVertex *v6, MVertex *v7,
	      const std::vector<MVertex*> &v, char order, int num=0, int part=0)
   : MHexahedron(v0, v1, v2, v3, v4, v5, v6, v7, num, part), _order(order), _vs(v)
  {
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
 MHexahedronN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
   :  MHexahedron(v[0], v[1], v[2], v[3],v[4], v[5], v[6], v[7], num, part), _order(order)
  {
    for(unsigned int i = 8; i < v.size(); i++) _vs.push_back(v[i]);
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
  ~MHexahedronN(){}
  virtual int getPolynomialOrder() const { return (int)_order; }
  virtual int getNumVertices() const { return 8 + _vs.size(); }
  virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
  virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }

  virtual int getNumEdgeVertices() const { return 12 * (_order - 1); }
  virtual int getNumFaceVertices() const
  {
    if (getIsAssimilatedSerendipity())
      return 0;
    else
      return 6 * (_order - 1)*(_order - 1);
  }
  virtual int getNumVolumeVertices() const
  {
    if (getIsAssimilatedSerendipity())
      return 0;
    else
      return (_order - 1) * (_order - 1) * (_order - 1);
  }
  virtual int getNumEdgesRep(bool curved);
  virtual void getEdgeRep(bool curved, 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);
    MHexahedron::_getEdgeVertices(num, v);
    for (int i = 0; i < _order - 1; i++) v[2+i] = _vs[(_order-1)*num+i];
  }
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    if (getIsAssimilatedSerendipity()) {
      v.resize(4 * _order);
    }
    else {
      v.resize((_order+1) * (_order+1));
    }

    MHexahedron::_getFaceVertices(num, v);
    int count = 3;

    int n = _order-1;
    for (int i = 0; i < 4; i++) {
      if(faces2edge_hexa(num, i) > 0){
        int edge_num = faces2edge_hexa(num, i) - 1;
        for (int j = 0; j < n; j++) v[++count] = _vs[n*edge_num + j];
      }
      else{
        int edge_num = -faces2edge_hexa(num, i) - 1;
        for (int j = n-1; j >= 0; j--) v[++count] = _vs[n*edge_num + j];
      }
    }

    if ((int)v.size() > count + 1) {
      int start = 12 * n + num * n*n;
      for (int i = 0; i < n*n; i++){
        v[++count] = _vs[start + i];
      }
    }
  }
  virtual int getTypeForMSH() const
  {
    // (p+1)^3
    if(_order == 1 && _vs.size() + 8 == 8)   return MSH_HEX_8;
    if(_order == 2 && _vs.size() + 8 == 27)  return MSH_HEX_27;
    if(_order == 3 && _vs.size() + 8 == 64 ) return MSH_HEX_64;
    if(_order == 4 && _vs.size() + 8 == 125) return MSH_HEX_125;
    if(_order == 5 && _vs.size() + 8 == 216) return MSH_HEX_216;
    if(_order == 6 && _vs.size() + 8 == 343) return MSH_HEX_343;
    if(_order == 7 && _vs.size() + 8 == 512) return MSH_HEX_512;
    if(_order == 8 && _vs.size() + 8 == 729) return MSH_HEX_729;
    if(_order == 9 && _vs.size() + 8 == 1000) return MSH_HEX_1000;

    if(_order == 2 && _vs.size() + 8 == 20)  return MSH_HEX_20;
    if(_order == 3 && _vs.size() + 8 == 32)  return MSH_HEX_32;
    if(_order == 4 && _vs.size() + 8 == 44)  return MSH_HEX_44;
    if(_order == 5 && _vs.size() + 8 == 56)  return MSH_HEX_56;
    if(_order == 6 && _vs.size() + 8 == 68)  return MSH_HEX_68;
    if(_order == 7 && _vs.size() + 8 == 80)  return MSH_HEX_80;
    if(_order == 8 && _vs.size() + 8 == 92)  return MSH_HEX_92;
    if(_order == 9 && _vs.size() + 8 == 104) return MSH_HEX_104;

    Msg::Error("no tag matches a p%d hexahedron with %d vertices", _order, 8+_vs.size());
    return 0;
  }
  virtual int getNumFacesRep(bool curved);
  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
  virtual void reverse()
  {
    Msg::Error("Reverse not implemented yet for MHexahedronN");
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

#endif
libgmsh-dev 2.10.1+dfsg1-1ubuntu4 / usr / include / gmsh / MHexahedron.h
