From f258c78ffc5762267d166c599e295cdb60bae119 Mon Sep 17 00:00:00 2001
From: Stefen Guzik <guzik2@llnl.gov>
Date: Fri, 11 Jul 2008 01:21:02 +0000
Subject: [PATCH] -added routines getEdgeVertices and getFaceVertices for
 returning all vertices on a face or edge.  Includes some modification to
 getEdge and getFace too. -added ASCII representations of each element type.

---
 Geo/MElement.h | 1090 ++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 973 insertions(+), 117 deletions(-)

diff --git a/Geo/MElement.h b/Geo/MElement.h
index 81438026fe..cf4e67e472 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -126,6 +126,10 @@ class MElement
   virtual int getNumEdgesRep() = 0;
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) = 0;
 
+  // get all the vertices on an edge
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+    = 0;
+
   // get the faces
   virtual int getNumFaces() = 0;
   virtual MFace getFace(int num) = 0;
@@ -134,6 +138,12 @@ class MElement
   virtual int getNumFacesRep() = 0;
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n) = 0;
 
+  // get all the vertices on a face
+  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(0);
+  }
+
   // get the max/min edge length
   virtual double maxEdge();
   virtual double minEdge();
@@ -235,9 +245,24 @@ class MElementFactory{
   MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0);
 };
 
+/*
+ * MLine
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   0----------1 --> u
+ *
+ */
 class MLine : public MElement {
  protected:
   MVertex *_v[2];
+  void _getEdgeVertices(std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[0];
+    v[1] = _v[1];
+  }
  public :
   MLine(MVertex *v0, MVertex *v1, int num=0, int part=0) 
     : MElement(num, part)
@@ -260,6 +285,11 @@ class MLine : public MElement {
   { 
     _getEdgeRep(_v[0], _v[1], x, y, z, n);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(2);
+    _getEdgeVertices(v);
+  }
   virtual int getNumFaces(){ return 0; }
   virtual MFace getFace(int num){ return MFace(); }
   virtual int getNumFacesRep(){ return 0; }
@@ -298,6 +328,16 @@ class MLine : public MElement {
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
 };
 
+/*
+ * MLine3
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   0-----2----1 --> u
+ *
+ */
 class MLine3 : public MLine {
  protected:
   MVertex *_vs[1];
@@ -332,12 +372,28 @@ class MLine3 : public MLine {
     };
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MLine::_getEdgeVertices(v);
+    v[2] = _vs[0];
+  }
   virtual int getTypeForMSH(){ return MSH_LIN_3; }
   virtual int getTypeForUNV(){ return 24; } // parabolic beam
   virtual int getTypeForVTK(){ return 21; }
   virtual const char *getStringForPOS(){ return "SL2"; }
 };
 
+/*
+ * MLineN
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   0---2---...-(N-1)-1 --> u
+ *
+ */
 class MLineN : public MLine {
  protected:
   std::vector<MVertex *> _vs;
@@ -368,6 +424,12 @@ class MLineN : public MLine {
                 getVertex((num == getNumEdgesRep() - 1) ? 1 : num + 2),
                 x, y, z, n);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(2 + _vs.size());
+    MLine::_getEdgeVertices(v);
+    for(int i = 0; i != _vs.size(); ++i) v[i+2] = _vs[i];
+  }
   virtual int getTypeForMSH(){ 
     if(_vs.size() == 2) return MSH_LIN_4; 
     if(_vs.size() == 3) return MSH_LIN_5; 
@@ -376,9 +438,35 @@ class MLineN : public MLine {
   }
 };
 
+/* 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)
@@ -410,12 +498,7 @@ class MTriangle : public MElement {
   virtual int getNumEdges(){ return 3; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_tri[3][2] = {
-      {0, 1},
-      {1, 2},
-      {2, 0}
-    };
-    return MEdge(_v[edges_tri[num][0]], _v[edges_tri[num][1]]);
+    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)
@@ -423,6 +506,11 @@ class MTriangle : public MElement {
     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)
   { 
@@ -433,6 +521,11 @@ class MTriangle : public MElement {
   { 
     _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 getTypeForMSH(){ return MSH_TRI_3; }
   virtual int getTypeForUNV(){ return 91; } // thin shell linear triangle
   virtual int getTypeForVTK(){ return 5; }
@@ -477,8 +570,34 @@ class MTriangle : public MElement {
     return true;
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ 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
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   2
+ *   |`\
+ *   |  `\
+ *   5    `4
+ *   |      `\
+ *   |        `\
+ *   0-----3----1 --> u
+ *
+ */
 class MTriangle6 : public MTriangle {
  protected:
   MVertex *_vs[3];
@@ -521,6 +640,12 @@ class MTriangle6 : public MTriangle {
     }; 
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
   } 
+  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(){ return 4; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -530,6 +655,14 @@ class MTriangle6 : public MTriangle {
     _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(6);
+    MTriangle::_getFaceVertices(v);
+    v[3] = _vs[0];
+    v[4] = _vs[1];
+    v[5] = _vs[2];
+  }
   virtual int getTypeForMSH(){ return MSH_TRI_6; }
   virtual int getTypeForUNV(){ return 92; } // thin shell parabolic triangle
   virtual int getTypeForVTK(){ return 22; }
@@ -551,6 +684,26 @@ class MTriangle6 : public MTriangle {
   }
 };
 
+/*
+ * MTriangleN
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   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 --> u
+ *
+ */
 class MTriangleN : public MTriangle {
  protected:
   std::vector<MVertex *> _vs;
@@ -591,11 +744,25 @@ class MTriangleN : public MTriangle {
     if(_order == 5 && _vs.size() == 18) return 6;
     return 0;
   }
-  virtual int getNumEdgeVertices(){ return _order - 1; }
+  virtual int getNumEdgeVertices() const { return _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(2 + getNumEdgeVertices());
+    MTriangle::_getEdgeVertices(num, v);
+    int j = 2;
+    const int ie = (num+1)*getNumEdgeVertices();
+    for(int i = num*getNumEdgeVertices(); 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(int i = 0; i != _vs.size(); ++i) v[i+3] = _vs[i];
+  }
   virtual int getTypeForMSH()
   {
     if(_order == 3 && _vs.size() == 6) return MSH_TRI_9; 
@@ -624,9 +791,37 @@ class MTriangleN : public MTriangle {
   }
 };
 
+/*
+ * MQuadrangle
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3----------2
+ *   |          |
+ *   |          |
+ *   |          |
+ *   |          |
+ *   |          |
+ *   0----------1 --> u
+ *
+ */
 class MQuadrangle : public MElement {
  protected:
   MVertex *_v[4];
+  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[edges_quad(num, 0)];
+    v[1] = _v[edges_quad(num, 1)];
+  }
+  void _getFaceVertices(std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[0];
+    v[1] = _v[1];
+    v[2] = _v[2];
+    v[3] = _v[3];
+  }
  public :
   MQuadrangle(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0) 
     : MElement(num, part)
@@ -650,13 +845,7 @@ class MQuadrangle : public MElement {
   virtual int getNumEdges(){ return 4; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_quad[4][2] = {
-      {0, 1},
-      {1, 2},
-      {2, 3},
-      {3, 0}
-    };
-    return MEdge(_v[edges_quad[num][0]], _v[edges_quad[num][1]]);
+    return MEdge(_v[edges_quad(num, 0)], _v[edges_quad(num, 1)]);
   }
   virtual int getNumEdgesRep(){ return 4; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -664,6 +853,11 @@ class MQuadrangle : public MElement {
     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], _v[3]); }
   virtual int getNumFacesRep(){ return 2; }
@@ -675,6 +869,11 @@ class MQuadrangle : public MElement {
     _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(v);
+  }
   virtual int getTypeForMSH(){ return MSH_QUA_4; }
   virtual int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
   virtual int getTypeForVTK(){ return 8; }
@@ -711,8 +910,35 @@ class MQuadrangle : public MElement {
     return true;
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ private:
+  int edges_quad(const int edge, const int vert) const
+  {
+    static const int e[4][2] = {
+      {0, 1},
+      {1, 2},
+      {2, 3},
+      {3, 0}
+    };
+    return e[edge][vert];
+  }
 };
 
+/*
+ * MQuadrangle8
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3-----6----2
+ *   |          |
+ *   |          |
+ *   7          5
+ *   |          |
+ *   |          |
+ *   0-----4----1 --> u
+ *
+ */
 class MQuadrangle8 : public MQuadrangle {
  protected:
   MVertex *_vs[4];
@@ -756,6 +982,12 @@ class MQuadrangle8 : public MQuadrangle {
     };
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MQuadrangle::_getEdgeVertices(num, v);
+    v[2] = _vs[num];
+  }
   virtual int getNumFacesRep(){ return 6; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -765,6 +997,15 @@ class MQuadrangle8 : public MQuadrangle {
     _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(8);
+    MQuadrangle::_getFaceVertices(v);
+    v[4] = _vs[0];
+    v[5] = _vs[1];
+    v[6] = _vs[2];
+    v[7] = _vs[3];
+  }
   virtual int getTypeForMSH(){ return MSH_QUA_8; }
   virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
   virtual int getTypeForVTK(){ return 23; }
@@ -778,6 +1019,22 @@ class MQuadrangle8 : public MQuadrangle {
   }
 };
 
+/*
+ * MQuadrangle9
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3-----6----2
+ *   |          |
+ *   |          |
+ *   7     8    5
+ *   |          |
+ *   |          |
+ *   0-----4----1 --> u
+ *
+ */
 class MQuadrangle9 : public MQuadrangle {
  protected:
   MVertex *_vs[5];
@@ -812,6 +1069,12 @@ class MQuadrangle9 : public MQuadrangle {
     };
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MQuadrangle::getEdgeVertices(num, v);
+    v[2] = _vs[num];
+  }
   virtual int getNumFacesRep(){ return 8; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -822,6 +1085,16 @@ class MQuadrangle9 : public MQuadrangle {
     _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(9);
+    MQuadrangle::_getFaceVertices(v);
+    v[4] = _vs[0];
+    v[5] = _vs[1];
+    v[6] = _vs[2];
+    v[7] = _vs[3];
+    v[8] = _vs[4];
+  }
   virtual int getTypeForMSH(){ return MSH_QUA_9; }
   virtual const char *getStringForPOS(){ return "SQ2"; }
   virtual void revert() 
@@ -833,9 +1106,43 @@ class MQuadrangle9 : public MQuadrangle {
   }
 };
 
+/*
+ * MTetrahedron
+ *
+ *                      v
+ *                    .
+ *                  ,/
+ *                 /
+ *              2
+ *            ,/|`\
+ *          ,/  |  `\
+ *        ,/    '.   `\
+ *      ,/       |     `\
+ *    ,/         |       `\
+ *   0-----------'+--------1 --> u
+ *    `\.         |      ,/
+ *       `\.      |    ,/
+ *          `\.   '. ,/
+ *             `\. |/
+ *                `3
+ *                   `\.
+ *                      ` w
+ *
+ */
 class MTetrahedron : public MElement {
  protected:
   MVertex *_v[4];
+  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[edges_tetra(num, 0)];
+    v[1] = _v[edges_tetra(num, 1)];
+  }
+  void _getFaceVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[faces_tetra(num, 0)];
+    v[1] = _v[faces_tetra(num, 1)];
+    v[2] = _v[faces_tetra(num, 2)];
+  }
  public :
   MTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0) 
     : MElement(num, part)
@@ -859,15 +1166,7 @@ class MTetrahedron : public MElement {
   virtual int getNumEdges(){ return 6; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_tetra[6][2] = {
-      {0, 1},
-      {1, 2},
-      {2, 0},
-      {3, 0},
-      {3, 2},
-      {3, 1}
-    };
-    return MEdge(_v[edges_tetra[num][0]], _v[edges_tetra[num][1]]);
+    return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]);
   }
   virtual int getNumEdgesRep(){ return 6; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -876,18 +1175,17 @@ class MTetrahedron : public MElement {
     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 4; }
   virtual MFace getFace(int num)
   {
-    static const int faces_tetra[4][3] = {
-      {0, 2, 1},
-      {0, 1, 3},
-      {0, 3, 2},
-      {3, 1, 2}
-    };
-    return MFace(_v[faces_tetra[num][0]],
-                 _v[faces_tetra[num][1]],
-                 _v[faces_tetra[num][2]]);
+    return MFace(_v[faces_tetra(num, 0)],
+		 _v[faces_tetra(num, 1)],
+		 _v[faces_tetra(num, 2)]);
   }
   virtual int getNumFacesRep(){ return 4; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -895,6 +1193,11 @@ class MTetrahedron : public MElement {
     MFace f(getFace(num));
     _getFaceRep(f.getVertex(0), f.getVertex(1), f.getVertex(2), x, y, z, n);
   }
+  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    _getFaceVertices(num, v);
+  }
   virtual int getTypeForMSH(){ return MSH_TET_4; }
   virtual int getTypeForUNV(){ return 111; } // solid linear tetrahedron
   virtual int getTypeForVTK(){ return 10; }
@@ -956,8 +1259,54 @@ class MTetrahedron : public MElement {
     return true;
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ private:
+  int edges_tetra(const int edge, const int vert) const
+  {
+    static const int e[6][2] = {
+      {0, 1},
+      {1, 2},
+      {2, 0},
+      {3, 0},
+      {3, 2},
+      {3, 1}
+    };
+    return e[edge][vert];
+  }
+  int faces_tetra(const int face, const int vert) const
+  {
+    static const int f[4][3] = {
+      {0, 2, 1},
+      {0, 1, 3},
+      {0, 3, 2},
+      {3, 1, 2}
+    };
+    return f[face][vert];
+  }
 };
 
+/*
+ * MTetrahedron10
+ *
+ *                      v
+ *                    .
+ *                  ,/
+ *                 /
+ *              2 
+ *            ,/|`\
+ *          ,/  |  `\
+ *        ,6    '.   `5
+ *      ,/       8     `\
+ *    ,/         |       `\
+ *   0--------4--'+--------1 --> u
+ *    `\.         |      ,/
+ *       `\.      |    ,9
+ *          `7.   '. ,/
+ *             `\. |/
+ *                `3
+ *                   `\.
+ *                      ` w
+ *
+ */
 class MTetrahedron10 : public MTetrahedron {
  protected:
   MVertex *_vs[6];
@@ -1010,6 +1359,12 @@ class MTetrahedron10 : public MTetrahedron {
     static const int f[12] = {0, 0, 0, 1, 2, 3};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MTetrahedron::_getEdgeVertices(num, v);
+    v[2] = _vs[num];
+  }
   virtual int getNumFacesRep(){ return 16; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1022,6 +1377,20 @@ class MTetrahedron10 : public MTetrahedron {
     _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(6);
+    MTetrahedron::_getFaceVertices(num, v);
+    static const int f[4][3] = {
+      {4, 5, 6},
+      {4, 7, 9},
+      {6, 7, 8},
+      {5, 8, 9}
+    };
+    v[3] = _vs[f[num][0]];
+    v[4] = _vs[f[num][1]];
+    v[5] = _vs[f[num][2]];
+  }
   virtual int getTypeForMSH(){ return MSH_TET_10; }
   virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
   virtual int getTypeForVTK(){ return 24; }
@@ -1045,7 +1414,29 @@ class MTetrahedron10 : public MTetrahedron {
   
 };
 
-
+/*
+ * MTetrahedronN
+ *
+ *                      v           E = order - 1
+ *                    .             C = 4 + 6*E
+ *                  ,/              F = ((order - 1)*(order - 2))/2
+ *                 /                N = total number of vertices
+ *              2
+ *            ,/|`\                 Interior vertex numbers
+ *          ,/  |  `\                 for edge 0 <= i <= 5: 4+i*E to 3+(i+1)*E
+ *        ,/    '.   `\               for face 0 <= j <= 3: C+j*F to C-1+(j+1)*F
+ *      ,/       |     `\             in volume           : C+4*F to N-1
+ *    ,/         |       `\
+ *   0-----------'+--------1 --> u
+ *    `\.         |      ,/
+ *       `\.      |    ,/
+ *          `\.   '. ,/
+ *             `\. |/
+ *                `3
+ *                   `\.
+ *                      ` w
+ *
+ */
 class MTetrahedronN : public MTetrahedron {
  protected:
   std::vector<MVertex *> _vs;
@@ -1067,16 +1458,27 @@ class MTetrahedronN : public MTetrahedron {
   virtual int getPolynomialOrder(){ return _order; }
   virtual int getNumVertices(){ return 4 + _vs.size(); }
   virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
-  virtual int getNumFaceVertices()
+  virtual int getNumEdgeVertices() const { return _order - 1; }
+  virtual int getNumFaceVertices() const
+  {
+    return ((_order - 1)*(_order - 2))/2;
+  }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(2 + getNumEdgeVertices());
+    MTetrahedron::_getEdgeVertices(num, v);
+    int j = 2;
+    const int ie = (num+1)*getNumEdgeVertices();
+    for(int i = num*getNumEdgeVertices(); i != ie; ++i) v[j++] = _vs[i];
+  }
+  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
   {
-    switch(getTypeForMSH()){
-    case MSH_TET_20 : return 1;
-    case MSH_TET_35 : return 3;
-    case MSH_TET_56 : return 6;
-    default : return 0;
-    }    
+    v.resize(3 + getNumFaceVertices());
+    MTetrahedron::_getFaceVertices(num, v);
+    int j = 3;
+    const int ie = (num+1)*getNumFaceVertices();
+    for(int i = num*getNumFaceVertices(); i != ie; ++i) v[j++] = _vs[i];
   }
-  virtual int getNumEdgeVertices(){ return _order - 1; }
   virtual int getTypeForMSH()
   {
     // (p+1)*(p+2)*(p+3)/6
@@ -1103,10 +1505,44 @@ class MTetrahedronN : public MTetrahedron {
   }
 };
 
-
+/*
+ * MHexahedron
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3----------2
+ *   |\         |\
+ *   | \        | \
+ *   |  \       |  \
+ *   |   7------+---6
+ *   |   |      |   |
+ *   0---+------1   | --> u
+ *    \  |       \  |
+ *     \ |        \ |
+ *      \|         \|
+ *       4----------5
+ *        \
+ *         \
+ *          w
+ *
+ */
 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) 
@@ -1132,21 +1568,7 @@ class MHexahedron : public MElement {
   virtual int getNumEdges(){ return 12; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_hexa[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 MEdge(_v[edges_hexa[num][0]], _v[edges_hexa[num][1]]);
+    return MEdge(_v[edges_hexa(num, 0)], _v[edges_hexa(num, 1)]);
   }
   virtual int getNumEdgesRep(){ return 12; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1155,21 +1577,18 @@ class MHexahedron : public MElement {
     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)
   {
-    static const int faces_hexa[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 MFace(_v[faces_hexa[num][0]],
-                 _v[faces_hexa[num][1]],
-                 _v[faces_hexa[num][2]],
-                 _v[faces_hexa[num][3]]);
+    return MFace(_v[faces_hexa(num, 0)],
+		 _v[faces_hexa(num, 1)],
+		 _v[faces_hexa(num, 2)],
+		 _v[faces_hexa(num, 3)]);
   }
   virtual int getNumFacesRep(){ return 12; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1185,6 +1604,11 @@ class MHexahedron : public MElement {
     _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 getTypeForMSH(){ return MSH_HEX_8; }
   virtual int getTypeForUNV(){ return 115; } // solid linear brick
   virtual int getTypeForVTK(){ return 12; }
@@ -1249,8 +1673,62 @@ class MHexahedron : public MElement {
     return true;
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ private:
+  int edges_hexa(const int edge, const int vert) const
+  {
+    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];
+  }
+  int faces_hexa(const int face, const int vert) const
+  {
+    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];
+  }
 };
 
+/*
+ * MHexahedron20
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3----13----2
+ *   |\         |\
+ *   |15        | 14
+ *   9  \       11 \
+ *   |   7----19+---6
+ *   |   |      |   |
+ *   0---+-8----1   | --> u
+ *    \ 17       \  18
+ *    10 |        12|
+ *      \|         \|
+ *       4----16----5
+ *        \
+ *         \
+ *          w
+ *
+ */
 class MHexahedron20 : public MHexahedron {
  protected:
   MVertex *_vs[12];
@@ -1316,6 +1794,12 @@ class MHexahedron20 : public MHexahedron {
     static const int f[12] = {0, 0, 1, 0, 1, 0, 3, 2, 1, 2, 3, 4};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  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(){ return 36; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1330,6 +1814,23 @@ class MHexahedron20 : public MHexahedron {
     _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(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(){ return MSH_HEX_20; }
   virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
   virtual int getTypeForVTK(){ return 25; }
@@ -1349,6 +1850,29 @@ class MHexahedron20 : public MHexahedron {
   }
 };
 
+/*
+ * MHexahedron27
+ *
+ *   v
+ *
+ *   ^
+ *   |
+ *   3----13----2
+ *   |\         |\
+ *   |15    24  | 14
+ *   9  \ 20    11 \
+ *   |   7----19+---6
+ *   |22 |  26  | 23|
+ *   0---+-8----1   | --> u
+ *    \ 17    25 \  18
+ *    10 |  21    12|
+ *      \|         \|
+ *       4----16----5
+ *        \
+ *         \
+ *          w
+ *
+ */
 class MHexahedron27 : public MHexahedron {
  protected:
   MVertex *_vs[19];
@@ -1400,6 +1924,12 @@ class MHexahedron27 : public MHexahedron {
     static const int f[12] = {0, 0, 1, 0, 1, 0, 3, 2, 1, 2, 3, 4};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  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(){ return 48; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1420,6 +1950,24 @@ class MHexahedron27 : public MHexahedron {
     _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(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(){ return MSH_HEX_27; }
   virtual const char *getStringForPOS(){ return "SH2"; }
   virtual void revert()
@@ -1441,9 +1989,44 @@ class MHexahedron27 : public MHexahedron {
   }
 };
 
+/*
+ * MPrism
+ *
+ *               w
+ *
+ *               ^
+ *               |
+ *               3
+ *             ,/|`\
+ *           ,/  |  `\
+ *         ,/    |    `\
+ *        4------+------5
+ *        |      |      |
+ *        |      0      |
+ *        |    ,/ `\    |
+ *        |  ,/     `\  |
+ *        |,/         `\|
+ *        1-------------2
+ *     ,/                 `\
+ *   ,/                     `\
+ *  u                          v
+ *
+ */
 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) 
@@ -1469,18 +2052,7 @@ class MPrism : public MElement {
   virtual int getNumEdges(){ return 9; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_prism[9][2] = {
-      {0, 1},
-      {0, 2},
-      {0, 3},
-      {1, 2},
-      {1, 4},
-      {2, 5},
-      {3, 4},
-      {3, 5},
-      {4, 5}
-    };
-    return MEdge(_v[edges_prism[num][0]], _v[edges_prism[num][1]]);
+    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)
@@ -1489,27 +2061,23 @@ class MPrism : public MElement {
     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 MFace getFace(int num)
   {
-    static const int trifaces_prism[2][3] = {
-      {0, 2, 1},
-      {3, 4, 5}
-    };
-    static const int quadfaces_prism[3][4] = {
-      {0, 1, 4, 3},
-      {0, 3, 5, 2},
-      {1, 2, 5, 4}
-    };
     if(num < 2)
-      return MFace(_v[trifaces_prism[num][0]],
-                   _v[trifaces_prism[num][1]],
-                   _v[trifaces_prism[num][2]]);
+      return MFace(_v[faces_prism(num, 0)],
+                   _v[faces_prism(num, 1)],
+                   _v[faces_prism(num, 2)]);
     else
-      return MFace(_v[quadfaces_prism[num - 2][0]],
-                   _v[quadfaces_prism[num - 2][1]],
-                   _v[quadfaces_prism[num - 2][2]],
-                   _v[quadfaces_prism[num - 2][3]]);
+      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)
@@ -1524,6 +2092,11 @@ class MPrism : public MElement {
     _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 getTypeForMSH(){ return MSH_PRI_6; }
   virtual int getTypeForUNV(){ return 112; } // solid linear wedge
   virtual int getTypeForVTK(){ return 13; }
@@ -1579,8 +2152,58 @@ class MPrism : public MElement {
       return false;
     return true;
   }
+ private:
+  int edges_prism(const int edge, const int vert) const
+  {
+    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];
+  }
+  int faces_prism(const int face, const int vert) const
+  {
+    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
+ *
+ *               w
+ *
+ *               ^
+ *               |
+ *               3
+ *             ,/|`\
+ *           12  |  13
+ *         ,/    8    `\
+ *        4------14-----5
+ *        |      |      |
+ *        |      0      |
+ *       10    ,/ `\    11
+ *        |  ,6     `7  |
+ *        |,/         `\|
+ *        1------9------2
+ *     ,/                 `\
+ *   ,/                     `\
+ *  u                          v
+ *
+ */
 class MPrism15 : public MPrism {
  protected:
   MVertex *_vs[9];
@@ -1638,6 +2261,12 @@ class MPrism15 : public MPrism {
     static const int f[9] = {0, 1, 2, 0, 2, 3, 1, 1, 1};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  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(){ return 26; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1651,6 +2280,23 @@ class MPrism15 : public MPrism {
     _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) ? 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(){ return MSH_PRI_15; }
   virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
   virtual const char *getStringForBDF(){ return "CPENTA"; }
@@ -1665,6 +2311,29 @@ class MPrism15 : public MPrism {
   }
 };
 
+/*
+ * MPrism18
+ *
+ *               w
+ *
+ *               ^
+ *               |
+ *               3
+ *             ,/|`\
+ *           12  |  13
+ *         ,/    8    `\
+ *        4------14-----5
+ *        |  15  |  16  |
+ *        |,/    0    `\|
+ *       10------17-----11
+ *        |  ,6     `7  |
+ *        |,/         `\|
+ *        1------9------2
+ *     ,/                 `\
+ *   ,/                     `\
+ *  u                          v
+ *
+ */
 class MPrism18 : public MPrism {
  protected:
   MVertex *_vs[12];
@@ -1709,6 +2378,12 @@ class MPrism18 : public MPrism {
     static const int f[9] = {0, 1, 2, 0, 2, 3, 1, 1, 1};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  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(){ return 32; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1725,6 +2400,26 @@ class MPrism18 : public MPrism {
     _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) ? 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(){ return MSH_PRI_18; }
   virtual const char *getStringForPOS(){ return "SI2"; }
   virtual void revert()
@@ -1741,9 +2436,50 @@ class MPrism18 : public MPrism {
   }
 };
 
+/*
+ * MPyramid
+ *
+ *                 4
+ *               ,/|\
+ *             ,/ .'|\
+ *           ,/   | | \
+ *  w      ,/    .' | `.
+ *       ,/      |  '. \
+ *  ^  ,/       .'   |  \
+ *  |,/         |    |   \
+ *  0----------.'----3    \ --> v
+ *   `\        |      `\  `.
+ *     `\     .'        `\ \
+ *       `\   |           `\\
+ *         `\.'             `\
+ *           1----------------2
+ *             `\
+ *               `\
+ *                  u
+ *
+ */
 class MPyramid : public MElement {
  protected:
   MVertex *_v[5];
+  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v[0] = _v[edges_pyramid(num, 0)];
+    v[1] = _v[edges_pyramid(num, 1)];
+  }
+  void _getFaceVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    if(num < 4) {
+      v[0] = _v[faces_pyramid(num, 0)];
+      v[1] = _v[faces_pyramid(num, 1)];
+      v[2] = _v[faces_pyramid(num, 2)];
+    }
+    else {
+      v[0] = _v[0];
+      v[1] = _v[3];
+      v[2] = _v[2];
+      v[3] = _v[1];
+    }
+  }
  public :
   MPyramid(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, 
            int num=0, int part=0) 
@@ -1768,17 +2504,7 @@ class MPyramid : public MElement {
   virtual int getNumEdges(){ return 8; }
   virtual MEdge getEdge(int num)
   {
-    static const int edges_pyramid[8][2] = {
-      {0, 1},
-      {0, 3},
-      {0, 4},
-      {1, 2},
-      {1, 4},
-      {2, 3},
-      {2, 4},
-      {3, 4}
-    };
-    return MEdge(_v[edges_pyramid[num][0]], _v[edges_pyramid[num][1]]);
+    return MEdge(_v[edges_pyramid(num, 0)], _v[edges_pyramid(num, 1)]);
   }
   virtual int getNumEdgesRep(){ return 8; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1787,19 +2513,18 @@ class MPyramid : public MElement {
     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 MFace getFace(int num)
   {
-    static const int faces_pyramid[4][3] = {
-      {0, 1, 4},
-      {3, 0, 4},
-      {1, 2, 4},
-      {2, 3, 4}
-    };
     if(num < 4)
-      return MFace(_v[faces_pyramid[num][0]],
-                   _v[faces_pyramid[num][1]],
-                   _v[faces_pyramid[num][2]]);
+      return MFace(_v[faces_pyramid(num, 0)],
+                   _v[faces_pyramid(num, 1)],
+                   _v[faces_pyramid(num, 2)]);
     else
       return MFace(_v[0], _v[3], _v[2], _v[1]);
   }
@@ -1816,6 +2541,11 @@ class MPyramid : public MElement {
     _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 < 4) ? 3 : 4);
+    _getFaceVertices(num, v);
+  }
   virtual int getTypeForMSH(){ return MSH_PYR_5; }
   virtual int getTypeForVTK(){ return 14; }
   virtual const char *getStringForPOS(){ return "SY"; }
@@ -1874,8 +2604,55 @@ class MPyramid : public MElement {
       return false;
     return true;
   }
+ private:
+  int edges_pyramid(const int edge, const int vert) const
+  {
+    static const int e[8][2] = {
+      {0, 1},
+      {0, 3},
+      {0, 4},
+      {1, 2},
+      {1, 4},
+      {2, 3},
+      {2, 4},
+      {3, 4}
+    };
+    return e[edge][vert];
+  }
+  int faces_pyramid(const int face, const int vert) const
+  {
+    static const int f[4][3] = {
+      {0, 1, 4},
+      {3, 0, 4},
+      {1, 2, 4},
+      {2, 3, 4}
+    };
+    return f[face][vert];
+  }
 };
 
+/*
+ * MPyramid13
+ *
+ *                 4
+ *               ,/|\
+ *             ,/ .'|\
+ *           ,/   | | \
+ *  w      ,7    .' | `.
+ *       ,/      |  12 \
+ *  ^  ,/       .'   |  \
+ *  |,/         9    |   11
+ *  0--------6-.'----3    \ --> v
+ *   `\        |      `\  `.
+ *     `5     .'        10 \
+ *       `\   |           `\\
+ *         `\.'             `\
+ *           1--------8-------2
+ *             `\
+ *               `\
+ *                  u
+ *
+ */
 class MPyramid13 : public MPyramid {
  protected:
   MVertex *_vs[8];
@@ -1921,6 +2698,12 @@ class MPyramid13 : public MPyramid {
     static const int f[8] = {0, 1, 1, 2, 0, 3, 2, 3};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MPyramid::_getEdgeVertices(num, v);
+    v[2] = _vs[num];
+  }
   virtual int getNumFacesRep(){ return 22; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -1934,6 +2717,28 @@ class MPyramid13 : public MPyramid {
     _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 < 4) ? 6 : 8);
+    MPyramid::_getFaceVertices(num, v);
+    static const int f[4][3] = {
+      {0, 4, 2},
+      {1, 2, 7},
+      {3, 6, 4},
+      {5, 7, 6}
+    };
+    if(num < 4) {
+      v[3] = _vs[f[num][0]];
+      v[4] = _vs[f[num][1]];
+      v[5] = _vs[f[num][2]];
+    }
+    else {
+      v[4] = _vs[1];
+      v[5] = _vs[5];
+      v[6] = _vs[3];
+      v[7] = _vs[0];
+    }
+  }
   virtual int getTypeForMSH(){ return MSH_PYR_13; }
   virtual void revert()
   {
@@ -1945,6 +2750,28 @@ class MPyramid13 : public MPyramid {
   }
 };
 
+/*
+ * MPyramid14
+ *
+ *                 4
+ *               ,/|\
+ *             ,/ .'|\
+ *           ,/   | | \
+ *  w      ,7    .' | `.
+ *       ,/      |  12 \
+ *  ^  ,/       .'   |  \
+ *  |,/         9    |   11
+ *  0--------6-.'----3    \ --> v
+ *   `\       `|      `\  `.
+ *     `5-----.'13------10 \
+ *       `\   |   `\      `\\
+ *         `\.'     `\      `\
+ *           1--------8-------2
+ *             `\
+ *               `\
+ *                  u
+ *
+ */
 class MPyramid14 : public MPyramid {
  protected:
   MVertex *_vs[9];
@@ -1987,6 +2814,12 @@ class MPyramid14 : public MPyramid {
     static const int f[8] = {0, 1, 1, 2, 0, 3, 2, 3};
     _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
   }
+  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+  {
+    v.resize(3);
+    MPyramid::_getEdgeVertices(num, v);
+    v[2] = _vs[num];
+  }
   virtual int getNumFacesRep(){ return 24; }
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
   { 
@@ -2001,6 +2834,29 @@ class MPyramid14 : public MPyramid {
     _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 < 4) ? 6 : 9);
+    MPyramid::_getFaceVertices(num, v);
+    static const int f[4][3] = {
+      {0, 4, 2},
+      {1, 2, 7},
+      {3, 6, 4},
+      {5, 7, 6}
+    };
+    if(num < 4) {
+      v[3] = _vs[f[num][0]];
+      v[4] = _vs[f[num][1]];
+      v[5] = _vs[f[num][2]];
+    }
+    else {
+      v[4] = _vs[1];
+      v[5] = _vs[5];
+      v[6] = _vs[3];
+      v[7] = _vs[0];
+      v[8] = _vs[8];
+    }
+  }
   virtual int getTypeForMSH(){ return MSH_PYR_14; }
   virtual const char *getStringForPOS(){ return "SY2"; }
   virtual void revert()
-- 
GitLab