From 5c19ac96c49d6829737a722dba46994a7710f584 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 3 Sep 2010 07:14:37 +0000
Subject: [PATCH] use inexact parFromPoint in mesh orientation

---
 Geo/ACISFace.cpp              |  2 +-
 Geo/ACISFace.h                |  2 +-
 Geo/GFace.cpp                 | 17 +++++++++--------
 Geo/GFace.h                   |  7 +++----
 Geo/GFaceCompound.cpp         |  2 +-
 Geo/GFaceCompound.h           |  2 +-
 Geo/MVertex.cpp               |  9 +++++----
 Geo/MVertex.h                 |  6 +++---
 Geo/OCCFace.cpp               |  2 +-
 Geo/OCCFace.h                 |  2 +-
 Geo/discreteFace.cpp          | 13 +------------
 Geo/discreteFace.h            |  2 +-
 Geo/fourierFace.cpp           |  2 +-
 Geo/fourierFace.h             |  2 +-
 Geo/fourierProjectionFace.cpp | 20 ++++++++++----------
 Geo/fourierProjectionFace.h   |  2 +-
 Geo/gmshFace.cpp              |  4 ++--
 Geo/gmshFace.h                |  2 +-
 Mesh/meshGFace.cpp            |  4 +++-
 Mesh/meshGRegion.cpp          |  6 +-----
 20 files changed, 48 insertions(+), 60 deletions(-)

diff --git a/Geo/ACISFace.cpp b/Geo/ACISFace.cpp
index f7bc632d3e..d5a5b5b67d 100644
--- a/Geo/ACISFace.cpp
+++ b/Geo/ACISFace.cpp
@@ -148,7 +148,7 @@ GPoint ACISFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) c
   return GPoint(fpt.x(),fpt.y(),fpt.z(),this,pt2);
 }
 
-SPoint2 ACISFace::parFromPoint(const SPoint3 &qp) const
+SPoint2 ACISFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
 {
   SURFACE *Surf = _f->geometry();
   SPAposition pt(qp.x(),qp.y(),qp.z());
diff --git a/Geo/ACISFace.h b/Geo/ACISFace.h
index 15b1bc579f..fa834b4106 100644
--- a/Geo/ACISFace.h
+++ b/Geo/ACISFace.h
@@ -37,7 +37,7 @@ class ACISFace : public GFace {
   virtual GEntity::GeomType geomType() const; 
   ModelType getNativeType() const { return AcisModel; }
   void *getNativePtr() const { return (void*)_f; }
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
   virtual double curvatureMax(const SPoint2 &param) const;
   virtual double curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
                             double *curvMax, double *curvMin) const;
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index f020352fc3..5b2a64f787 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -779,12 +779,11 @@ void GFace::getMetricEigenVectors(const SPoint2 &param,
   }
 }
 
-void GFace::XYZtoUV(const double X, const double Y, const double Z,
-                    double &U, double &V, const double relax,
-                    const bool onSurface) const
+void GFace::XYZtoUV(double X, double Y, double Z, double &U, double &V,
+                    double relax, bool onSurface) const
 {
-  const double Precision = 1.e-8;
-  const int MaxIter = 25;
+  const double Precision = onSurface ? 1.e-8 : 1.e-3;
+  const int MaxIter = onSurface ? 25 : 10;
   const int NumInitGuess = 11;
 
   double Unew = 0., Vnew = 0., err, err2;
@@ -855,6 +854,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
     }
   }
 
+  if(!onSurface) return;
+
   if(relax < 1.e-6)
     Msg::Error("Could not converge: surface mesh will be wrong");
   else {
@@ -863,14 +864,14 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
   }
 }
 
-SPoint2 GFace::parFromPoint(const SPoint3 &p) const
+SPoint2 GFace::parFromPoint(const SPoint3 &p, bool onSurface) const
 {
   double U = 0., V = 0.;
-  XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0);
+  XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0, onSurface);
   return SPoint2(U, V);
 }
 
-GPoint GFace::closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const
+GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
 {
   Msg::Error("Closest point not implemented for this type of surface");
   return GPoint(0, 0, 0);
diff --git a/Geo/GFace.h b/Geo/GFace.h
index f7ae16463a..e1c39e80ba 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -105,9 +105,8 @@ class GFace : public GEntity
   virtual void setVisibility(char val, bool recursive=false);
 
   // compute the parameters UV from a point XYZ
-  void XYZtoUV(const double X, const double Y, const double Z,
-               double &U, double &V, const double relax,
-               const bool onSurface=true) const;
+  void XYZtoUV(double X, double Y, double Z, double &U, double &V,
+               double relax, bool onSurface=true) const;
 
   // get the bounding box
   virtual SBoundingBox3d bounds() const;
@@ -158,7 +157,7 @@ class GFace : public GEntity
 
   // return the parmater location on the face given a point in space
   // that is on the face
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
 
   // true if the parameter value is interior to the face
   virtual bool containsParam(const SPoint2 &pt) const;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index f3411af952..292deca07e 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1418,7 +1418,7 @@ double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
 
 }
 
-SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p) const
+SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
 {
   if(!oct) parametrize();
 
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 095309577e..1e43c0fd7c 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -106,7 +106,7 @@ class GFaceCompound : public GFace {
   { return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
   virtual GPoint point(double par1, double par2) const; 
   typeOfMapping getTypeOfMapping() { return _mapping;}
-  SPoint2 parFromPoint(const SPoint3 &p) const;
+  SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
   virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const;
   virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const; 
   virtual GEntity::GeomType geomType() const { return CompoundSurface; }
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 5f5981e9ee..ac1a6a256b 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -326,7 +326,8 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
   }
 }
 
-bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
+bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param,
+                             bool onSurface)
 {
   if (gf->geomType() == GEntity::CompoundSurface &&
       v->onWhat()->dim() < 2){
@@ -337,7 +338,7 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
 
   if(v->onWhat()->geomType() == GEntity::DiscreteCurve ||        
      v->onWhat()->geomType() == GEntity::BoundaryLayerCurve){    
-    param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+    param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()), onSurface);
     return true;
   }
 
@@ -345,7 +346,7 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
     GVertex *gv = (GVertex*)v->onWhat();
     // hack for bug in periodic curves
     if (gv->getNativeType() == GEntity::GmshModel && gf->geomType() == GEntity::Plane)
-      param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+      param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()), onSurface);
     else
       param = gv->reparamOnFace(gf, 1);
     // shout, we could be on a seam
@@ -370,7 +371,7 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
     }
     else {
       // brute force!
-      param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+      param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()), onSurface);
     }
   }
   return true;
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 12abb066e4..40c8d4fd93 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -127,8 +127,7 @@ class MEdgeVertex : public MVertex{
   {
   }
   virtual ~MEdgeVertex(){}
-  virtual bool getParameter(int i, double &par) const { 
-    par = _u; return true; }
+  virtual bool getParameter(int i, double &par) const { par = _u; return true; }
   virtual bool setParameter(int i, double par){ _u = par; return true; }
   double getLc() const { return _lc; }
 };
@@ -155,7 +154,8 @@ class MFaceVertex : public MVertex{
 
 bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
                            SPoint2 &param1, SPoint2 &param2);
-bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param);
+bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param,
+                             bool onSurface=true);
 bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double &param);
 
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 10213ec439..b72cdb0241 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -174,7 +174,7 @@ GPoint OCCFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) co
   return GPoint(pnt.X(), pnt.Y(), pnt.Z(), this, pp);
 }
 
-SPoint2 OCCFace::parFromPoint(const SPoint3 &qp) const
+SPoint2 OCCFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
 {
   gp_Pnt pnt(qp.x(), qp.y(), qp.z());
   GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 6a9227f5f5..453602116f 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -39,7 +39,7 @@ class OCCFace : public GFace {
   virtual GEntity::GeomType geomType() const; 
   ModelType getNativeType() const { return OpenCascadeModel; }
   void *getNativePtr() const { return (void*)&s; }
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
   virtual double curvatureMax(const SPoint2 &param) const;
   virtual double curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMin,
                             double *curvMax, double *curvMin) const;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index d07aeb3e8a..f85b47ce99 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -23,7 +23,6 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
 
 void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_edges)
 {
-
   std::set<MEdge, Less_Edge> bound_edges;
   for (unsigned int iFace = 0; iFace  < getNumMeshElements() ; iFace++) {
     MElement *e = getMeshElement(iFace);
@@ -50,7 +49,6 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
       itmap->second = tagFaces;
     }
   }
-
 }
 
 void discreteFace::setBoundEdges(std::vector<int> tagEdges)
@@ -69,9 +67,8 @@ GPoint discreteFace::point(double par1, double par2) const
   return GPoint();
 }
 
-SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
+SPoint2 discreteFace::parFromPoint(const SPoint3 &p, bool onSurface) const
 {
-
   if (getCompound()){
     return getCompound()->parFromPoint(p);
   }
@@ -83,7 +80,6 @@ SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
 
 SVector3 discreteFace::normal(const SPoint2 &param) const
 {
-
   if (getCompound()){
     return getCompound()->normal(param);
   }
@@ -91,12 +87,10 @@ SVector3 discreteFace::normal(const SPoint2 &param) const
     Msg::Error("Cannot evaluate normal on discrete face");
     return SVector3();
   }
-
 }
 
 double discreteFace::curvatureMax(const SPoint2 &param) const
 {
-
   if (getCompound()){
     return getCompound()->curvatureMax(param);
   }
@@ -104,12 +98,10 @@ double discreteFace::curvatureMax(const SPoint2 &param) const
     Msg::Error("Cannot evaluate curvature on discrete face");
     return false;
   }
-
 }
 
 Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 &param) const
 {
-
   if (getCompound()){
     return getCompound()->firstDer(param);
   }
@@ -117,13 +109,11 @@ Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 &param) const
     Msg::Error("Cannot evaluate derivative on discrete face");
     return Pair<SVector3, SVector3>();
   }
-
 }
 
 void discreteFace::secondDer(const SPoint2 &param,
                              SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
 {
-
   if (getCompound()){
     return getCompound()->secondDer(param, dudu, dvdv, dudv);
   }
@@ -131,5 +121,4 @@ void discreteFace::secondDer(const SPoint2 &param,
     Msg::Error("Cannot evaluate second derivative on discrete face");
     return;
   }
-
 }
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 14952b1d7f..db16cfa9a4 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -16,7 +16,7 @@ class discreteFace : public GFace {
   discreteFace(GModel *model, int num);
   virtual ~discreteFace() {}
   GPoint point(double par1, double par2) const;
-  SPoint2 parFromPoint(const SPoint3 &p) const;
+  SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
   SVector3 normal(const SPoint2 &param) const;
   double curvatureMax(const SPoint2 &param) const;
   GEntity::GeomType geomType() const { return DiscreteSurface; }
diff --git a/Geo/fourierFace.cpp b/Geo/fourierFace.cpp
index bf771cfb8b..569f600bcc 100644
--- a/Geo/fourierFace.cpp
+++ b/Geo/fourierFace.cpp
@@ -36,7 +36,7 @@ GPoint fourierFace::point(double par1, double par2) const
   return GPoint(x, y, z, this, pp);
 }
 
-SPoint2 fourierFace::parFromPoint(const SPoint3 &p) const
+SPoint2 fourierFace::parFromPoint(const SPoint3 &p, bool onSurface) const
 {
   double u, v, x, y, z;
   x = p.x(); y = p.y(); z = p.z();
diff --git a/Geo/fourierFace.h b/Geo/fourierFace.h
index e6d31db6a9..06bdd29f23 100644
--- a/Geo/fourierFace.h
+++ b/Geo/fourierFace.h
@@ -24,7 +24,7 @@ class fourierFace : public GFace {
   virtual ~fourierFace() {}
   Range<double> parBounds(int i) const; 
   virtual GPoint point(double par1, double par2) const; 
-  virtual SPoint2 parFromPoint(const SPoint3 &p) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
   virtual bool containsParam(const SPoint2 &pt) const; 
   virtual SVector3 normal(const SPoint2 &param) const; 
   virtual GEntity::GeomType geomType() const;
diff --git a/Geo/fourierProjectionFace.cpp b/Geo/fourierProjectionFace.cpp
index ee3eb3e3e0..20cdb7c114 100644
--- a/Geo/fourierProjectionFace.cpp
+++ b/Geo/fourierProjectionFace.cpp
@@ -21,24 +21,24 @@ fourierProjectionFace::~fourierProjectionFace() {}
 GPoint fourierProjectionFace::point(double par1, double par2) const
 {
   SVector3 p;
-  ps_->F(par1,par2,p[0],p[1],p[2]);
-  return GPoint(p[0],p[1],p[2]);
+  ps_->F(par1, par2, p[0], p[1], p[2]);
+  return GPoint(p[0], p[1], p[2]);
 }
 
-SPoint2 fourierProjectionFace::parFromPoint(const SPoint3 &p) const
+SPoint2 fourierProjectionFace::parFromPoint(const SPoint3 &p, bool onSurface) const
 {  
   double u,v;
-  ps_->Inverse(p[0],p[1],p[2],u,v);
-  return SPoint2(u,v);
+  ps_->Inverse(p[0], p[1], p[2], u, v);
+  return SPoint2(u, v);
 }
 
 Pair<SVector3,SVector3> fourierProjectionFace::firstDer(const SPoint2 &param) const
 {
   SVector3 du;
   SVector3 dv;
-  ps_->Dfdu(param.x(),param.y(),du[0],du[1],du[2]);
-  ps_->Dfdv(param.x(),param.y(),dv[0],dv[1],dv[2]);
-  return Pair<SVector3,SVector3>(du,dv);
+  ps_->Dfdu(param.x(), param.y(), du[0], du[1], du[2]);
+  ps_->Dfdv(param.x(), param.y(), dv[0], dv[1], dv[2]);
+  return Pair<SVector3,SVector3>(du, dv);
 } 
 
 void fourierProjectionFace::secondDer(const SPoint2 &param, 
@@ -50,8 +50,8 @@ void fourierProjectionFace::secondDer(const SPoint2 &param,
 SVector3 fourierProjectionFace::normal(const SPoint2 &param) const
 {       
   double x, y, z;
-  ps_->GetUnitNormal(param.x(),param.y(),x,y,z);
-  return SVector3(x,y,z);
+  ps_->GetUnitNormal(param.x(), param.y(), x, y, z);
+  return SVector3(x, y, z);
 }
 
 Range<double> fourierProjectionFace::parBounds(int i) const
diff --git a/Geo/fourierProjectionFace.h b/Geo/fourierProjectionFace.h
index fca1a5abaa..5987a7428b 100644
--- a/Geo/fourierProjectionFace.h
+++ b/Geo/fourierProjectionFace.h
@@ -25,7 +25,7 @@ class fourierProjectionFace : public GFace {
   SVector3 normal(const SPoint2 &param) const; 
   Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const; 
   void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const; 
-  SPoint2 parFromPoint(const SPoint3 &) const;
+  SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
   virtual GEntity::GeomType geomType() const { return GEntity::ProjectionFace; }
   ModelType getNativeType() const { return UnknownModel; }
   void *getNativePtr() const { return ps_; } 
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 0b152f8df4..1f47b024e2 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -258,7 +258,7 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
   return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
 }
 
-SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
+SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
 {
   if(s->Typ == MSH_SURF_PLAN){
     double x, y, z, VX[3], VY[3];
@@ -269,7 +269,7 @@ SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
     return SPoint2(u, v); 
   }
   else{
-    return GFace::parFromPoint(qp);
+    return GFace::parFromPoint(qp, onSurface);
   }
 }
 
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index ad6c31abc1..fa229909c9 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -32,7 +32,7 @@ class gmshFace : public GFace {
   virtual GEntity::GeomType geomType() const; 
   ModelType getNativeType() const { return GmshModel; }
   void *getNativePtr() const { return s; }
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
+  virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
   virtual void resetMeshAttributes();
 };
 
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 31730f2ca0..59b3724d76 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1814,6 +1814,7 @@ void orientMeshGFace::operator()(GFace *gf)
   //   do not seem to be consistent with the orientation of the
   //   bounding edges
 
+
   // first, try to find an element with one vertex categorized on the
   // surface and for which we have valid surface parametric
   // coordinates
@@ -1844,7 +1845,7 @@ void orientMeshGFace::operator()(GFace *gf)
     bool ok = true;
     for(int j = 0; j < e->getNumVertices(); j++){
       SPoint2 p;
-      bool ok = reparamMeshVertexOnFace(e->getVertex(j), gf, p);
+      bool ok = reparamMeshVertexOnFace(e->getVertex(j), gf, p, false);
       if(!ok) break;
       param += p;
     }
@@ -1860,5 +1861,6 @@ void orientMeshGFace::operator()(GFace *gf)
       return;
     }
   }
+
   Msg::Warning("Could not orient mesh in face %d", gf->tag());
 }
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index a14b8b59cb..de64974bef 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -326,14 +326,12 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
   // everything
   gr->model()->destroyMeshCaches();
   std::list<GFace*> faces = gr->faces();
-  std::list<GFace*>::iterator it = faces.begin();
-  while(it != faces.end()){
+  for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
     GFace *gf = (*it); 
     for(unsigned int i = 0; i < gf->triangles.size(); i++)
       delete gf->triangles[i];
     gf->triangles.clear();
     gf->deleteVertexArrays();
-    ++it;
   }
 
   // TODO: re-create 1D mesh
@@ -341,7 +339,6 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
     MVertex *v[2];
     v[0] = numberedV[out.edgelist[i * 2 + 0] - 1];
     v[1] = numberedV[out.edgelist[i * 2 + 1] - 1];
-
     //implement here the 1D mesh ...
   }
 
@@ -515,7 +512,6 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
       gr->set(faces);
       return;
     }
-    printf("transfer tetgen \n");
     TransferTetgenMesh(gr, in, out, numberedV);
   }
 
-- 
GitLab