From 39437758e45dcc211f4543e24b49da4e89f7103f Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 22 Mar 2015 07:44:12 +0000
Subject: [PATCH] more fixes...

---
 Mesh/BGMBase.cpp          | 115 +++++++++++-------------
 Mesh/BGMBase.h            |  26 +++---
 Mesh/BackgroundMesh2D.cpp | 184 ++++++++++++++------------------------
 Mesh/pointInsertion.cpp   |  63 +++++--------
 Mesh/pointInsertion.h     |  66 ++++----------
 5 files changed, 166 insertions(+), 288 deletions(-)

diff --git a/Mesh/BGMBase.cpp b/Mesh/BGMBase.cpp
index 436c47430c..c4bc70d0e0 100644
--- a/Mesh/BGMBase.cpp
+++ b/Mesh/BGMBase.cpp
@@ -8,8 +8,8 @@
 #include "GmshDefines.h"
 #include "MElementOctree.h"
 
-void BGMBase::export_scalar(const string &filename, const DoubleStorageType &_whatToPrint) const{
-//  cout << "BGMBase::size of sizeField: " << sizeField.size()  << endl;
+void BGMBase::export_scalar(const std::string &filename, const DoubleStorageType &_whatToPrint) const
+{
   FILE *f = Fopen (filename.c_str(),"w");
   fprintf(f,"View \"Background Mesh\"{\n");
 
@@ -36,7 +36,7 @@ void BGMBase::export_scalar(const string &filename, const DoubleStorageType &_wh
 
     fprintf(f,"%s(",s);
     const MVertex *v;
-    vector<double> values(nvertex);
+    std::vector<double> values(nvertex);
     for (int iv=0;iv<nvertex;iv++){
       v = elem->getVertex(iv);
       values[iv] = get_nodal_value(v,_whatToPrint);
@@ -54,11 +54,10 @@ void BGMBase::export_scalar(const string &filename, const DoubleStorageType &_wh
   }
   fprintf(f,"};\n");
   fclose(f);
-//  cout << "export_scalar DONE " << endl;
 }
 
-
-void BGMBase::export_vector(const string &filename, const VectorStorageType &_whatToPrint)const{
+void BGMBase::export_vector(const std::string &filename, const VectorStorageType &_whatToPrint) const
+{
   FILE *f = Fopen (filename.c_str(),"w");
   fprintf(f,"View \"Background Mesh\"{\n");
 
@@ -85,10 +84,10 @@ void BGMBase::export_vector(const string &filename, const VectorStorageType &_wh
 
     fprintf(f,"%s(",s);
     const MVertex *v;
-    vector<double> values(nvertex*3);
+    std::vector<double> values(nvertex*3);
     for (int iv=0;iv<nvertex;iv++){
       v = elem->getVertex(iv);
-      vector<double> temp = get_nodal_value(v,_whatToPrint);
+      std::vector<double> temp = get_nodal_value(v,_whatToPrint);
       for (int j=0;j<3;j++)
         values[iv*3+j] = temp[j];
       GPoint p = get_GPoint_from_MVertex(v);
@@ -109,7 +108,8 @@ void BGMBase::export_vector(const string &filename, const VectorStorageType &_wh
 }
 
 
-void BGMBase::export_tensor_as_vectors(const string &filename, const TensorStorageType &_whatToPrint)const{
+void BGMBase::export_tensor_as_vectors(const std::string &filename, const TensorStorageType &_whatToPrint)const
+{
   FILE *f = Fopen (filename.c_str(),"w");
   fprintf(f,"View \"Background Mesh\"{\n");
 
@@ -129,91 +129,85 @@ void BGMBase::export_tensor_as_vectors(const string &filename, const TensorStora
 }
 
 
-BGMBase::BGMBase(int dim,GEntity *_gf):octree(NULL),gf(_gf), DIM(dim), order(1),debug(false){
-  if(debug) cout << "BGMBase::constructor " << endl;
+BGMBase::BGMBase(int dim,GEntity *_gf):octree(NULL),gf(_gf), DIM(dim), order(1)
+{
 }
 
-
 BGMBase::~BGMBase(){
 }
 
-
-bool BGMBase::inDomain (double u, double v, double w){
+bool BGMBase::inDomain (double u, double v, double w)
+{
   return (findElement(u, v, w) != NULL);
 }
 
-
-const MElement* BGMBase::findElement(double u, double v, double w, bool strict){
+const MElement* BGMBase::findElement(double u, double v, double w, bool strict)
+{
   return (getOctree()->find(u, v, w, DIM, strict));
 }
 
-
-vector<double> BGMBase::get_field_value(double u, double v, double w, const VectorStorageType &data){
+std::vector<double> BGMBase::get_field_value(double u, double v, double w, const VectorStorageType &data)
+{
   MElement *e = const_cast<MElement*>(findElement(u, v, w ));
-  if (!e) return vector<double>(3,-1000.);
-  vector<vector<double> > val = get_nodal_values(e,data);
-  vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,w);
+  if (!e) return std::vector<double>(3,-1000.);
+  std::vector<std::vector<double> > val = get_nodal_values(e,data);
+  std::vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,w);
 
-  vector<double> res(3);
+  std::vector<double> res(3);
   for (int j=0;j<3;j++){
-    double values[e->getNumVertices()];
+    std::vector<double> values(e->getNumVertices());
     for (int i=0;i<e->getNumVertices();i++) values[i]=val[i][j];
-    res[j] = e->interpolate(values, element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
+    res[j] = e->interpolate(&values[0], element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
   }
   return res;
 }
 
-
-double BGMBase::get_field_value(double u, double v, double w, const DoubleStorageType &data){
+double BGMBase::get_field_value(double u, double v, double w, const DoubleStorageType &data)
+{
   MElement *e = const_cast<MElement*>(findElement(u, v, w));
   if (!e) return -1000.;
-  vector<double> val = get_nodal_values(e,data);
-  vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,w);
-  double values[e->getNumVertices()];
+  std::vector<double> val = get_nodal_values(e,data);
+  std::vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,w);
+  std::vector<double> values(e->getNumVertices());
   for (int i=0;i<e->getNumVertices();i++)
     values[i]=val[i];
   
-  return e->interpolate(values, element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
+  return e->interpolate(&values[0], element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
 }
 
-
-double BGMBase::size(double u, double v, double w){
+double BGMBase::size(double u, double v, double w)
+{
   return get_field_value(u,v,w,sizeField);
 }
 
-
-double BGMBase::size(const MVertex *v){
+double BGMBase::size(const MVertex *v)
+{
   return get_nodal_value(v,sizeField);
 }
 
-
-vector<double> BGMBase::get_nodal_value(const MVertex *v,const VectorStorageType &data)const{
-  if(debug) cout << "BGMBase::get_nodal_value(const MVertex *v,const map<MVertex*,vector<double> > &data)" << endl;
+std::vector<double> BGMBase::get_nodal_value(const MVertex *v,const VectorStorageType &data)const
+{
   VectorStorageType::const_iterator itfind = data.find(const_cast<MVertex*>(v));
   if (itfind==data.end()){
-    cout << "WARNING: BGMBase::get_nodal_value (vector): unknown vertex ! " << v << endl;
-    throw;
-    return vector<double>(3,0.);
+    Msg::Error("Unknown vertex %d in BGMBase::get_nodal_value", v->getNum());
+    return std::vector<double>(3,0.);
   }
   return itfind->second;
 }
 
-
-double BGMBase::get_nodal_value(const MVertex *v,const DoubleStorageType &data)const{
-  if(debug) cout << "BGMBase::get_nodal_value(const MVertex *v,const map<MVertex*,double> &data)" << endl;
+double BGMBase::get_nodal_value(const MVertex *v,const DoubleStorageType &data)const
+{
   DoubleStorageType::const_iterator itfind = data.find(const_cast<MVertex*>(v));
   if (itfind==data.end()){
-    cout << "WARNING: BGMBase::get_nodal_value: unknown vertex ! " << endl;
-    throw;
+    Msg::Error("Unknown vertex %d in BGMBase::get_nodal_value", v->getNum());
     return 0.;
   }
   return itfind->second;
 }
 
-
-vector<vector<double> > BGMBase::get_nodal_values(const MElement *e,const VectorStorageType &data)const{
-  if(debug) cout << "BGMBase::get_nodal_values(const MElement *e,const map<MVertex*,vector<double> > &data)" << endl;
-  vector<vector<double> > res(e->getNumVertices());
+std::vector<std::vector<double> > BGMBase::get_nodal_values(const MElement *e,const VectorStorageType &data)const
+{
+  std::vector<std::vector<double> > res(e->getNumVertices());
 
   for (int i=0;i<e->getNumVertices();i++){
     VectorStorageType::const_iterator itfind = data.find(const_cast<MVertex*>(e->getVertex(i)));
@@ -223,31 +217,30 @@ vector<vector<double> > BGMBase::get_nodal_values(const MElement *e,const Vector
   return res;
 }
 
-
-vector<double> BGMBase::get_nodal_values(const MElement *e,const DoubleStorageType &data)const{
-  if(debug) cout << "BGMBase::get_nodal_values(const MElement *e,const map<MVertex*,double> &data)" << endl;
-  vector<double> res(e->getNumVertices(),0.);
+std::vector<double> BGMBase::get_nodal_values(const MElement *e,const DoubleStorageType &data)const
+{
+  std::vector<double> res(e->getNumVertices(),0.);
 
   for (int i=0;i<e->getNumVertices();i++)
     res[i] = (data.find(const_cast<MVertex*>(e->getVertex(i))))->second;
   return res;
 }
 
-
-vector<double> BGMBase::get_element_uvw_from_xyz (const MElement *e, double x, double y,double z) const{
+std::vector<double> BGMBase::get_element_uvw_from_xyz (const MElement *e, double x, double y,double z) const
+{
   double element_uvw[3];
   double xyz[3] = {x, y, z};
   e->xyz2uvw(xyz, element_uvw);
-  vector<double> res(3,0.);
+  std::vector<double> res(3,0.);
   for (int i=0;i<3;i++) {
     res[i] = element_uvw[i];
   }
   return res;
 }
 
-
-set<MVertex*> BGMBase::get_vertices_of_maximum_dim(int dim){
-  set<MVertex*> bnd_vertices;
+std::set<MVertex*> BGMBase::get_vertices_of_maximum_dim(int dim)
+{
+  std::set<MVertex*> bnd_vertices;
   for(unsigned int i=0;i<gf->getNumMeshElements();i++){
     MElement* element = gf->getMeshElement(i);
     for(int j=0;j<element->getNumVertices();j++){
@@ -258,8 +251,8 @@ set<MVertex*> BGMBase::get_vertices_of_maximum_dim(int dim){
   return bnd_vertices;
 }
 
-
-GEntity* BGMBase::getBackgroundGEntity(){
+GEntity* BGMBase::getBackgroundGEntity()
+{
   return gf;
 }
 
diff --git a/Mesh/BGMBase.h b/Mesh/BGMBase.h
index 1ba8e99125..5cf5deadbb 100644
--- a/Mesh/BGMBase.h
+++ b/Mesh/BGMBase.h
@@ -30,8 +30,6 @@ class MElement;
 class MVertex;
 class GEntity;
 
-using namespace std;
-
 class BGMBase{
   public:
     typedef MVertex* hash_key_ptr;
@@ -40,7 +38,7 @@ class BGMBase{
 //    typedef tr1::unordered_map<hash_key_ptr, vector<double> > VectorStorageType;
     typedef std::map<hash_key_ptr, STensor3 > TensorStorageType;
     typedef std::map<hash_key_ptr, double > DoubleStorageType;
-    typedef std::map<hash_key_ptr, vector<double> > VectorStorageType;
+    typedef std::map<hash_key_ptr, std::vector<double> > VectorStorageType;
 
   protected:
     mutable MElementOctree *octree;
@@ -52,9 +50,9 @@ class BGMBase{
 
     int DIM,order;
 
-    virtual void export_scalar(const string &filename, const DoubleStorageType&) const;
-    virtual void export_vector(const string &filename, const VectorStorageType&) const;
-    virtual void export_tensor_as_vectors(const string &filename, const TensorStorageType &_whatToPrint)const;
+    virtual void export_scalar(const std::string &filename, const DoubleStorageType&) const;
+    virtual void export_vector(const std::string &filename, const VectorStorageType&) const;
+    virtual void export_tensor_as_vectors(const std::string &filename, const TensorStorageType &_whatToPrint)const;
     
     virtual void propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane = false)=0;
     virtual void computeSizeField()=0;
@@ -62,16 +60,16 @@ class BGMBase{
     virtual const MElement* getElement(unsigned int i)const=0;
     virtual unsigned int getNumMeshElements()const=0;
 
-    virtual vector<double> get_nodal_values(const MElement *e,const DoubleStorageType &data)const;
-    virtual vector<vector<double> > get_nodal_values(const MElement *e,const VectorStorageType &data)const;
+    virtual std::vector<double> get_nodal_values(const MElement *e,const DoubleStorageType &data)const;
+    virtual std::vector<std::vector<double> > get_nodal_values(const MElement *e,const VectorStorageType &data)const;
 
-    virtual vector<double> get_nodal_value(const MVertex *v,const VectorStorageType &data)const;
+    virtual std::vector<double> get_nodal_value(const MVertex *v,const VectorStorageType &data)const;
     virtual double get_nodal_value(const MVertex *v,const DoubleStorageType &data)const;
 
-    virtual vector<double> get_element_uvw_from_xyz (const MElement *e, double x, double y,double z) const;
+    virtual std::vector<double> get_element_uvw_from_xyz (const MElement *e, double x, double y,double z) const;
     
     virtual double get_field_value(double u, double v, double w, const DoubleStorageType &data);
-    virtual vector<double> get_field_value(double u, double v, double w, const VectorStorageType &data);
+    virtual std::vector<double> get_field_value(double u, double v, double w, const VectorStorageType &data);
 
   
   public:
@@ -80,8 +78,6 @@ class BGMBase{
     
     virtual MElementOctree* getOctree()=0;
 
-    const bool debug;
-
     virtual GEntity* getBackgroundGEntity();
 
     virtual double size(double u, double v, double w=0.);// get the size field, element interpolation
@@ -89,10 +85,10 @@ class BGMBase{
 
     virtual bool inDomain (double u, double v, double w=0.);
 
-    virtual inline void exportSizeField(const string &filename) const{export_scalar(filename,sizeField);};
+    virtual inline void exportSizeField(const std::string &filename) const{export_scalar(filename,sizeField);};
 
     // warning: these are "3D", "real" vertices, not copies in a parametric domain
-    virtual set<MVertex*> get_vertices_of_maximum_dim(int dim);
+    virtual std::set<MVertex*> get_vertices_of_maximum_dim(int dim);
     virtual const MElement* findElement(double u, double v, double w=0.,bool strict=true);
 
 };
diff --git a/Mesh/BackgroundMesh2D.cpp b/Mesh/BackgroundMesh2D.cpp
index 3170f88841..c8f3bb3ef0 100644
--- a/Mesh/BackgroundMesh2D.cpp
+++ b/Mesh/BackgroundMesh2D.cpp
@@ -31,10 +31,6 @@
 #include "linearSystemPETSc.h"
 #endif
 
-
-using namespace std;
-
-
 class evalDiffusivityFunction : public simpleFunction<double>{
   public:
     evalDiffusivityFunction(frameFieldBackgroundMesh2D *_bgm, double t=0.95):bgm(_bgm),threshold(t){};
@@ -50,16 +46,16 @@ class evalDiffusivityFunction : public simpleFunction<double>{
 //TODO: move this fct ???
 /* applies rotations of amplitude pi to set the 
    angle in the first quadrant (in [0,pi/2[ ) */
-void normalizeAngle(double &angle) {
+void normalizeAngle(double &angle) 
+{
   if (angle < 0) 
     while ( angle <  0 ) angle += (M_PI * .5);
   else if (angle >= M_PI * .5) 
     while ( angle >= M_PI * .5 ) angle -= (M_PI * .5);
 }
 
-
-void backgroundMesh2D::create_face_mesh(){
-
+void backgroundMesh2D::create_face_mesh()
+{
   quadsToTriangles(dynamic_cast<GFace*>(gf), 100000);
 
   // storing the initial mesh from GFace
@@ -83,7 +79,8 @@ void backgroundMesh2D::create_face_mesh(){
 }
 
 
-MElementOctree* backgroundMesh2D::getOctree(){
+MElementOctree* backgroundMesh2D::getOctree()
+{
   if(!octree){
     Msg::Debug("Rebuilding BackgroundMesh element octree");
     octree = new MElementOctree(elements);
@@ -91,13 +88,13 @@ MElementOctree* backgroundMesh2D::getOctree(){
   return octree;
 }
 
-
-const MElement* backgroundMesh2D::getElement(unsigned int i)const{
+const MElement* backgroundMesh2D::getElement(unsigned int i)const
+{
   return elements[i];
 }
 
-
-void backgroundMesh2D::reset(bool erase_2D3D){
+void backgroundMesh2D::reset(bool erase_2D3D)
+{
   unset();
 
   // create face mesh - this was previously done for old backgroundmesh in buildBackGroundMesh !
@@ -108,7 +105,7 @@ void backgroundMesh2D::reset(bool erase_2D3D){
     computeSizeField();
   }
   else 
-    for (map<MVertex*, MVertex*>::iterator itv2 = _2Dto3D.begin() ; itv2 != _2Dto3D.end(); ++itv2)
+    for (std::map<MVertex*, MVertex*>::iterator itv2 = _2Dto3D.begin() ; itv2 != _2Dto3D.end(); ++itv2)
       sizeField[itv2->first] = CTX::instance()->mesh.lcMax;
 
   // ensure that other criteria are fullfilled
@@ -131,14 +128,14 @@ void backgroundMesh2D::unset(){
 
 void backgroundMesh2D::create_mesh_copy(){
   // TODO: useful to extend it to other elements ???
-  //set<SPoint2> myBCNodes;
+  //std::set<SPoint2> myBCNodes;
   GFace *face = dynamic_cast<GFace*>(gf);
   for (unsigned int i = 0; i < face->triangles.size(); i++){
     MTriangle *e = face->triangles[i];
     MVertex *news[3];
     for (int j=0;j<3;j++){
       MVertex *v = e->getVertex(j);
-      map<MVertex*,MVertex*>::iterator it = _3Dto2D.find(v);
+      std::map<MVertex*,MVertex*>::iterator it = _3Dto2D.find(v);
       MVertex *newv =0;
       if (it == _3Dto2D.end()){
         SPoint2 p;
@@ -162,20 +159,20 @@ GPoint backgroundMesh2D::get_GPoint_from_MVertex(const MVertex *v)const{
 }
 
 
-backgroundMesh2D::backgroundMesh2D(GFace *_gf, bool erase_2D3D):BGMBase(2,_gf),sizeFactor(1.){
-  if(debug) cout << "backgroundMesh2D::constructor " << endl;
+backgroundMesh2D::backgroundMesh2D(GFace *_gf, bool erase_2D3D):BGMBase(2,_gf),sizeFactor(1.)
+{
   reset(erase_2D3D);
 
   if (erase_2D3D){
-    // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh from GFace, back to the previous one !
+    // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh 
+    // from GFace, back to the previous one ! 
     GFace *face = dynamic_cast<GFace*>(gf);
     face->triangles = tempTR;
   }
-
 }
 
-
-backgroundMesh2D::~backgroundMesh2D(){
+backgroundMesh2D::~backgroundMesh2D()
+{
   unset();
 }
 
@@ -203,9 +200,8 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
     myAssembler.fixVertex(itv->first, 0, 1, itv->second);
   }
 
-
   // Number vertices
-  set<MVertex*> vs;
+  std::set<MVertex*> vs;
   GFace *face = dynamic_cast<GFace*>(gf);
   for (unsigned int k = 0; k < face->triangles.size(); k++)
     for (int j=0;j<3;j++)vs.insert(face->triangles[k]->getVertex(j));
@@ -213,9 +209,9 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
     for (int j=0;j<4;j++)vs.insert(face->quadrangles[k]->getVertex(j));
 
 
-  map<MVertex*,SPoint3> theMap;
+  std::map<MVertex*,SPoint3> theMap;
   if ( in_parametric_plane) {
-    for (set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+    for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
       SPoint2 p;
       reparamMeshVertexOnFace ( *it, face, p);
       theMap[*it] = SPoint3((*it)->x(),(*it)->y(),(*it)->z());
@@ -223,7 +219,7 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
     }
   }
 
-  for (set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it)
+  for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it)
     myAssembler.numberVertex(*it, 0, 1);
 
   // Assemble
@@ -240,12 +236,12 @@ void backgroundMesh2D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
   }
 
   // save solution
-  for (set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+  for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
     myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
   }
 
   if ( in_parametric_plane) {
-    for (set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+    for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
       SPoint3 p = theMap[(*it)];
       (*it)->setXYZ(p.x(),p.y(),p.z());
     }
@@ -287,7 +283,7 @@ void backgroundMesh2D::computeSizeField(){
   simpleFunction<double> ONE(1.0);
   propagateValues(sizes,ONE);
 
-  map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+  std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
   for ( ; itv2 != _2Dto3D.end(); ++itv2){
     MVertex *v_2D = itv2->first;
     MVertex *v_3D = itv2->second;
@@ -330,7 +326,7 @@ void backgroundMesh2D::updateSizes(){
   // do not allow large variations in the size field
   // (Int. J. Numer. Meth. Engng. 43, 1143-1165 (1998) MESH GRADATION
   // CONTROL, BOROUCHAKI, HECHT, FREY)
-  set<MEdge,Less_Edge> edges;
+  std::set<MEdge,Less_Edge> edges;
   for (unsigned int i = 0; i < getNumMeshElements(); i++){
     for (int j = 0; j < getElement(i)->getNumEdges(); j++){
       edges.insert(getElement(i)->getEdge(j));
@@ -338,7 +334,7 @@ void backgroundMesh2D::updateSizes(){
   }
   const double _beta = 1.3;
   for (int i=0;i<0;i++){
-    set<MEdge,Less_Edge>::iterator it = edges.begin();
+    std::set<MEdge,Less_Edge>::iterator it = edges.begin();
     for ( ; it != edges.end(); ++it){
       MVertex *v0 = it->getVertex(0);
       MVertex *v1 = it->getVertex(1);
@@ -356,27 +352,21 @@ void backgroundMesh2D::updateSizes(){
 
 
 
-frameFieldBackgroundMesh2D::frameFieldBackgroundMesh2D(GFace *_gf):backgroundMesh2D(_gf,false){
-
-  if(debug) cout << "frameFieldBackgroundMesh2D::constructor " << endl;
-
+frameFieldBackgroundMesh2D::frameFieldBackgroundMesh2D(GFace *_gf):backgroundMesh2D(_gf,false)
+{
   reset();
 
-  // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh from GFace, back to the previous one !
+  // now, the new mesh has been copied in local in backgroundMesh2D, deleting the mesh
+  // from GFace, back to the previous one !
   GFace *face = dynamic_cast<GFace*>(gf);
   face->triangles = tempTR;
 }
 
-
 frameFieldBackgroundMesh2D::~frameFieldBackgroundMesh2D(){}
 
-
-void frameFieldBackgroundMesh2D::reset(bool erase_2D3D){
-
+void frameFieldBackgroundMesh2D::reset(bool erase_2D3D)
+{
   // computes cross field
-  if(debug) cout << "frameFieldBackgroundMesh2D::reset " << endl;
-
-
   simpleFunction<double> ONE(1.0);
   computeCrossField(ONE);
   computeSmoothness();
@@ -402,46 +392,42 @@ void frameFieldBackgroundMesh2D::reset(bool erase_2D3D){
   }
 }
 
-
-double frameFieldBackgroundMesh2D::get_smoothness(MVertex *v){
-  if(debug) cout << "frameFieldBackgroundMesh2D::get_smoothness(MVertex *v)" << endl;
+double frameFieldBackgroundMesh2D::get_smoothness(MVertex *v)
+{
   return get_nodal_value(v,smoothness);
 }
 
-
-double frameFieldBackgroundMesh2D::get_smoothness(double u, double v){
-  if(debug) cout << "frameFieldBackgroundMesh2D::get_smoothness(double u, double v)" << endl;
+double frameFieldBackgroundMesh2D::get_smoothness(double u, double v)
+{
   return get_field_value(u,v,0.,smoothness);
 }
 
-
-double frameFieldBackgroundMesh2D::angle(MVertex *v){
-  if(debug) cout << "frameFieldBackgroundMesh2D::angle(MVertex *v)" << endl;
+double frameFieldBackgroundMesh2D::angle(MVertex *v)
+{
   return get_nodal_value(v,angles);
 }
 
-
-double frameFieldBackgroundMesh2D::angle(double u, double v){
+double frameFieldBackgroundMesh2D::angle(double u, double v)
+{
   MElement *e = const_cast<MElement*>(findElement(u, v));
   if (!e) return -1000.0;
-  vector<double> val = get_nodal_values(e,angles);
-  vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,0.);
-  double cosvalues[e->getNumVertices()],sinvalues[e->getNumVertices()];
+  std::vector<double> val = get_nodal_values(e,angles);
+  std::vector<double> element_uvw = get_element_uvw_from_xyz(e,u,v,0.);
+  std::vector<double> cosvalues(e->getNumVertices()), sinvalues(e->getNumVertices());
   for (int i=0;i<e->getNumVertices();i++){
     cosvalues[i]=cos(4*val[i]);
     sinvalues[i]=sin(4*val[i]);
   }
-  double cos4 = e->interpolate(cosvalues, element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
-  double sin4 = e->interpolate(sinvalues, element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
+  double cos4 = e->interpolate(&cosvalues[0], element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
+  double sin4 = e->interpolate(&sinvalues[0], element_uvw[0], element_uvw[1], element_uvw[2], 1, order);
   double a = atan2(sin4,cos4)/4.0;
   normalizeAngle (a);
   return a;
 }
 
-
-void frameFieldBackgroundMesh2D::computeCrossField(simpleFunction<double> &eval_diffusivity){
+void frameFieldBackgroundMesh2D::computeCrossField(simpleFunction<double> &eval_diffusivity)
+{
   angles.clear();
-  if(debug) cout << "frameFieldBackgroundMesh2D::computeCrossField " << endl;
 
   DoubleStorageType _cosines4,_sines4;
 
@@ -488,9 +474,7 @@ void frameFieldBackgroundMesh2D::computeCrossField(simpleFunction<double> &eval_
   propagateValues(_cosines4,eval_diffusivity,false);
   propagateValues(_sines4,eval_diffusivity,false);
 
-
-
-  map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+  std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
   for ( ; itv2 != _2Dto3D.end(); ++itv2){
     MVertex *v_2D = itv2->first;
     MVertex *v_3D = itv2->second;
@@ -500,8 +484,8 @@ void frameFieldBackgroundMesh2D::computeCrossField(simpleFunction<double> &eval_
   }
 }
 
-
-void frameFieldBackgroundMesh2D::eval_crossfield(double u, double v, STensor3 &cf){
+void frameFieldBackgroundMesh2D::eval_crossfield(double u, double v, STensor3 &cf)
+{
   double quadAngle  = angle(u,v);
   Pair<SVector3, SVector3> dirs = compute_crossfield_directions(u,v,quadAngle);
   SVector3 n = crossprod(dirs.first(),dirs.second());
@@ -512,8 +496,6 @@ void frameFieldBackgroundMesh2D::eval_crossfield(double u, double v, STensor3 &c
     cf(i,2) = n[i];
   }
 
-
-
   //  SVector3 t1,t2,n;
   //  GFace *face = dynamic_cast<GFace*>(gf);
   //  Pair<SVector3, SVector3> der = face->firstDer(SPoint2(u,v));
@@ -531,26 +513,22 @@ void frameFieldBackgroundMesh2D::eval_crossfield(double u, double v, STensor3 &c
   //    cf(i,1) = t2[i];
   //    cf(i,2) = n[i];
   //  }
-
-  return;
 }
 
-
-void frameFieldBackgroundMesh2D::eval_crossfield(MVertex *vert, STensor3 &cf){
+void frameFieldBackgroundMesh2D::eval_crossfield(MVertex *vert, STensor3 &cf)
+{
   SPoint2 parampoint;
   GFace *face = dynamic_cast<GFace*>(gf);
   reparamMeshVertexOnFace(vert, face, parampoint);
   return eval_crossfield(parampoint[0], parampoint[1], cf);
 }
 
-
-void frameFieldBackgroundMesh2D::computeSmoothness(){
+void frameFieldBackgroundMesh2D::computeSmoothness()
+{
   smoothness.clear();
 
-  if(debug) cout << "frameFieldBackgroundMesh2D::computeSmoothness " << endl;
-
   // build vertex -> neighbors table
-  multimap<MVertex*,MVertex*> vertex2vertex;
+  std::multimap<MVertex*,MVertex*> vertex2vertex;
   for (std::vector<MElement*>::iterator it = beginelements();it!=endelements();it++){
     MElement *e = *it;
     for (int i=0;i<e->getNumVertices();i++){
@@ -568,10 +546,10 @@ void frameFieldBackgroundMesh2D::computeSmoothness(){
     MVertex *v = *it;
     double angle_current = angle(v);
     // compare to all neighbors...
-    pair<multimap<MVertex*,MVertex*>::iterator, multimap<MVertex*,MVertex*>::iterator> range = vertex2vertex.equal_range(v);
+    std::pair<std::multimap<MVertex*,MVertex*>::iterator, std::multimap<MVertex*,MVertex*>::iterator> range = vertex2vertex.equal_range(v);
     double minangle,totalangle=0.;
     int N=0;
-    for (multimap<MVertex*,MVertex*>::iterator itneighbor = range.first;itneighbor!=range.second;itneighbor++){
+    for (std::multimap<MVertex*,MVertex*>::iterator itneighbor = range.first;itneighbor!=range.second;itneighbor++){
       N++;
       minangle=M_PI/2;
       MVertex *v_nb = itneighbor->second;
@@ -587,11 +565,11 @@ void frameFieldBackgroundMesh2D::computeSmoothness(){
   }
 }
 
-
-void frameFieldBackgroundMesh2D::exportCrossField(const string &filename){
+void frameFieldBackgroundMesh2D::exportCrossField(const std::string &filename)
+{
   FILE *f = Fopen (filename.c_str(),"w");
   fprintf(f,"View \"Cross Field\"{\n");
-  vector<double> deltas(2);
+  std::vector<double> deltas(2);
   deltas[0] = 0.;
   deltas[1] = M_PI;
 
@@ -609,10 +587,9 @@ void frameFieldBackgroundMesh2D::exportCrossField(const string &filename){
   fclose(f);
 }
 
-
 // returns the cross field as a pair of othogonal vectors (NOT in parametric coordinates, but real 3D coordinates)
-Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directions(double u,double v, double angle_current){
-
+Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directions(double u,double v, double angle_current)
+{
   // get the unit normal at that point
   GFace *face = dynamic_cast<GFace*>(gf);
   Pair<SVector3, SVector3> der = face->firstDer(SPoint2(u,v));
@@ -636,8 +613,8 @@ Pair<SVector3, SVector3> frameFieldBackgroundMesh2D::compute_crossfield_directio
       SVector3(t2[0],t2[1],t2[2]));
 }
 
-
-bool frameFieldBackgroundMesh2D::compute_RK_infos(double u,double v, double x, double y, double z, RK_form &infos){
+bool frameFieldBackgroundMesh2D::compute_RK_infos(double u,double v, double x, double y, double z, RK_form &infos)
+{
 
   // check if point is in domain
   if (!inDomain(u,v)) return false;
@@ -746,33 +723,4 @@ bool frameFieldBackgroundMesh2D::compute_RK_infos(double u,double v, double x, d
   infos.normal = n;
 
   return true;
-
 }
-
-
-
-
-//  if (use_previous_basis){
-//    std::map<double, double> angles;
-//    SVector3 temp = crossprod(previous_t1, t1);
-//    double a = atan2(dot(t1, previous_t1), sign(dot(temp,n))*temp.norm() );
-//    angles.insert(std::make_pair(abs(a),a));
-//    temp = crossprod(previous_t2, t1);
-//    a = atan2(dot(t1, previous_t2), sign(dot(temp,n))*temp.norm());
-//    angles.insert(std::make_pair(abs(a),a));
-//    temp = crossprod(-1.*previous_t1, t1);
-//    a = atan2(dot(t1, -1.*previous_t1), sign(dot(temp,n))*temp.norm());
-//    angles.insert(std::make_pair(abs(a),a));
-//    temp = crossprod(-1.*previous_t2, t1);
-//    a = atan2(dot(t1, -1.*previous_t2), sign(dot(temp,n))*temp.norm());
-//    angles.insert(std::make_pair(abs(a),a));
-//    //    std::if(debug) cout << "angles: " << std::endl;
-//    //    for (int i=0;i<4;i++)  std::if(debug) cout << angles[i] << "  " << std::endl;
-//    double min_angle = -(angles.begin()->second); 
-//    //    std::if(debug) cout << "min angle = " << min_angle << std::endl;
-//    t1 =  cos(min_angle)*previous_t1 + sin(min_angle)*previous_t2;
-//    t2 = -sin(min_angle)*previous_t1 + cos(min_angle)*previous_t2;
-//    //    std::if(debug) cout << "new corrected t        : (" << t1(0) << "," <<  t1(1) << ") (" << t2(0) << "," << t2(1) << std::endl;
-//  }
-
-
diff --git a/Mesh/pointInsertion.cpp b/Mesh/pointInsertion.cpp
index 6060a240e8..9fbcab3656 100644
--- a/Mesh/pointInsertion.cpp
+++ b/Mesh/pointInsertion.cpp
@@ -1,3 +1,7 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
 
 #include "pointInsertion.h"
 #include "BackgroundMeshManager.h"
@@ -5,35 +9,23 @@
 #include "BackgroundMesh3D.h"
 #include "GFace.h"
 #include "GRegion.h"
-
-
 #include "OS.h"
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <string>
-
 #include "Context.h"
 #include "meshGRegion.h"
-
 #include "pointInsertionRTreeTools.h"
-
 #include "intersectCurveSurface.h"
 
 //#include "google/profiler.h"
 
-
 using namespace std;
 
-
 bool old_algo_hexa(){
   return true;
 }
 
-
 template<typename T>
-void print_nodal_info(string filename, map<MVertex*, T> &mapp){
+void print_nodal_info(string filename, map<MVertex*, T> &mapp)
+{
   ofstream out(filename.c_str());
 
   out << "View \"\"{" << endl;
@@ -283,21 +275,19 @@ void Filler2D::pointInsertion2D(GFace* gf,  vector<MVertex*> &packed, vector<SMe
 
   if (debug) cout << "ENTERING POINTINSERTION2D" << endl;
   
-  clock_t a;
+  double a;
 
   // acquire background mesh
   if(debug) cout << "pointInsertion2D: recover BGM" << endl;
-  a=clock();
+  a=Cpu();
   frameFieldBackgroundMesh2D *bgm = dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(gf));
-  time_bgm_and_smoothing += (clock() - a) / (double)CLOCKS_PER_SEC;
-
+  time_bgm_and_smoothing += (Cpu() - a);
 
   if (!bgm){
-    cout << "pointInsertion2D:: BGM dynamic cast failed ! " << endl;
-    throw;
+    Msg::Error("BGM dynamic cast failed in filler2D::pointInsertion2D");
+    return;
   }
 
-
   // export BGM size field
   if(export_stuff){
     cout << "pointInsertion2D: export size field " << endl;
@@ -319,7 +309,7 @@ void Filler2D::pointInsertion2D(GFace* gf,  vector<MVertex*> &packed, vector<SMe
 
 
   // point insertion algorithm:
-  a=clock();
+  a=Cpu();
 
   // for debug check...
   int priority_counter=0;
@@ -417,7 +407,7 @@ void Filler2D::pointInsertion2D(GFace* gf,  vector<MVertex*> &packed, vector<SMe
     }
     if(debug) cout << "////////// nbre of added point: " << count_nbaddedpt << endl;
   }
-  time_insertion += (clock() - a) / (double)CLOCKS_PER_SEC;
+  time_insertion += (Cpu() - a);
 
 
   if (debug){
@@ -496,16 +486,16 @@ bool Filler3D::treat_region(GRegion *gr){
 
   const bool debug=false;
   const bool export_stuff=true;
-  clock_t a;
+  double a;
 
   cout << "ENTERING POINTINSERTION3D" << endl;
 
 
   // acquire background mesh
   cout << "pointInsertion3D: recover BGM" << endl;
-  a = clock();
+  a = Cpu();
   frameFieldBackgroundMesh3D *bgm = dynamic_cast<frameFieldBackgroundMesh3D*>(BGMManager::get(gr));
-  time_smoothing += (clock() - a) / (double)CLOCKS_PER_SEC;
+  time_smoothing += (Cpu() - a);
 
   if (!bgm){
     cout << "pointInsertion3D:: BGM dynamic cast failed ! " << endl;
@@ -544,7 +534,7 @@ bool Filler3D::treat_region(GRegion *gr){
   cout << "pointInsertion3D : inserting points in region " << gr->tag()  << endl;
 
   //ProfilerStart("/home/bernard/profile");
-  a = clock();
+  a = Cpu();
 
   // ----- initialize fifo list -----
 
@@ -741,7 +731,7 @@ bool Filler3D::treat_region(GRegion *gr){
 
   //ProfilerStop();
 
-  time_insert_points += (clock() - a) / (double)CLOCKS_PER_SEC;
+  time_insert_points += (Cpu() - a);
 
 
   // --- output ---
@@ -760,7 +750,7 @@ bool Filler3D::treat_region(GRegion *gr){
   // ------- meshing using new points
   cout << "tets in gr before= " << gr->tetrahedra.size() << endl;
   cout << "nb new vertices= " << new_vertices.size() << endl;
-  a=clock();
+  a=Cpu();
 
   int option = CTX::instance()->mesh.algo3d;
   CTX::instance()->mesh.algo3d = ALGO_3D_DELAUNAY;
@@ -773,7 +763,7 @@ bool Filler3D::treat_region(GRegion *gr){
   meshGRegion mesher(regions); //?
   mesher(gr); //?
   MeshDelaunayVolume(regions);
-  time_meshing += (clock() - a) / (double)CLOCKS_PER_SEC;
+  time_meshing += (Cpu() - a);
 
   cout << "tets in gr after= " << gr->tetrahedra.size() << endl;
   cout << "gr tag=" << gr->tag() << endl;
@@ -812,21 +802,8 @@ Filler3D::~Filler3D(){
   cout << "  ------- CUMULATIVE TOTAL 3D TIME (new)   : " << time_meshing+time_smoothing+time_insert_points << " s." << endl;
 }
 
-
 std::vector<MVertex*> Filler3D::new_vertices;
 
-
 double Filler3D::time_smoothing = 0.;
 double Filler3D::time_insert_points = 0.;
 double Filler3D::time_meshing = 0.;
-
-
-
-
-
-
-
-
-
-
-
diff --git a/Mesh/pointInsertion.h b/Mesh/pointInsertion.h
index 74f3c7557e..ecda4439e7 100644
--- a/Mesh/pointInsertion.h
+++ b/Mesh/pointInsertion.h
@@ -7,8 +7,6 @@
 #ifndef _POINTINSERTION_H_
 #define _POINTINSERTION_H_
 
-
-
 #include "STensor3.h"
 #include <vector>
 #include <string>
@@ -19,61 +17,27 @@ class GFace;
 class GRegion;
 class MVertex;
 
-using namespace std;
-
 extern bool old_algo_hexa();
 
-
 class Filler2D{
-  public:
-    Filler2D();
-    ~Filler2D();
-    void pointInsertion2D(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
-  private:
-    static double time_bgm_and_smoothing,time_insertion;
+ public:
+  Filler2D();
+  ~Filler2D();
+  void pointInsertion2D(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
+ private:
+  static double time_bgm_and_smoothing,time_insertion;
 };
 
-
 class Filler3D{
-  private:
-    static std::vector<MVertex*> new_vertices;// these are used in meshGRegion.cpp using static !!!
-    static double time_smoothing,time_insert_points, time_meshing;
-  public:
-    Filler3D();
-    ~Filler3D();
-    virtual bool treat_region(GRegion*);
-    static int get_nbr_new_vertices();
-    static MVertex* get_new_vertex(int);
+ private:
+  static std::vector<MVertex*> new_vertices;// these are used in meshGRegion.cpp using static !!!
+  static double time_smoothing,time_insert_points, time_meshing;
+ public:
+  Filler3D();
+  ~Filler3D();
+  virtual bool treat_region(GRegion*);
+  static int get_nbr_new_vertices();
+  static MVertex* get_new_vertex(int);
 };
 
-
-
-
-
-
-//template<class T>
-//bool readValue(string filename, string keystr, T &value){
-//  ifstream in(filename.c_str());
-//  if (!in.is_open()){
-//    cerr << "Can't open file " << filename << endl;
-//    throw;
-//    return false;
-//  }
-//  string name;
-//
-//  bool found=false;
-//  while (!(in.eof())){
-//    in.ignore(100,'$');
-//    in >> name;
-//    if (in.eof()) break;
-//    if (!name.compare(keystr.c_str())){// strings are equals
-//      found=true;
-//      in >> value;
-//    }
-//  }
-//  in.close();
-//  return found;
-//}
-
-
 #endif
-- 
GitLab