fix bug reversed tetrahedron / prism / pyramid (fix #293)

Geo/MPrism.cpp

@@ -7,11 +7,14 @@
 #include "Numeric.h"
 #include "BasisFactory.h"
 #include "Context.h"
+#include "pointsGenerators.h"
 
 #if defined(HAVE_MESH)
 #include "qualityMeasures.h"
 #endif
 
+std::map<int, indicesReversed> MPrismN::_order2indicesReversed;
+
 void MPrism::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
                         SVector3 *n)
 {
@@ -546,3 +549,49 @@ double MPrism::gammaShapeMeasure()
   return 0.;
 #endif
 }
+
+void _getIndicesReversed(int order, indicesReversed &indices)
+{
+  fullMatrix<double> ref = gmshGenerateMonomialsPrism(order);
+
+  indices.resize(ref.size1());
+  for (int i = 0; i < ref.size1(); ++i) {
+    const double u = ref(i, 0);
+    const double v = ref(i, 1);
+    const double w = ref(i, 2);
+    for (int j = 0; j < ref.size1(); ++j) {
+      if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
+        indices[i] = j;
+        break;
+      }
+    }
+  }
+}
+
+void MPrismN::reverse()
+{
+  std::map<int, indicesReversed>::iterator it;
+  it = _order2indicesReversed.find(_order);
+  if (it == _order2indicesReversed.end()) {
+    indicesReversed indices;
+    _getIndicesReversed(_order, indices);
+    _order2indicesReversed[_order] = indices;
+    it = _order2indicesReversed.find(_order);
+  }
+
+  indicesReversed &indices = it->second;
+
+  // copy vertices
+  std::vector<MVertex*> oldv(6 + _vs.size());
+  std::copy(_v, _v+6, oldv.begin());
+  std::copy(_vs.begin(), _vs.end(), oldv.begin()+6);
+
+  // reverse
+  for (int i = 0; i < 6; ++i) {
+    _v[i] = oldv[indices[i]];
+  }
+  for (int i = 0; i < _vs.size(); ++i) {
+    _vs[i] = oldv[indices[6+i]];
+  }
+}
+

Geo/MPrism.h

@@ -397,7 +397,12 @@ class MPrism18 : public MPrism {
 /*
  * MPrismN
  */
+
+typedef std::vector<int> indicesReversed;
+
 class MPrismN : public MPrism {
+  static std::map<int, indicesReversed> _order2indicesReversed;
+
  protected:
   std::vector<MVertex *> _vs;
   const char _order;
@@ -504,10 +509,7 @@ class MPrismN : public MPrism {
     }
     return "";
   }
-  virtual void reverse()
-  {
-    Msg::Error("Reverse not implemented yet for MPrismN");
-  }
+  virtual void reverse();
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
     const fullMatrix<double> &p = getFunctionSpace()->points;

Geo/MPyramid.cpp

@@ -8,6 +8,9 @@
 #include "Context.h"
 #include "BergotBasis.h"
 #include "BasisFactory.h"
+#include "pointsGenerators.h"
+
+std::map<int, indicesReversed> MPyramidN::_order2indicesReversed;
 
 int MPyramid::getVolumeSign()
 {
@@ -331,3 +334,49 @@ void MPyramidN::getFaceRep(bool curved, int num,
     _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
   else MPyramid::getFaceRep(false, num, x, y, z, n);
 }
+
+void _getIndicesReversed(int order, indicesReversed &indices)
+{
+  fullMatrix<double> ref = gmshGenerateMonomialsPyramid(order);
+
+  indices.resize(ref.size1());
+  for (int i = 0; i < ref.size1(); ++i) {
+    const double u = ref(i, 0);
+    const double v = ref(i, 1);
+    const double w = ref(i, 2);
+    for (int j = 0; j < ref.size1(); ++j) {
+      if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
+        indices[i] = j;
+        break;
+      }
+    }
+  }
+}
+
+void MPyramidN::reverse()
+{
+  std::map<int, indicesReversed>::iterator it;
+  it = _order2indicesReversed.find(_order);
+  if (it == _order2indicesReversed.end()) {
+    indicesReversed indices;
+    _getIndicesReversed(_order, indices);
+    _order2indicesReversed[_order] = indices;
+    it = _order2indicesReversed.find(_order);
+  }
+
+  indicesReversed &indices = it->second;
+
+  // copy vertices
+  std::vector<MVertex*> oldv(5 + _vs.size());
+  std::copy(_v, _v+5, oldv.begin());
+  std::copy(_vs.begin(), _vs.end(), oldv.begin()+5);
+
+  // reverse
+  for (int i = 0; i < 5; ++i) {
+    _v[i] = oldv[indices[i]];
+  }
+  for (int i = 0; i < _vs.size(); ++i) {
+    _vs[i] = oldv[indices[5+i]];
+  }
+}
+

Geo/MPyramid.h

@@ -219,7 +219,10 @@ class MPyramid : public MElement {
 
 //------------------------------------------------------------------------------
 
+typedef std::vector<int> indicesReversed;
+
 class MPyramidN : public MPyramid {
+  static std::map<int, indicesReversed> _order2indicesReversed;
 
  protected:
   std::vector<MVertex*> _vs;
@@ -317,19 +320,7 @@ class MPyramidN : public MPyramid {
     Msg::Error("no tag matches a p%d pyramid with %d vertices", _order, 5+_vs.size());
     return 0;
   }
-  virtual void reverse()
-  {
-    /*
-    MVertex *tmp;
-    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
-    std::vector<MVertex*> inv(_vs.size());
-    std::vector<int> reverseIndices = _getReverseIndices(_order);
-    for (unsigned int i = 0; i< _vs.size(); i++)
-      inv[i] = _vs[reverseIndices[i + 4] - 4];
-    _vs = inv;
-    */
-    Msg::Error("Reverse not implemented yet for MPyramidN");
-  }
+  virtual void reverse();
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z,
                           SVector3 *n);
   virtual int getNumEdgesRep(bool curved);

Geo/MTetrahedron.cpp

@@ -8,6 +8,7 @@
 #include "Numeric.h"
 #include "Context.h"
 #include "BasisFactory.h"
+#include "pointsGenerators.h"
 
 #if defined(HAVE_MESH)
 #include "qualityMeasures.h"
@@ -17,6 +18,8 @@
 
 #define SQU(a)      ((a)*(a))
 
+std::map<int, indicesReversed> MTetrahedronN::_order2indicesReversed;
+
 void MTetrahedron::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
                               SVector3 *n)
 {
@@ -315,104 +318,47 @@ void MTetrahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &
   Msg::Error("Could not get face information for tetrahedron %d", getNum());
 }
 
-static std::vector<std::vector<int> > tetReverseIndices(20);
-
-const std::vector<int> &MTetrahedronN::_getReverseIndices (int order)
+void _getIndicesReversed(int order, indicesReversed &indices)
 {
-  if(order >= (int)tetReverseIndices.size())
-    tetReverseIndices.resize(order + 1);
-  std::vector<int> &r = tetReverseIndices[order];
-  if (r.size() != 0) return r;
-  //
-  // not the funniest code ever ... (guaranteed correct only up to order 5)
-  //
-  int nb = (order+1)*(order+2)*(order+3)/6;
-  r.resize(nb);
-  int p=0;
-  for (int layerOrder = order; layerOrder>=0; layerOrder-=4) {
-    //principal vertices
-    r[p+0] = p+0;
-    if (layerOrder ==0) break;
-    r[p+1] = p+2;
-    r[p+2] = p+1;
-    r[p+3] = p+3;
-    p+=4;
-    for (int i = 0; i<layerOrder-1; i++) {
-      //E2 reversed switches with E0
-      r[p+i] = p+3*(layerOrder-1)-(i+1);
-      r[p+3*(layerOrder-1)-(i+1)] = p+i;
-      //E1 is reversed
-      r[p+(layerOrder-1)+i] = p+2*(layerOrder-1)-(i+1);
-      //E3 is preserved
-      r[p+3*(layerOrder-1)+i] = p+3*(layerOrder-1)+i;
-      //E4 switches with E5
-      r[p+4*(layerOrder-1)+i] = p+5*(layerOrder-1)+i;
-      r[p+5*(layerOrder-1)+i] = p+4*(layerOrder-1)+i;
-    }
-    p+=6*(layerOrder-1);
-    //F0(=012) switches its nodes 1 and 2
-    for (int of = layerOrder-3; of >= 0; of -= 3) {
-      r[p] = p;
-      if (of == 0) {
-        p+=1;
+  fullMatrix<double> ref = gmshGenerateMonomialsTetrahedron(order);
+
+  indices.resize(ref.size1());
+  for (int i = 0; i < ref.size1(); ++i) {
+    const double u = ref(i, 0);
+    const double v = ref(i, 1);
+    const double w = ref(i, 2);
+    for (int j = 0; j < ref.size1(); ++j) {
+      if (u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
+        indices[i] = j;
         break;
       }
-      r[p+1] = p+2;
-      r[p+2] = p+1;
-      for (int i = 0; i < of-1; i++) {
-        //switch edges 0 and 2
-        r[p+3+i] = p+3+3*(of-1)-(i+1);
-        r[p+3+3*(of-1)-(i+1)] = p+3+i;
-        //reverse edge 1
-        r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
-      }
-      p += 3*of;
-    }
-    //F1 (=013) reversed switches with F2 (=032)
-    int nf = (layerOrder-2)*(layerOrder-1)/2;
-    for (int of = layerOrder-3; of >= 0; of -= 3) {
-      r[p] = p+nf;
-      r[p+nf] = p;
-      if (of == 0) {
-        p += 1;
-        break;
     }
-      r[p+1] = p+nf+2;
-      r[p+nf+2] = p+1;
-      r[p+2] = p+nf+1;
-      r[p+nf+1] = p+2;
-      for (int i = 0; i < of-1; i++) {
-        //switch edges 0 and 2
-        r[p+3+i] = p+3+3*(of-1)-(i+1)+nf;
-        r[p+3+3*(of-1)-(i+1)] = p+3+i+nf;
-        r[p+3+i+nf] = p+3+3*(of-1)-(i+1);
-        r[p+3+3*(of-1)-(i+1)+nf] = p+3+i;
-        //reverse edge 1
-        r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1)+nf;
-        r[p+3+(of-1)+i+nf] = p+3+2*(of-1)-(i+1);
-      }
-      p += 3*of;
-    }
-    p+=nf;
-
-    //F3(=312) switches its nodes 1 and 2
-    for (int of = layerOrder-3; of >= 0; of -= 3) {
-      r[p] = p;
-      if (of == 0) {
-        p += 1;
-        break;
   }
-      r[p+1] = p+2;
-      r[p+2] = p+1;
-      for (int i = 0; i < of-1; i++) {
-        //switch edges 0 and 2
-        r[p+3+i] = p+3+3*(of-1)-(i+1);
-        r[p+3+3*(of-1)-(i+1)] = p+3+i;
-        //reverse edge 1
-        r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
 }
-      p += 3*of;
+
+void MTetrahedronN::reverse()
+{
+  std::map<int, indicesReversed>::iterator it;
+  it = _order2indicesReversed.find(_order);
+  if (it == _order2indicesReversed.end()) {
+    indicesReversed indices;
+    _getIndicesReversed(_order, indices);
+    _order2indicesReversed[_order] = indices;
+    it = _order2indicesReversed.find(_order);
+  }
+
+  indicesReversed &indices = it->second;
+
+  // copy vertices
+  std::vector<MVertex*> oldv(4 + _vs.size());
+  std::copy(_v, _v+4, oldv.begin());
+  std::copy(_vs.begin(), _vs.end(), oldv.begin()+4);
+
+  // reverse
+  for (int i = 0; i < 4; ++i) {
+    _v[i] = oldv[indices[i]];
   }
+  for (int i = 0; i < _vs.size(); ++i) {
+    _vs[i] = oldv[indices[4+i]];
   }
-  return r;
 }
\ No newline at end of file

Geo/MTetrahedron.h

@@ -309,8 +309,11 @@ class MTetrahedron10 : public MTetrahedron {
  *
  */
 
+typedef std::vector<int> indicesReversed;
+
 class MTetrahedronN : public MTetrahedron {
-  static const std::vector<int> &_getReverseIndices(int order);
+  static std::map<int, indicesReversed> _order2indicesReversed;
+
  protected:
   std::vector<MVertex *> _vs;
   const char _order;
@@ -412,16 +415,7 @@ class MTetrahedronN : public MTetrahedron {
     Msg::Error("no tag matches a p%d tetrahedron with %d vertices", _order, 4+_vs.size());
     return 0;
   }
-  virtual void reverse()
-  {
-    MVertex *tmp;
-    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
-    std::vector<MVertex*> inv(_vs.size());
-    std::vector<int> reverseIndices = _getReverseIndices(_order);
-    for (unsigned int i = 0; i< _vs.size(); i++)
-      inv[i] = _vs[reverseIndices[i + 4] - 4];
-    _vs = inv;
-  }
+  virtual void reverse();
   virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
   virtual int getNumEdgesRep(bool curved);
   virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);