Improve P2 prisms visualisation

Geo/MPrism.cpp

@@ -6,6 +6,7 @@
 #include "MPrism.h"
 #include "Numeric.h"
 #include "BasisFactory.h"
+#include "Context.h"
 
 int MPrism::getVolumeSign()
 { 
@@ -159,3 +160,295 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
   }
   Msg::Error("Could not get face information for prism %d", getNum());
 }
+
+static void _myGetEdgeRep(MPrism *pri, int num, double *x, double *y, double *z,
+                          SVector3 *n, int numSubEdges) {
+
+  //const int numSubEdges = CTX::instance()->mesh.numSubEdges;
+  static double pp[6][3] = {
+    {0,0,-1},{1,0,-1},{0,1,-1},
+    {0,0,1},{1,0,1},{0,1,1} };
+  static int ed [9][2] = {
+    {0,1},{0,2},{0,3},{1,2},{1,4},{2,5},
+    {3,4},{3,5},{4,5}
+  };
+  int iEdge = num / numSubEdges;
+  int iSubEdge = num % numSubEdges;
+
+  int iVertex1 = ed [iEdge][0];
+  int iVertex2 = ed [iEdge][1];
+  double t1 = (double) iSubEdge / (double) numSubEdges;
+  double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+  double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+  double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+
+  double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+  double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+  double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+  double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+
+  SPoint3 pnt1, pnt2;
+  pri->pnt(u1, v1, w1, pnt1);
+  pri->pnt(u2, v2, w2, pnt2);
+  x[0] = pnt1.x(); x[1] = pnt2.x();
+  y[0] = pnt1.y(); y[1] = pnt2.y();
+  z[0] = pnt1.z(); z[1] = pnt2.z();
+
+  // not great, but better than nothing
+  //static const int f[6] = {0, 0, 0, 1, 2, 3};
+  n[0] = n[1] = 1 ;
+}
+
+
+void MPrism15::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) {
+  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+void MPrism18::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) {
+  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+int MPrism15::getNumEdgesRep() {
+  return 9 * CTX::instance()->mesh.numSubEdges;
+}
+
+int MPrism18::getNumEdgesRep() {
+  return 9 * CTX::instance()->mesh.numSubEdges;
+}
+
+
+void _myGetFaceRep(MPrism *pri, int num, double *x, double *y, double *z,
+                          SVector3 *n, int numSubEdges)
+{
+  static double pp[6][3] = {
+    {0,0,-1},{1,0,-1},{0,1,-1},
+    {0,0,1},{1,0,1},{0,1,1} };
+
+  int iFace    = num / (numSubEdges * numSubEdges);
+  int iSubFace = num % (numSubEdges * numSubEdges);
+
+  if (iFace > 1) {
+    iFace = num / (2*numSubEdges * numSubEdges) + 1;
+    iSubFace = num % (2*numSubEdges * numSubEdges);
+  }
+
+  int iVertex1 = pri->faces_prism(iFace,0);
+  int iVertex2 = pri->faces_prism(iFace,1);
+  int iVertex3 = pri->faces_prism(iFace,2);
+  int iVertex4 = pri->faces_prism(iFace,3);
+
+  SPoint3 pnt1, pnt2, pnt3;
+ // double J1[3][3], J2[3][3], J3[3][3];
+
+  /*
+    0
+    0 1
+    0 1 2
+    0 1 2 3
+    0 1 2 3 4
+    0 1 2 3 4 5
+  */
+
+  // on each layer, we have (numSubEdges) * 2 triangles
+  // ix and iy are the coordinates of the sub-quadrangle
+
+  if (iFace > 1) {
+    int io = iSubFace%2;
+    int ix = (iSubFace/2)/numSubEdges;
+    int iy = (iSubFace/2)%numSubEdges;
+
+    const double d = 2. / numSubEdges;
+    double ox = -1. + d*ix;
+    double oy = -1. + d*iy;
+
+    if (io == 0){
+      double U1 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V1 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W1 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      ox += d;
+
+      double U2 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V2 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W2 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      oy += d;
+
+      double U3 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V3 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W3 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      pri->pnt(U1, V1, W1, pnt1);
+      pri->pnt(U2, V2, W2, pnt2);
+      pri->pnt(U3, V3, W3, pnt3);
+    }
+    else{
+      double U1 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V1 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W1 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      ox += d;
+      oy += d;
+
+      double U2 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V2 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W2 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      ox -= d;
+
+      double U3 =
+        pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+      double V3 =
+        pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+      double W3 =
+        pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+        pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+        pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+        pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+      pri->pnt(U1, V1, W1, pnt1);
+      pri->pnt(U2, V2, W2, pnt2);
+      pri->pnt(U3, V3, W3, pnt3);
+    }
+
+  } else {
+    int ix = 0, iy = 0;
+    int nbt = 0;
+    for (int i = 0; i < numSubEdges; i++){
+      int nbl = (numSubEdges - i - 1) * 2 + 1;
+      nbt += nbl;
+      if (nbt > iSubFace){
+        iy = i;
+        ix = nbl - (nbt - iSubFace);
+        break;
+      }
+    }
+
+    const double d = 1. / numSubEdges;
+
+    double u1, v1, u2, v2, u3, v3;
+    if (ix % 2 == 0){
+      u1 = ix / 2 * d; v1= iy*d;
+      u2 = (ix / 2 + 1) * d ; v2 =  iy * d;
+      u3 = ix / 2 * d ; v3 =  (iy+1) * d;
+    }
+    else{
+      u1 = (ix / 2 + 1) * d; v1= iy * d;
+      u2 = (ix / 2 + 1) * d; v2= (iy + 1) * d;
+      u3 = ix / 2 * d ; v3 =  (iy + 1) * d;
+    }
+
+    double U1 = pp[iVertex1][0] * (1.-u1-v1) + pp[iVertex2][0] * u1 + pp[iVertex3][0] * v1;
+    double U2 = pp[iVertex1][0] * (1.-u2-v2) + pp[iVertex2][0] * u2 + pp[iVertex3][0] * v2;
+    double U3 = pp[iVertex1][0] * (1.-u3-v3) + pp[iVertex2][0] * u3 + pp[iVertex3][0] * v3;
+
+    double V1 = pp[iVertex1][1] * (1.-u1-v1) + pp[iVertex2][1] * u1 + pp[iVertex3][1] * v1;
+    double V2 = pp[iVertex1][1] * (1.-u2-v2) + pp[iVertex2][1] * u2 + pp[iVertex3][1] * v2;
+    double V3 = pp[iVertex1][1] * (1.-u3-v3) + pp[iVertex2][1] * u3 + pp[iVertex3][1] * v3;
+
+    double W1 = pp[iVertex1][2] * (1.-u1-v1) + pp[iVertex2][2] * u1 + pp[iVertex3][2] * v1;
+    double W2 = pp[iVertex1][2] * (1.-u2-v2) + pp[iVertex2][2] * u2 + pp[iVertex3][2] * v2;
+    double W3 = pp[iVertex1][2] * (1.-u3-v3) + pp[iVertex2][2] * u3 + pp[iVertex3][2] * v3;
+
+    pri->pnt(U1, V1, W1, pnt1);
+    pri->pnt(U2, V2, W2, pnt2);
+    pri->pnt(U3, V3, W3, pnt3);
+  }
+
+
+
+  x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
+  y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
+  z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
+
+  SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
+  SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
+  n[0] = crossprod(d1, d2);
+  n[0].normalize();
+  n[1] = n[0];
+  n[2] = n[0];
+}
+
+void MPrism15::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+void MPrism18::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+int MPrism15::getNumFacesRep() {
+  return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
+
+int MPrism18::getNumFacesRep() {
+  return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
\ No newline at end of file

Geo/MPrism.h

@@ -241,42 +241,16 @@ class MPrism15 : public MPrism {
   }
   virtual MVertex *getVertexINP(int num){ return getVertexBDF(num); }
   virtual int getNumEdgeVertices() const { return 9; }
-  virtual int getNumEdgesRep(){ return 18; }
-  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
-  {
-    static const int e[18][2] = {
-      {0, 6}, {6, 1},
-      {0, 7}, {7, 2},
-      {0, 8}, {8, 3},
-      {1, 9}, {9, 2},
-      {1, 10}, {10, 4},
-      {2, 11}, {11, 5},
-      {3, 12}, {12, 4},
-      {3, 13}, {13, 5},
-      {4, 14}, {14, 5}
-    };
-    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 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(){ return 26; }
-  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
-  {
-    static const int f[26][3] = {
-      {0, 7, 6}, {2, 9, 7}, {1, 6, 9}, {6, 7, 9},
-      {3, 12, 13}, {4, 14, 12}, {5, 13, 14}, {12, 14, 13},
-      {0, 6, 8}, {1, 10, 6}, {4, 12, 10}, {3, 8, 12}, {6, 10, 12}, {6, 12, 8},
-      {0, 8, 7}, {3, 13, 8}, {5, 11, 13}, {2, 7, 11}, {7, 8, 13}, {7, 13, 11},
-      {1, 9, 10}, {2, 11, 9}, {5, 14, 11}, {4, 10, 14}, {9, 11, 14}, {9, 14, 10}
-    };
-    _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
-                x, y, z, n);
-  }
+  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);
@@ -361,45 +335,16 @@ class MPrism18 : public MPrism {
   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(){ return 18; }
-  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
-  {
-    static const int e[18][2] = {
-      {0, 6}, {6, 1},
-      {0, 7}, {7, 2},
-      {0, 8}, {8, 3},
-      {1, 9}, {9, 2},
-      {1, 10}, {10, 4},
-      {2, 11}, {11, 5},
-      {3, 12}, {12, 4},
-      {3, 13}, {13, 5},
-      {4, 14}, {14, 5}
-    };
-    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 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(){ return 32; }
-  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
-  {
-    static const int f[32][3] = {
-      {0, 7, 6}, {2, 9, 7}, {1, 6, 9}, {6, 7, 9},
-      {3, 12, 13}, {4, 14, 12}, {5, 13, 14}, {12, 14, 13},
-      {0, 6, 15}, {0, 15, 8}, {1, 10, 15}, {1, 15, 6},
-      {4, 12, 15}, {4, 15, 10}, {3, 8, 15}, {3, 15, 12},
-      {0, 8, 16}, {0, 16, 7}, {3, 13, 16}, {3, 16, 8},
-      {5, 11, 16}, {5, 16, 13}, {2, 7, 16}, {2, 16, 11},
-      {1, 9, 17}, {1, 17, 10}, {2, 11, 17}, {2, 17, 9},
-      {5, 14, 17}, {5, 17, 11}, {4, 10, 17}, {4, 17, 14}
-    };
-    _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
-                x, y, z, n);
-  }
+  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);