From ac5628ccba7bcf0ef6f320d6373c4083968c3545 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 4 May 2014 05:14:50 +0000
Subject: [PATCH] fixes coverity errors

---
 Geo/GEdge.cpp                                 |  2 +-
 Geo/boundaryLayersData.cpp                    | 18 +++--
 Geo/closestPoint.cpp                          | 18 +++--
 Geo/closestPoint.h                            |  1 -
 Mesh/meshGRegion.cpp                          |  5 +-
 Plugin/CurvedBndDist.cpp                      | 31 ++++++--
 Plugin/CurvedBndDist.h                        | 24 +++----
 .../OptHomFastCurving.cpp                     |  4 +-
 .../OptHomIntegralBoundaryDist.cpp            | 17 ++---
 .../OptHomIntegralBoundaryDist.h              | 13 ++--
 contrib/HighOrderMeshOptimizer/OptHomRun.cpp  | 15 ++--
 contrib/HighOrderMeshOptimizer/SuperEl.cpp    | 70 ++++++++-----------
 contrib/HighOrderMeshOptimizer/SuperEl.h      | 22 ++----
 13 files changed, 119 insertions(+), 121 deletions(-)

diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 66dd62cc9c..ddbaec1283 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -17,7 +17,7 @@
 #include "closestPoint.h"
 
 GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
-  : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1), compound(0),_cp(0)
+  : GEntity(model, tag), _tooSmall(false), _cp(0), v0(_v0), v1(_v1), compound(0)
 {
   if(v0) v0->addEdge(this);
   if(v1 && v1 != v0) v1->addEdge(this);
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 72432af3f7..36c2638bd8 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -350,10 +350,12 @@ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
 	  _dirs.push_back(x);
 	}
       }
+      /*
       else {
 	_dirs.push_back(N1[SIDE]);
 	_dirs.push_back(N2[SIDE]);
 	}
+      */
     }
   }
 }
@@ -972,11 +974,13 @@ fanTopology :: fanTopology (GRegion * gr, const std::set<GEdge*> &detectedFans,
 // V = (x,y,z) / (x^2+y^2+z^2)^{1/2} is a point on the unit sphere
 // the n[i]'s are points on the unit sphere
 // V maximizes  min_i (V * n[i])
-// This means I'd like to find point V that is 
+// This means I'd like to find point V that is
 
-static void filterVectors(std::vector<SVector3> &n){
+/*
+static void filterVectors(std::vector<SVector3> &n)
+{
   std::vector<SVector3> temp;
-  temp.push_back(n[0]);  
+  temp.push_back(n[0]);
   for (unsigned int i = 1 ; i<n.size() ; i++){
     bool found = false;
     for (unsigned int j = 0 ; j<temp.size() ; j++){
@@ -987,8 +991,10 @@ static void filterVectors(std::vector<SVector3> &n){
   }
   n = temp;
 }
+*/
 
-static SVector3 computeBestNormal(std::vector<SVector3> &n){
+static SVector3 computeBestNormal(std::vector<SVector3> &n)
+{
   //  filterVectors (n);
   SVector3 V;
   if (n.size() == 1)V = n[0];
@@ -1056,13 +1062,13 @@ static int createColumnsBetweenFaces(GRegion *gr,
     joints.push_back(joint);
     //    SVector3 avg (0,0,0);
     std::vector<SVector3> ns;
-    
+
     for (unsigned int i=0;i<joint.size(); i++){
       ns.push_back( n[inv[joint[i]]] );
       //      avg += n[inv[joint[i]]];
     }
     SVector3 avg = computeBestNormal(ns);
-    
+
 
     std::vector<MVertex*> _column;
     std::vector<SMetric3> _metrics;
diff --git a/Geo/closestPoint.cpp b/Geo/closestPoint.cpp
index 07d00dac61..2b3816d7f3 100644
--- a/Geo/closestPoint.cpp
+++ b/Geo/closestPoint.cpp
@@ -21,18 +21,22 @@ static void oversample (std::vector<SPoint3> &s, double tol){
   s = t;
 }
 
-closestPointFinder :: closestPointFinder (GEntity *ge, double e) : _ge (ge) , _tolerance (e)
+closestPointFinder :: closestPointFinder (GEntity *ge, double e) : _tolerance (e)
 {
 #if defined(HAVE_ANN)
   std::vector<SPoint3> pts;
   if (ge->dim() == 1){
     GEdge *edge = ge->cast2Edge();
-    if (!edge)Msg::Fatal("in the constructor of closestPoint");
-    std::vector<double> ts;
-    edge->discretize(_tolerance, pts,ts);
-    //    printf("%d points\n",pts.size());
-    oversample (pts, _tolerance);
-    //    printf("%d points after oversampling\n",pts.size());
+    if (edge){
+      std::vector<double> ts;
+      edge->discretize(_tolerance, pts,ts);
+      //    printf("%d points\n",pts.size());
+      oversample (pts, _tolerance);
+      //    printf("%d points after oversampling\n",pts.size());
+    }
+    else{
+      Msg::Error("Can get edge in closestPointFinder");
+    }
   }
   index = new ANNidx[1];
   dist = new ANNdist[1];
diff --git a/Geo/closestPoint.h b/Geo/closestPoint.h
index 24ecc73fa4..b80b6f9850 100644
--- a/Geo/closestPoint.h
+++ b/Geo/closestPoint.h
@@ -14,7 +14,6 @@ class closestPointFinder
   ANNidxArray index;
   ANNdistArray dist;
 #endif
-  GEntity* _ge;
   double _tolerance;
  public :
   closestPointFinder (GEntity*, double);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 74fe7e0f27..95fcdffc24 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1040,7 +1040,8 @@ static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr, spl
 	}
       }
       //      printf("counter = %d\n",counter);
-      opposite /= (double)counter;
+      if(counter)
+        opposite /= (double)counter;
 
       SVector3 dir = center - opposite;
       MTriangle temp (it->first.getVertex(0),it->first.getVertex(1),it->first.getVertex(2));
@@ -1833,7 +1834,7 @@ GFace *findInFaceSearchStructure(MVertex *p1, MVertex *p2, MVertex *p3,
   return it->second;
 }
 
-GFace *findInFaceSearchStructure(const MFace &ff, 
+GFace *findInFaceSearchStructure(const MFace &ff,
                                  const fs_cont &search)
 {
   fs_cont::const_iterator it = search.find(ff);
diff --git a/Plugin/CurvedBndDist.cpp b/Plugin/CurvedBndDist.cpp
index 6de56777d3..9b7b90752e 100644
--- a/Plugin/CurvedBndDist.cpp
+++ b/Plugin/CurvedBndDist.cpp
@@ -4,13 +4,17 @@
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "Gmsh.h"
+#include "GmshConfig.h"
 #include "GModel.h"
-#include "CurvedBndDist.h"
 #include "MElement.h"
 #include "PViewDataList.h"
-#include "OptHomIntegralBoundaryDist.h"
 #include "BasisFactory.h"
 
+#if defined(HAVE_OPTHOM)
+#include "CurvedBndDist.h"
+#include "OptHomIntegralBoundaryDist.h"
+#endif
+
 StringXNumber CurvedBndDistOptions_Number[] = {
  /* {GMSH_FULLRC, "Dim", NULL, -1},*/
 };
@@ -56,6 +60,7 @@ static void addLine(PViewDataList *data, const SVector3 &p0, const SVector3 &p1,
   data->SQ.push_back(v); data->SQ.push_back(v); data->SQ.push_back(v); data->SQ.push_back(v);
 }
 */
+
 static void addPoint(PViewDataList *data, const SVector3 &p0, double v0)
 {
   data->NbSP ++;
@@ -76,8 +81,11 @@ static void addVP(PViewDataList *data, const SVector3 p0, SVector3 v)
   data->VP.push_back(v.z());
 }
 
+#if defined(HAVE_OPTHOM)
+
 #include <limits>
-static void drawElementDist(PViewDataList *data, GEdge *edge, const std::vector<MVertex *>&vertices, const nodalBasis &basis)
+static void drawElementDist(PViewDataList *data, GEdge *edge,
+                            const std::vector<MVertex *>&vertices, const nodalBasis &basis)
 {
   std::vector<double> gradient;
   std::vector<double> param(vertices.size());
@@ -98,7 +106,8 @@ static void drawElementDist(PViewDataList *data, GEdge *edge, const std::vector<
   }
 }
 
-static void drawElement(PViewDataList *data, GEdge *edge, const std::vector<MVertex *>&vertices, const nodalBasis &basis)
+static void drawElement(PViewDataList *data, GEdge *edge,
+                        const std::vector<MVertex *>&vertices, const nodalBasis &basis)
 {
   std::vector<double> gradient;
   std::vector<double> param(vertices.size());
@@ -124,13 +133,16 @@ static void drawElement(PViewDataList *data, GEdge *edge, const std::vector<MVer
   }
 }
 
+#endif
+
 PView *GMSH_CurvedBndDistPlugin::execute(PView *v)
 {
+#if defined(HAVE_OPTHOM)
   PView *pv = new PView();
   PViewDataList *data = getDataList(pv);
   data->Time.push_back(0.);
-  _m = GModel::current();
-  for (GModel::fiter iface = _m->firstFace(); iface != _m->lastFace(); ++iface) {
+  GModel *m = GModel::current();
+  for (GModel::fiter iface = m->firstFace(); iface != m->lastFace(); ++iface) {
     GFace *face = *iface;
     for (size_t iElement = 0; iElement < face->getNumMeshElements(); ++iElement) {
       MElement *element = face->getMeshElement(iElement);
@@ -148,12 +160,17 @@ PView *GMSH_CurvedBndDistPlugin::execute(PView *v)
           }
         }
         if (edge){
-	  drawElementDist(data, edge, vertices, *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)));
+	  drawElementDist(data, edge, vertices,
+                          *BasisFactory::getNodalBasis(elbasis.getClosureType(clId)));
 	}
       }
     }
   }
   data->finalize();
   return pv;
+#else
+  Msg::Error("Plugin requires OPTHOM module");
+  return v;
+#endif
 }
 
diff --git a/Plugin/CurvedBndDist.h b/Plugin/CurvedBndDist.h
index b346dd788e..416435bc30 100644
--- a/Plugin/CurvedBndDist.h
+++ b/Plugin/CurvedBndDist.h
@@ -15,17 +15,17 @@ extern "C"
 
 class GMSH_CurvedBndDistPlugin : public GMSH_PostPlugin
 {
-    GModel *_m;
-  public :
-    GMSH_CurvedBndDistPlugin(){}
-    std::string getName() const { return "CurvedBndDist"; }
-    std::string getShortHelp() const {
-      return "Check validity of elements and/or compute min&max of the jacobian";
-    }
-    std::string getHelp() const;
-    std::string getAuthor() const { return "Jonathan Lambrechts"; }
-    int getNbOptions() const;
-    StringXNumber *getOption(int);  
-    PView *execute(PView *);
+ public :
+  GMSH_CurvedBndDistPlugin(){}
+  std::string getName() const { return "CurvedBndDist"; }
+  std::string getShortHelp() const
+  {
+    return "Compute distance to curved boundary";
+  }
+  std::string getHelp() const;
+  std::string getAuthor() const { return "Jonathan Lambrechts"; }
+  int getNbOptions() const;
+  StringXNumber *getOption(int);
+  PView *execute(PView *);
 };
 #endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 20315c69fc..15474d5041 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -197,7 +197,7 @@ MElement *getFirstEl(const std::map<MVertex*, std::vector<MElement*> > &vertex2e
 
   const std::vector<MElement*> &nb0 = vertex2elements.find(vb0)->second;
   const std::vector<MElement*> &nb1 = vertex2elements.find(vb1)->second;
-  MElement *firstEl;
+  MElement *firstEl = 0;
   for (int iEl=0; iEl<nb0.size(); iEl++) {
     if (find(nb1.begin(), nb1.end(), nb0[iEl]) != nb1.end()) {
       firstEl = nb0[iEl];
@@ -316,7 +316,7 @@ bool getTopPrimVertices(std::map<MVertex*, std::vector<MElement *> > &vertex2ele
       std::map<MVertex*, SVector3>::const_iterator itNorm = normVert.find(baseVert[i]);
       if (itNorm == normVert.end()) {
         Msg::Error("Normal to vertex not found in getTopPrimVertices");
-        itNorm == normVert.begin();
+        itNorm = normVert.begin();
       }
       std::vector<MVertex*> col;
       bool colFound = findColumn(vertex2elements, baseVert[i], itNorm->second,
diff --git a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
index 5f3c02f2d6..954302a0a0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.cpp
@@ -7,9 +7,9 @@
 #include "discreteFrechetDistance.h"
 #include "MLine.h"
 
-
-parametricLineNodalBasis::parametricLineNodalBasis(const nodalBasis &basis, const std::vector<SPoint3> &xyz):
-  _basis(basis), _xyz(xyz)  {};
+parametricLineNodalBasis::parametricLineNodalBasis(const nodalBasis &basis,
+                                                   const std::vector<SPoint3> &xyz):
+  _basis(basis), _xyz(xyz) {};
 
 SPoint3 parametricLineNodalBasis::operator()(double xi) const
 {
@@ -50,8 +50,6 @@ SVector3 parametricLineNodalBasis::secondDerivative(double xi) const
   return p;
 }
 
-
-
 parametricLineGEdge::parametricLineGEdge(const GEdge *edge, double t0, double t1):
   _edge(edge), _t0(t0), _t1(t1) {}
 
@@ -60,6 +58,7 @@ SPoint3 parametricLineGEdge::operator()(double xi) const
   GPoint gp = _edge->point(_t0 + (_t1 - _t0) * xi);
   return SPoint3 (gp.x(), gp.y(), gp.z());
 }
+
 SVector3 parametricLineGEdge::derivative(double xi) const
 {
   return _edge->firstDer(_t0 + (_t1 - _t0) * xi);
@@ -70,9 +69,9 @@ SVector3 parametricLineGEdge::secondDerivative(double xi) const
   return _edge->secondDer(_t0 + (_t1 - _t0) * xi);
 }
 
-static void oversample (std::vector<SPoint3> &s, double tol){
+static void oversample (std::vector<SPoint3> &s, double tol)
+{
   std::vector<SPoint3> t;
-
   for (unsigned int i=1;i<s.size();i++){
     SPoint3 p0 = s[i-1];
     SPoint3 p1 = s[i];
@@ -193,7 +192,6 @@ double parametricLine::hausdorffDistance(const parametricLine &l, SPoint3 &p1, S
   }
 }
 
-
 // DISCRETE FRECHET DISTANCE
 double computeBndDistF(GEdge *edge,
 		       std::vector<double> & params, // the model edge
@@ -208,7 +206,6 @@ double computeBndDistF(GEdge *edge,
   return l1.frechetDistance(l2,p1,p2,tolerance);
 }
 
-
 // GMSH's DISTANCE
 /*
  */
@@ -255,8 +252,6 @@ double computeBndDistG_(GEdge *edge, std::vector<double> & p, // the model edge
 
   //  printf("computing diustance with tolerance %g\n",tolerac);
 
-
-
   double D = 0.0;
   const double U0 = basis.points(0,0);
   const double U1 = basis.points(1,0);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
index a0e22b57fe..3a9e5db2ce 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
@@ -7,10 +7,10 @@ class SPoint3;
 class MElement;
 class MVertex;
 double computeBndDist(MElement *element, int distanceDefinition, double tolerance);
-double computeBndDistAndGradient(GEdge *edge, 
-				 std::vector<double> &param, 
-				 const std::vector<MVertex*> &vs, 
-				 const nodalBasis &basis, std::vector<SPoint3> &xyz, 
+double computeBndDistAndGradient(GEdge *edge,
+				 std::vector<double> &param,
+				 const std::vector<MVertex*> &vs,
+				 const nodalBasis &basis, std::vector<SPoint3> &xyz,
 				 std::vector<bool> &onEdge, std::vector<double> &grad, double tolerance);
 
 class parametricLine {
@@ -22,10 +22,10 @@ class parametricLine {
   void discretize(std::vector<SPoint3> &dpts,std::vector<double> &ts,double Prec,
 		  double t0 = 0.0, double t1 = 1.0) const;
   void recur_discretize(const double &t1, const double &t2,
-			const SPoint3 &p1, const SPoint3 &p2, 
+			const SPoint3 &p1, const SPoint3 &p2,
 			std::vector<SPoint3> &dpts,
 			std::vector<double> &ts,
-			double Prec, int depth) const;  
+			double Prec, int depth) const;
   double frechetDistance(const parametricLine &l, SPoint3 &p1, SPoint3 &p2, double tol = 1.e-6) const;
   double hausdorffDistance(const parametricLine &l, SPoint3 &p1, SPoint3 &p2, double tol = 1.e-6) const;
 };
@@ -34,7 +34,6 @@ class parametricLineNodalBasis : public parametricLine
 {
   const nodalBasis &_basis;
   const std::vector<SPoint3> &_xyz;
-  double _u0, _u1;
  public :
   parametricLineNodalBasis(const nodalBasis &basis, const std::vector<SPoint3> &xyz);
   virtual SPoint3 operator()(double xi) const;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index ab0a0f7921..3630ae260b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -625,14 +625,15 @@ static void optimizeOneByOne
 #endif
 
 #include "OptHomIntegralBoundaryDist.h"
-double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double tolerance){
+double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double tolerance)
+{
   double maxd = 0.0;
   double sum = 0.0;
-  int NUM;
+  int NUM = 0;
   for (int iEl = 0; iEl < ge->getNumMeshElements();iEl++) {
     MElement *el = ge->getMeshElement(iEl);
     if (ge->dim() == el->getDim()){
-      const double DISTE =computeBndDist(el,distanceDefinition, tolerance); 
+      const double DISTE =computeBndDist(el,distanceDefinition, tolerance);
       if (DISTE != 0.0){
 	NUM++;
 	//	if(distanceDefinition == 1)printf("%d %12.5E\n",iEl,DISTE);
@@ -641,7 +642,7 @@ double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double t
       }
     }
   }
-  if (distanceDefinition == 2)return sum/NUM;
+  if (distanceDefinition == 2 && NUM) return sum / (double)NUM;
   return maxd;
 }
 
@@ -672,17 +673,17 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
     for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
       double jmin, jmax;
       MElement *el = entity->getMeshElement(iEl);
-      const double DISTE =computeBndDist(el,2,fabs(p.discrTolerance)); 
+      const double DISTE =computeBndDist(el,2,fabs(p.discrTolerance));
       //      printf("Element %d Distance %12.5E\n",iEl,DISTE);
       maxdist = std::max(DISTE, maxdist);
       if (el->getDim() == p.dim) {
 	if (p.optCAD && DISTE > p.optCADDistMax)
-	  badasses.insert(el); 
+	  badasses.insert(el);
 
 	el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0);
 	if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
 	if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX){
-	  badasses.insert(el); 
+	  badasses.insert(el);
         }
       }
     }
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.cpp b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
index 59382a442a..d7ec149660 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.cpp
+++ b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
@@ -29,6 +29,7 @@
 
 #include <iostream>
 #include <sstream>
+#include "GmshMessage.h"
 #include "GEdge.h"
 #include "MLine.h"
 #include "MQuadrangle.h"
@@ -37,13 +38,10 @@
 #include "BasisFactory.h"
 #include "SuperEl.h"
 
-
-
 std::map<int,SuperEl::superInfoType> SuperEl::_superInfo;
 
-
-SuperEl::superInfoType::superInfoType(int type, int order) {
-
+SuperEl::superInfoType::superInfoType(int type, int order)
+{
   int iBaseFace = 0, iTopFace = 0;
   switch (type) {
     case TYPE_QUA:
@@ -56,7 +54,7 @@ SuperEl::superInfoType::superInfoType(int type, int order) {
       iBaseFace = 0; iTopFace = 5;
       break;
     default:
-      std::cout << "ERROR: SuperEl not implemented for element of type " << type << std::endl;
+      Msg::Error("SuperEl not implemented for element of type %d", type);
       nV = 0;
       return;
   }
@@ -64,27 +62,28 @@ SuperEl::superInfoType::superInfoType(int type, int order) {
   // Get HO nodal basis
   const int tag = ElementType::getTag(type, order, true);                 // Get tag for incomplete element type
 //  const int tag = ElementType::getTag(type, order);                     // Get tag for complete element type
-  const nodalBasis *basis = tag ? BasisFactory::getNodalBasis(tag) : 0;
-
-  nV = basis->getNumShapeFunctions();
-//  _superInfo[type].nV1 = basis->getNumShapeFunctions();
-  points = basis->points;
-
-  baseInd = basis->getClosure(basis->getClosureId(iBaseFace,1,0));
-  topInd = basis->getClosure(basis->getClosureId(iTopFace,0,0));
-  otherInd.reserve(nV-baseInd.size()-topInd.size());
-  for(int i=0; i<nV; ++i) {
-    const std::vector<int>::iterator inBaseFace = find(baseInd.begin(),baseInd.end(),i);
-    const std::vector<int>::iterator inTopFace = find(topInd.begin(),topInd.end(),i);
-    if (inBaseFace == baseInd.end() && inTopFace == topInd.end()) otherInd.push_back(i);
+  if(tag){
+    const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
+
+    nV = basis->getNumShapeFunctions();
+    //  _superInfo[type].nV1 = basis->getNumShapeFunctions();
+    points = basis->points;
+
+    baseInd = basis->getClosure(basis->getClosureId(iBaseFace,1,0));
+    topInd = basis->getClosure(basis->getClosureId(iTopFace,0,0));
+    otherInd.reserve(nV-baseInd.size()-topInd.size());
+    for(int i=0; i<nV; ++i) {
+      const std::vector<int>::iterator inBaseFace = find(baseInd.begin(),baseInd.end(),i);
+      const std::vector<int>::iterator inTopFace = find(topInd.begin(),topInd.end(),i);
+      if (inBaseFace == baseInd.end() && inTopFace == topInd.end()) otherInd.push_back(i);
+    }
   }
-
 }
 
 
 
 SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
-                 const std::vector<MVertex*> &topPrimVert)
+                 const std::vector<MVertex*> &topPrimVert) : _superEl(0), _superEl0(0)
 {
 
   // Get useful info on meta-element type if not already there
@@ -120,10 +119,8 @@ SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
       _superEl0 = new MHexahedron(_superVert);
       break;
     default:
-      std::cout << "ERROR: SuperEl not implemented for element of type " << type << std::endl;
-      _superEl0 = 0;
+      Msg::Error("SuperEl not implemented for element of type %d", type);
       return;
-      break;
   }
 
   // Add HO vertices in top face
@@ -183,30 +180,23 @@ SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
 
 }
 
-
-
 SuperEl::~SuperEl()
 {
   for (int i = 0; i < _superVert.size(); i++) delete _superVert[i];
   _superVert.clear();
   delete _superEl;
-  delete _superEl0;
+  if(_superEl0) delete _superEl0;
 }
 
-
-
-bool SuperEl::isPointIn(const SPoint3 p) const {
-
+bool SuperEl::isPointIn(const SPoint3 p) const
+{
   double xyz[3] = {p.x(), p.y(), p.z()}, uvw[3];
   _superEl0->xyz2uvw(xyz,uvw);
   return _superEl0->isInside(uvw[0],uvw[1],uvw[2]);
-
 }
 
-
-
-bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const {
-
+bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const
+{
   double uvw[3];
   _superEl0->xyz2uvw(xyzS,uvw);
   if (!_superEl0->isInside(uvw[0],uvw[1],uvw[2])) return false;
@@ -218,13 +208,10 @@ bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const {
   xyzC[2] = pC[2];
 
   return true;
-
 }
 
-
-
-std::string SuperEl::printPOS() {
-
+std::string SuperEl::printPOS()
+{
   std::vector<MVertex*> verts;
   _superEl->getVertices(verts);
 //  std::string posStr = _superEl->getStringForPOS();
@@ -243,5 +230,4 @@ std::string SuperEl::printPOS() {
   oss << "};\n";
 
   return oss.str();
-
 }
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.h b/contrib/HighOrderMeshOptimizer/SuperEl.h
index 54ba6991a6..b360495176 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.h
+++ b/contrib/HighOrderMeshOptimizer/SuperEl.h
@@ -33,45 +33,35 @@
 #include <string>
 #include "MElement.h"
 
-
-
 class SuperEl
 {
-public:
-
+ public:
   SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
           const std::vector<MVertex*> &topPrimVert);
   ~SuperEl();
-
   bool isOK() const { return _superEl; }
   bool isPointIn(const SPoint3 p) const;
   bool straightToCurved(double *xyzS, double *xyzC) const;
-
   std::string printPOS();
-
-  void printCoord() {
+  void printCoord()
+  {
     std::cout << "DBGTT: superEl -> ";
     for(int i = 0; i < _superVert.size(); i++){
-      std::cout << "v" << i << " = (" << _superVert[i]->x() << "," <<  _superVert[i]->y() << "," <<  _superVert[i]->z() << ")";
+      std::cout << "v" << i << " = (" << _superVert[i]->x() << ","
+                <<  _superVert[i]->y() << "," <<  _superVert[i]->z() << ")";
       if (i == _superVert.size()-1) std::cout << "\n"; else std::cout << ", ";
     }
-
   }
-
-private:
-
+ private:
   struct superInfoType {
     int nV;
     fullMatrix<double> points;
     std::vector<int> baseInd, topInd, otherInd;
     superInfoType(int type, int order);
   };
-
   static std::map<int,superInfoType> _superInfo;
-
   std::vector<MVertex*> _superVert;
   MElement *_superEl, *_superEl0;
 };
 
-
 #endif
-- 
GitLab