Curvature.cpp

// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#include "Curvature.h"
#include "MElement.h"
#include "MTriangle.h"
#include "GEntity.h"
#include "GFaceCompound.h"
#include "MLine.h"
#include "GRbf.h"
#include "OS.h"
#include "SBoundingBox3d.h"
#include "discreteEdge.h"

#include<iostream>
#include<fstream>
#include<cmath>

#define NEXT(i) ((i)<2 ? (i)+1 : (i)-2)
#define PREV(i) ((i)>0 ? (i)-1 : (i)+2)
//========================================================================================================

//Initialization of the static variables:
Curvature* Curvature::_instance = 0;
bool Curvature::_destroyed = false;
bool Curvature::_alreadyComputedCurvature = false;

//========================================================================================================

//CONSTRUCTOR
Curvature::Curvature() : _isMapInitialized(false)
{
}

// Curvature::Curvature(std::list<GFace*> &myFaces){

//   std::list<GFace*>::const_iterator it = myFaces.begin();
//   for( ; it != myFaces.end() ; ++it){
//          _ptFinalEntityList.push_back(*it);
//   }

// }

//========================================================================================================

Curvature::~Curvature()
{
   _instance = 0;
  _destroyed = true;

}
//========================================================================================================
void Curvature::onDeadReference()
{
  std::cout << "Dead reference of Curvature detected" << std::endl;
}
//========================================================================================================

Curvature& Curvature::getInstance()
{
  if (!_instance)  {
    if(_destroyed){
      onDeadReference();
    }
    else{
      create();
    }
  }
  return *_instance;
}
//========================================================================================================
  bool Curvature::valueAlreadyComputed()
  {
    return _alreadyComputedCurvature;
  }

//========================================================================================================

 void Curvature::create()
 {
   static Curvature instance;
   _instance = & instance;
 }

 //========================================================================================================
 void Curvature::retrieveCompounds()  {
#if defined(HAVE_SOLVER)

   /// -------------------------------------------------------------------------------------------
   // Loop over all faces. Check if the face is a compound. If it is, store all of its discrete
   // faces in _EntityArray

   std::list<GFace*> global_face_list;

   for(GModel::fiter face_iter = _model->firstFace(); face_iter != _model->lastFace(); ++face_iter)
   {
     if ( (*face_iter)->geomType() == GEntity::CompoundSurface )
     {

       GFaceCompound* compound = dynamic_cast<GFaceCompound*>(*face_iter);
       std::list<GFace*> tempcompounds = compound->getCompounds();
       std::list<GFace*>::iterator surfIterator;

       for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator)
       {
         if ((*surfIterator)->geomType() == GEntity::DiscreteSurface)
         {
           global_face_list.push_back(*surfIterator);
         }
       }
     }

   } // Loop over faces of the model

   global_face_list.sort();
   global_face_list.unique();
   _EntityArray.resize(global_face_list.size());
   std::copy(global_face_list.begin(),global_face_list.end(),_EntityArray.begin());

//   std::cout << "Retrieve compounds: some preliminary info:" << std::endl;
//   std::cout << "The size of _ptFinalEntityList = " << _EntityArray.size() << std::endl;
//   std::cout << "Number of elements in all faces: " << std::endl;
//   for(int i = 0; i < _EntityArray.size(); ++i)
//   {
//     std::cout << "\t" << _EntityArray[i]->getNumMeshElements() << " elements" << std::endl;
//   }

#endif
 }

 //========================================================================================================

 // Method to detect boundary edges of the mesh. We want to find which edges are on the boundary in order
 // to correct the curvature field close to boundaries

 void Curvature::correctOnEdges()
 {
#if defined(HAVE_SOLVER)

   if (! _alreadyComputedCurvature )
   {
     Msg::Error("Cannot correct the curvature values at the edges because the curvatures weren't computed yet");
     return;
   }
   // Remark: this can only be used after the call to initializeMap() !

  buildEdgeList();

  std::list<MeshEdgeInfo>::iterator VertToEdgeListIter;

   _isOnBoundary.resize(_VertexToInt.size());
   for (int i = 0; i < _VertexToInt.size(); ++i)
   {
     _isOnBoundary[i] = 0;
   }

   // To detect the nodes on the egdes of a geometry, we create a list of all edges on the mesh. The
   // edges which are shared by only one mesh element are boundary edges. Their nodes are tagged by 1
   for (int i = 0; i < _VertexToInt.size(); ++i)
   {
     for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
                                                                                    ++VertToEdgeListIter)
     {
       if ((*VertToEdgeListIter).NbElemNeighbour == 1)
       {
         _isOnBoundary[(*VertToEdgeListIter).StartV] = 1;
         _isOnBoundary[(*VertToEdgeListIter).EndV] = 1;
       }
     }

   }

   // We want to find the nodes that are the immediate neighbours of the boundary nodes. Those nodes are
   // considered nodes with 'level 2'. The neighbours of neighbours have 'level 3' etc.
   // We want to construct levels 1,2,3. Nodes with level 0 are internal nodes of the mesh
   //Loop over the vector of edgelist
   for (int level = 1; level < 3; ++level)
   {
     for (int i = 0; i < _VertexToEdgeList.size(); ++i)
     {
       for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
                                                                                ++VertToEdgeListIter)
       {
         if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level && _isOnBoundary[(*VertToEdgeListIter).EndV] == 0)
         {
          _isOnBoundary[(*VertToEdgeListIter).EndV] = level+1;
         }
         if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level && _isOnBoundary[(*VertToEdgeListIter).StartV] == 0)
         {
           _isOnBoundary[(*VertToEdgeListIter).StartV] = level+1;
         }
       }
     }
   }//Loop over the level of the ring


   // Check test to see the tag of the nodes on the boundary:
//   for (int i = 0; i< _EntityArray.size(); ++i)
//   {
//     GFace* face = _EntityArray[i]; //List of all the discrete face on which we compute the curvature

//     for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
//     {
//       MElement *e = face->getMeshElement(iElem);

//       // The NEW tag of the corresponding element
//       const int E = _ElementToInt[e->getNum()];

//       // Pointers to vertices of triangle
//       MVertex* A = e->getVertex(0);
//       MVertex* B = e->getVertex(1);
//       MVertex* C = e->getVertex(2);

//       V[0] = _VertexToInt[A->getNum()];  //The new number of the vertex
//       V[1] = _VertexToInt[B->getNum()];
//       V[2] = _VertexToInt[C->getNum()];

//       if (_isOnBoundary[V[0]] == 3)
//       {
//         std::cout << "Vertex: " << A->getNum() << " -- " << A->x() << "; " << A->y() << "; " << A->z() << std::endl;
//       }
//       if (_isOnBoundary[V[1]] == 3)
//       {
//         std::cout << "Vertex: " << B->getNum() << " -- " << B->x() << "; " << B->y() << "; " << B->z() << std::endl;
//       }
//       if (_isOnBoundary[V[2]] == 3)
//       {
//         std::cout << "Vertex: " << C->getNum() << " -- " << C->x() << "; " << C->y() << "; " << C->z() << std::endl;
//       }
//     }
//   }


   // Now we'll propagate the cuvature values from inside nodes with level = 3 close to the boundary - first
   // to nodes with level 2, then from nodes with level 2 to nodes with level 1 (on the boundary)
   _NbNeighbour.resize(_VertexToInt.size());
   for(int i = 0; i < _NbNeighbour.size(); ++i)
   {
     _NbNeighbour[i] = 0;
   }

   for (int level = 2; level > 0 ; --level)
   {

     for (int i = 0; i< _NbNeighbour.size(); ++i)
     {
       _NbNeighbour[i] = 0;
       if (_isOnBoundary[i] == level)
       {
        _VertexCurve[i] = 0.0;
       }
     }


     for (int i = 0; i < _VertexToEdgeList.size(); ++i)
     {
       for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
                                                                                ++VertToEdgeListIter)
       {

         if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level && _isOnBoundary[(*VertToEdgeListIter).EndV] == level+1)
         {
            _VertexCurve[(*VertToEdgeListIter).StartV] += _VertexCurve[(*VertToEdgeListIter).EndV];
            _NbNeighbour[(*VertToEdgeListIter).StartV] ++;

         }
         if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level && _isOnBoundary[(*VertToEdgeListIter).StartV] == level+1)
         {
            _VertexCurve[(*VertToEdgeListIter).EndV] += _VertexCurve[(*VertToEdgeListIter).StartV];
            _NbNeighbour[(*VertToEdgeListIter).EndV] ++;
         }

       }
     }

     // Correction for a degenerate case when a node has neighbours with the same or lower level, but zero
     // neighbours with level+1

     for (int i = 0; i < _VertexToEdgeList.size(); ++i)
     {
       for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
                                                                                ++VertToEdgeListIter)
       {

         if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level
             && _isOnBoundary[(*VertToEdgeListIter).EndV] == level
             && _NbNeighbour[(*VertToEdgeListIter).StartV] == 0)
         {
           _VertexCurve[(*VertToEdgeListIter).StartV] += _VertexCurve[(*VertToEdgeListIter).EndV];
           _NbNeighbour[(*VertToEdgeListIter).StartV] = _NbNeighbour[(*VertToEdgeListIter).EndV];
         }
         if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level
             && _isOnBoundary[(*VertToEdgeListIter).StartV] == level
             && _NbNeighbour[(*VertToEdgeListIter).EndV] == 0)
         {
           _VertexCurve[(*VertToEdgeListIter).EndV] += _VertexCurve[(*VertToEdgeListIter).StartV];
           _NbNeighbour[(*VertToEdgeListIter).EndV] = _NbNeighbour[(*VertToEdgeListIter).StartV];
         }
       }
     }

     for (int i = 0; i < _isOnBoundary.size(); ++i)
     {
       if (_isOnBoundary[i] == level)
       {
         if (_NbNeighbour[i] == 0)
         {
           std::cout << "ERROR: VERTEX " << i << " WITH LEVEL " << level << " HAS 0 NEIGHBOURS" << std::endl;
         }
         _VertexCurve[i] = _VertexCurve[i]/_NbNeighbour[i];
       }
     }

   }//Loop over the levels of the ring

//   for (int i = 0; i < _isOnBoundary.size(); ++i)
//   {
//     _VertexCurve[i] = _isOnBoundary[i];
//   }

#endif
 }

//========================================================================================================

//INITIALIZATION OF THE MAP  AND  RENUMBERING OF THE SELECTED ENTITIES:

void Curvature::initializeMap()
{
  if (! _isMapInitialized)
  {
    for (int i = 0; i< _EntityArray.size(); ++i)
    {
      GFace* face = _EntityArray[i];

      for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
      {
        MElement *e = face->getMeshElement(iElem);
        const int E = e->getNum();
        _ElementToInt[E] = 1;

        const int A = e->getVertex(0)->getNum();   //Pointer to 1st vertex of the triangle
        const int B = e->getVertex(1)->getNum();
        const int C = e->getVertex(2)->getNum();

        _VertexToInt[A] = 1;
        _VertexToInt[B] = 1;
        _VertexToInt[C] = 1;
      }
    }

    /// Set up a new numbering of chosen vertices and triangles
    int idx = 0;

    // map between the pointer to vertex and the new numbering of the vertex
    std::map<int,int>::iterator vertex_iterator;
    // map between the pointer to element and the new numbering of the element
    std::map<int,int>::iterator element_iterator;

    for (vertex_iterator = _VertexToInt.begin() ; vertex_iterator !=_VertexToInt.end() ; ++ vertex_iterator, ++idx)
      (*vertex_iterator).second = idx;

    idx = 0;

    for (element_iterator = _ElementToInt.begin() ; element_iterator !=_ElementToInt.end() ; ++ element_iterator, ++idx)
      (*element_iterator).second = idx;

   _isMapInitialized = true;

  }

}

//========================================================================================================

//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND

void Curvature::computeVertexNormals()
{
  SVector3 vector_AB;
  SVector3 vector_AC;

  _TriangleArea.resize(_ElementToInt.size() );
  _VertexArea.resize(_ElementToInt.size() );
  _VertexNormal.resize(_VertexToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    // face is a pointer to one surface of the group "FinalEntityList"
    GFace* face = _EntityArray[i];

    //Loop over the element all the element of the "myTag"-surface
    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      // Pointer to one element
      MElement *e = face->getMeshElement(iElem);
      // The NEW tag of the corresponding element
      const int E = _ElementToInt[e->getNum()];
      // std::cout << "We are now looking at element Nr: " << E << std::endl;

      // Pointers to vertices of triangle
      MVertex* A = e->getVertex(0);
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int V0 = _VertexToInt[A->getNum()];  //The new number of the vertex
      const int V1 = _VertexToInt[B->getNum()];
      const int V2 = _VertexToInt[C->getNum()];

      vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
      vector_AC = SVector3(C->x() - A->x(), C->y() - A->y(), C->z() - A->z() );

      const SVector3 cross = crossprod(vector_AB, vector_AC);

      // Area of the triangles:
      _TriangleArea[E] = 0.5*cross.norm();
      // std::cout << "The area of the triangle nr: " << e->getNum() << " is: "<< TriangleArea[E] << std::endl;

      _VertexArea[V0] += _TriangleArea[E];
      _VertexArea[V1] += _TriangleArea[E];
      _VertexArea[V2] += _TriangleArea[E];

      _VertexNormal[V0] += cross;  //here we are actually computing the unit normal vector per vertex
      _VertexNormal[V1] += cross;
      _VertexNormal[V2] += cross;

    } // end of loop over elements of one face

  } //Loop over _ptFinalEntityList

    ///////Normalize the vertex-normals.
    for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
    {
      _VertexNormal[n].normalize();
    }
}

//========================================================================================================

void Curvature::curvatureTensor()
{

  STensor3 TempTensor;
  STensor3 TijABTensorProduct;
  SVector3 TijAB;
  SVector3 vector_AB;

  _CurveTensor.resize(_VertexToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
    {
      MElement *e = face->getMeshElement(iElem);  //Pointer to one element
      const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element

      for (unsigned int i = 0; i<3; ++i)  // Loop over the 3 edges of each element
      {

        MVertex* A = e->getVertex(i);                   //Pointer to 1st vertex of the edge A-to-B
        MVertex* B = e->getVertex((i+1)%3);             //Pointer to 2nd vertex of the edge A-to-B

        const int V0 = _VertexToInt[A->getNum()];       //Tag of the 1st vertex of the edge A-to-B
        const int V1 = _VertexToInt[B->getNum()];       //Tag of the 2nd vertex of the edge A-to-B

        //Weight for triangle-i-th's contribution to the shape tensor:
        const double Wij0 = _TriangleArea[E] / (2 * _VertexArea[V0]);
        const double Wij1 = _TriangleArea[E] / (2 * _VertexArea[V1]);

        //Approximate Curvature "kappa" along some tangential vector T:
        vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
        const double k_nominator0 =  dot(_VertexNormal[V0], vector_AB); //Dot-product of the 2 vectors
        const double k_nominator1 = -dot(_VertexNormal[V1], vector_AB);

        const double coeff   = 2.0/vector_AB.normSq(); //normSq is the norm to the power 2
        const double KijAB_0 = coeff*k_nominator0;
        const double KijAB_1 = coeff*k_nominator1;

        //Projection of the edge vector AB on the tangential plane:
        //** For Vertex 0
        tensprod(_VertexNormal[V0], _VertexNormal[V0], TempTensor);
        TempTensor      *= -1.0; //-N*Nt
        TempTensor(0,0) +=  1.0; //I-N*Nt
        TempTensor(1,1) +=  1.0;
        TempTensor(2,2) +=  1.0;

        for (int m = 0; m<3; ++m)
        {
          TijAB(m) = 0.0;
          for (int n = 0; n<3; ++n)
          {
            TijAB(m) += TempTensor(m,n) * vector_AB(n);
          }
        }

        TijAB.normalize();
        tensprod(TijAB, TijAB, TijABTensorProduct);
        _CurveTensor[V0] += Wij0*KijAB_0*TijABTensorProduct;

        //** For Vertex 1
        tensprod(_VertexNormal[V1], _VertexNormal[V1], TempTensor);
        TempTensor      *= -1.0; //-N*Nt
        TempTensor(0,0) +=  1.0; //I-N*Nt
        TempTensor(1,1) +=  1.0;
        TempTensor(2,2) +=  1.0;

         for (int m = 0; m<3; ++m)
        {
          TijAB(m) = 0.0;
          for (int n = 0; n<3; ++n)
          {
            TijAB(m) += TempTensor(m,n) * vector_AB(n);
          }
        }

        TijAB.normalize();
        tensprod(TijAB, TijAB, TijABTensorProduct);
        _CurveTensor[V1] += Wij1*KijAB_1*TijABTensorProduct;

      }//end of loop over the edges

    }//end of loop over all elements of one GFace

  }//End of loop over ptFinalEntityList

}//End of method

//========================================================================================================

void Curvature::computeCurvature_Simple()
{
  SVector3 vector_E;
  SVector3 vector_A;
  SVector3 vector_B;
  SVector3 vector_Wvi;
  STensor3 Qvi;
  STensor3 QviT;
  STensor3 Holder;

  initializeMap();
  computeVertexNormals();
  curvatureTensor();

  _VertexCurve.resize(_VertexToInt.size());

  for (unsigned int n = 0; n < _VertexToInt.size(); ++ n) //Loop over the vertex
  {
    vector_E = SVector3(1,0,0);
    vector_A = vector_E + _VertexNormal[n];
    vector_B = vector_E - _VertexNormal[n];

    const double MagA = vector_A.norm();
    const double MagB = vector_B.norm();

    if (MagB > MagA)
    {
      vector_Wvi = vector_B;
    }
    else
    {
      vector_Wvi = vector_A;
    }

    vector_Wvi.normalize();
    //to obtain the Qvi = Id -2*Wvi*Wvi^t
    tensprod(vector_Wvi, vector_Wvi, Qvi); //Qvi = Wvi*Wvi^t
    Qvi      *= -2.0;                      //-2*Wvi*Wvi^t
    Qvi(0,0) +=  1.0;                      //I - 2*Wvi*Wvi^t  ==> Householder transformation
    Qvi(1,1) +=  1.0;
    Qvi(2,2) +=  1.0;

    //Transpose the matrix:
    for (int i = 0; i<3; ++i)
    {
      for (int j = 0; j<3; ++j)
      {
        QviT(i,j) = Qvi(j,i);
       }
     }

     QviT *= _CurveTensor[n];
     QviT *=Qvi;
     Holder = QviT;

     //Eigenvalues of the 2*2 minor from the Householder matrix
     const double A = 1.0;
     const double B = -(Holder(1,1) + Holder(2,2));
     const double C = Holder(1,1)*Holder(2,2) - Holder(1,2)*Holder(2,1);
     const double Delta = std::sqrt(B*B-4*A*C);

     if((B*B-4.*A*C) < 0.0)
     {
       std::cout << "WARNING: negative discriminant: " << B*B-4.*A*C << std::endl;
     }

     const double m11 = (-B + Delta)/(2*A);  //Eigenvalue of Householder submatrix
     const double m22 = (-B - Delta)/(2*A);

     //_VertexCurve[n] = (3*m11-m22)*(3*m22-m11);  //Gaussian Curvature
     _VertexCurve[n] = ((3*m11-m22) + (3*m22-m11))*0.5; //Mean Curvature

  }
}

//========================================================================================================

//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND

void Curvature::computeRusinkiewiczNormals()
{
  SVector3 vector_AB;
  SVector3 vector_AC;
  SVector3 vector_BC;

  _TriangleArea.resize(_ElementToInt.size() );
  _VertexNormal.resize(_VertexToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    // face is a pointer to one surface of the group "FinalEntityList"
    GFace* face = _EntityArray[i];

    //Loop over the element all the element of the "myTag"-surface
    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      // Pointer to one element
      MElement *e = face->getMeshElement(iElem);
      const int E = _ElementToInt[e->getNum()];
    
      // Pointers to vertices of triangle
      MVertex* A = e->getVertex(0);
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int V0 = _VertexToInt[A->getNum()];  //The new number of the vertex
      const int V1 = _VertexToInt[B->getNum()];
      const int V2 = _VertexToInt[C->getNum()];

      vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
      vector_AC = SVector3(C->x() - A->x(), C->y() - A->y(), C->z() - A->z() );
      vector_BC = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z() );

      const SVector3 cross = crossprod(vector_AB, vector_AC);

      // Area of the triangles:
      _TriangleArea[E] = 0.5*cross.norm();

      const double l_AB = vector_AB.normSq();
      const double l_AC = vector_AC.normSq();
      const double l_BC = vector_BC.normSq();

      _VertexNormal[V0] += cross * (1.0/ (l_AB*l_AC));
      _VertexNormal[V1] += cross * (1.0/ (l_AB*l_BC));
      _VertexNormal[V2] += cross * (1.0/ (l_AC*l_BC));

    } // end of loop over elements of one face

  } //Loop over _ptFinalEntityList

    ///////Normalize the vertex-normals.
    for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
    {
      _VertexNormal[n].normalize();
    }
}

//========================================================================================================
// Compute per-vertex point areas
void Curvature::computePointareas(){

  SVector3 e[3];
  SVector3 l2;
  SVector3 ew;
  double factor[3];

  _pointareas.resize(_VertexToInt.size());
  _cornerareas.resize(_ElementToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
    {
      MElement *E = face->getMeshElement(iElem);  //Pointer to one element
      // The NEW tag of the corresponding element
      const int EIdx = _ElementToInt[E->getNum()];

      MVertex* A = E->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = E->getVertex(1);
      MVertex* C = E->getVertex(2);

      const int V0 = _VertexToInt[A->getNum()];
      const int V1 = _VertexToInt[B->getNum()];
      const int V2 = _VertexToInt[C->getNum()];

      if (_isOnBoundary[V0])
      {
          factor[0] = 1.0;
      }
      else {factor[0] = 1.0;}

      if (_isOnBoundary[V1])
      {
          factor[1] = 1.0;
      }
      else {factor[1] = 1.0;}

      if (_isOnBoundary[V2])
      {
          factor[2] = 1.0;
      }
      else {factor[2] = 1.0;}

      //Edges
      e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangular element
      e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
      e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());

      // Compute corner weights
      //SVector3 len = crossprod(e[0], e[1]);
      //len = norm
      //len2 = normSq
      double area = 0.5 * norm(crossprod(e[0], e[1])); //area of a triangle
      l2 = SVector3( normSq(e[0]), normSq(e[1]), normSq(e[2]) );
      ew = SVector3( l2[0] * (l2[1] + l2[2] - l2[0]),
                     l2[1] * (l2[2] + l2[0] - l2[1]),
                     l2[2] * (l2[0] + l2[1] - l2[2]) );

      if (ew[0] <= 0.0)
      {
        _cornerareas[EIdx][1] = -0.25 * l2[2] * area / dot(e[0], e[2]);
        _cornerareas[EIdx][2] = -0.25 * l2[1] * area / dot(e[0], e[1]);
        _cornerareas[EIdx][0] = area - _cornerareas[EIdx][1] - _cornerareas[EIdx][2];
      }
      else if (ew[1] <= 0.0)
      {
        _cornerareas[EIdx][2] = -0.25 * l2[0] * area / dot(e[1], e[0]);
        _cornerareas[EIdx][0] = -0.25 * l2[2] * area / dot(e[1], e[2]);
        _cornerareas[EIdx][1] = area - _cornerareas[EIdx][2] - _cornerareas[EIdx][0];
      }
      else if (ew[2] <= 0.0)
      {
        _cornerareas[EIdx][0] = -0.25 * l2[1] * area / dot(e[2], e[1]);
        _cornerareas[EIdx][1] = -0.25 * l2[0] * area / dot(e[2], e[0]);
        _cornerareas[EIdx][2] = area - _cornerareas[EIdx][0] - _cornerareas[EIdx][1];
      }
      else
      {
        float ewscale = 0.5 * area / (ew[0] + ew[1] + ew[2]);
        for (int j = 0; j < 3; j++)
          _cornerareas[EIdx][j] = ewscale * (ew[(j+1)%3] + ew[(j+2)%3]);
      }

      _pointareas[V0] += _cornerareas[EIdx][0];

      _pointareas[V1] += _cornerareas[EIdx][1];

      _pointareas[V2] += _cornerareas[EIdx][2];

      for (int j = 0; j < 3; j++)
      {
          _cornerareas[EIdx][j] = factor[j]*_cornerareas[EIdx][j];
      }


    } //End of loop over iElem

//      std::cout << "_pointareas.size = " << _pointareas.size() << std::endl;
//      std::cout << "_cornerareas.size = " << _cornerareas.size() << std::endl;

  } //End of loop over _ptFinalEntityList

} //End of the method "computePointareas"


//========================================================================================================
//Rotate a coordinate system to be perpendicular to the given normal
void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v){

  new_u = old_u;
  new_v = old_v;
  SVector3 old_norm = crossprod(old_u, old_v);
  double ndot = dot(old_norm, new_norm);
//  if (unlikely(ndot <= -1.0f))
  if (ndot <= -1.0f)
  {
    new_u = -1.0*new_u;
    new_v = -1.0*new_v;
    return;
  }

  SVector3 perp_old = new_norm - ndot*old_norm;
  SVector3 dperp = 1.0f/(1 + ndot) * (old_norm + new_norm);
  new_u -= dperp*dot(new_u, perp_old);
  new_v -= dperp*dot(new_v, perp_old);
}

//========================================================================================================

//Project a curvature tensor from the basis spanned by old_u and old_v
//(which are assumed to be unit-length and perpendicular) to the new_u
//and new_v basis

void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
                          double old_ku, double old_kuv, double old_kv,
                          const SVector3  &new_u, const SVector3 &new_v,
                          double &new_ku, double &new_kuv, double &new_kv){
  SVector3 r_new_u;
  SVector3 r_new_v;
  rot_coord_sys(new_u, new_v, crossprod(old_u,old_v), r_new_u, r_new_v);

  const double u1 = dot(r_new_u, old_u);
  const double v1 = dot(r_new_u, old_v);
  const double u2 = dot(r_new_v, old_u);
  const double v2 = dot(r_new_v, old_v);

  new_ku  =   old_ku*u1*u1 + old_kuv*(2.0f * u1*v1)   + old_kv*v1*v1;
  new_kuv  =  old_ku*u1*u2 + old_kuv*(u1*v2 * u2*v1)  + old_kv*v1*v2;
  new_kv  =   old_ku*u2*u2 + old_kuv*(2.0f * u2*v2)   + old_kv*v2*v2;
}


//========================================================================================================

//Given a curvature tensor, find principal directions and curvatures
//Makes sure that pdir1 and pdir2 are perpendicular to normal

void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
                      double ku, double kuv, double kv,
                      const SVector3 &new_norm,
                      SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2){
  SVector3 r_old_u;
  SVector3 r_old_v;

  double c = 1.0;
  double s = 0.0;
  double tt = 0.0;

  rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);

//  if(unlikely(kuv !=0.0f))
  if(kuv !=0.0)  {
    //Jacobi rotation to diagonalize
    double h= 0.5*(kv -ku)/ kuv;
    tt = (h<0.0) ?
          1.0 / (h - std::sqrt(1.0 + h*h)):
          1.0 / (h + std::sqrt(1.0 + h*h));
    c = 1.0 / std::sqrt(1.0 + tt*tt);
    s = tt*c;
  }

  k1 = ku - tt *kuv;
  k2 = kv + tt *kuv;

  if(std::abs(k1) >= std::abs(k2))  {
    pdir1 = c*r_old_u - s*r_old_v;
  }
  else  {
    std::swap(k1,k2);
    pdir1 = s*r_old_u + c*r_old_v;
  }
  pdir2 = crossprod(new_norm, pdir1);
}

//========================================================================================================
void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
{

  _model = model;

  double t0 = Cpu();
  Msg::StatusBar(2, true, "(C) Computing Curvature");
  if (typ == RUSIN)
    computeCurvature_Rusinkiewicz(0);
  else if (typ == RBF)
    computeCurvature_RBF();
  else if (typ == SIMPLE)
    computeCurvature_Simple();

  double t1 = Cpu();
  Msg::StatusBar(2, true, "(C) Done Computing Curvature (%g s)", t1-t0);

  writeToMshFile("curvature.msh");
  writeToPosFile("curvature.pos");
  writeToVtkFile("curvature.vtk");

}

//========================================================================================================

void Curvature::buildEdgeList()
{

  int V[3];
  bool edgefound;

  _VertexToEdgeList.clear();

  std::list<MeshEdgeInfo>::iterator VertToEdgeListIter;
  MeshEdgeInfo TempEdge;

  _VertexToEdgeList.resize(_VertexToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i]; //List of all the discrete face on which we compute the curvature

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      MElement *e = face->getMeshElement(iElem);

      // The NEW tag of the corresponding element
      const int E = _ElementToInt[e->getNum()];

      // Pointers to vertices of triangle
      MVertex* A = e->getVertex(0);
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      V[0] = _VertexToInt[A->getNum()];  //The new number of the vertex
      V[1] = _VertexToInt[B->getNum()];
      V[2] = _VertexToInt[C->getNum()];

      // Try to add edge [V0,V1] to the global list
      for (int j = 0; j < 3; ++j)
      {
          const int StartV = std::min(V[j], V[(j+1)%3]);
          const int EndV = std::max(V[j], V[(j+1)%3]);

          edgefound = false;
          for (VertToEdgeListIter = _VertexToEdgeList[StartV].begin(); VertToEdgeListIter != _VertexToEdgeList[StartV].end();
                                                                                   ++VertToEdgeListIter)
          {
            if(StartV == (*VertToEdgeListIter).StartV && EndV == (*VertToEdgeListIter).EndV)
            {
              (*VertToEdgeListIter).NbElemNeighbour ++;
              edgefound = true;
            }

          }
          if (false == edgefound)
          {
            TempEdge.StartV = StartV;
            TempEdge.EndV = EndV;
            TempEdge.NbElemNeighbour = 1;
            _VertexToEdgeList[StartV].push_back(TempEdge);
          }


      } // Loop over j


    } // Loop over the elements (triangles) of the face
  }

  int NbEdges = 0;
  for(int i = 0; i < _VertexToEdgeList.size(); ++i)
  {
    NbEdges += _VertexToEdgeList[i].size();
  }

  //std::cout << "Euler formula:" << std::endl;
  //std::cout << "Edges + 2 =        " << NbEdges + 2 << std::endl;
  //std::cout << "Elements + Nodes = " << _VertexToInt.size() + _ElementToInt.size() << std::endl;

}

//========================================================================================================

void Curvature::smoothCurvatureField(const int NbIter)
{
  if ( _VertexToEdgeList.size() == 0 ) { buildEdgeList(); }

  std::vector<double> smoothedCurvature;
  smoothedCurvature.resize( _VertexToInt.size() );

  _NbNeighbour.resize(_VertexToInt.size());
  for(int i = 0; i < _VertexToInt.size(); ++i)
  {
    _NbNeighbour[i] = 0;
  }


  // Smoothing iterations
  for(int iter = 0; iter < NbIter; ++iter)
  {

    for(int i = 0; i < smoothedCurvature.size(); ++i)
    {
      smoothedCurvature[i] = 0.0;
    }

    std::list<MeshEdgeInfo>::const_iterator edgeIterator;


    for(int i = 0; i < _VertexToEdgeList.size(); ++i)
    {
      for(edgeIterator = _VertexToEdgeList[i].begin(); edgeIterator != _VertexToEdgeList[i].end(); ++edgeIterator)
      {
        const int N0 = (*edgeIterator).StartV;
        const int N1 = (*edgeIterator).EndV;

        const int V0 = _VertexToInt[N0];
        const int V1 = _VertexToInt[N1];

        smoothedCurvature[V0] += _VertexCurve[V1];
        smoothedCurvature[V1] += _VertexCurve[V0];
        _NbNeighbour[V0]++;
        _NbNeighbour[V1]++;
      }
    }



    const double Lambda = 0.3;
    for(int i = 0; i < _VertexCurve.size(); ++i)
    {
      _VertexCurve[i] = Lambda*_VertexCurve[i] + (1-Lambda)*smoothedCurvature[i] / _NbNeighbour[i];
    }
  }

}

//========================================================================================================



void Curvature::computeCurvature_Rusinkiewicz(int isMax)
{
  retrieveCompounds();
  initializeMap();

  computeRusinkiewiczNormals();
  computePointareas();

  SVector3 e[3];
  STensor3 w;
  SVector3 t;
  SVector3 n;
  SVector3 b;
  SVector3 m;
  SVector3 dn;

  double u;
  double v;
  double dnu;
  double dnv;
  double c1;
  double c2;
  double c12;
  double wt;

  _pdir1.resize(_VertexToInt.size());
  _pdir2.resize(_VertexToInt.size());

  _curv1.resize(_VertexToInt.size());
  _curv2.resize(_VertexToInt.size());
  _curv12.resize(_VertexToInt.size());

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
    {
      MElement *E = face->getMeshElement(iElem);  //Pointer to one element

      MVertex* A = E->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = E->getVertex(1);
      MVertex* C = E->getVertex(2);

      const int V0 = _VertexToInt[A->getNum()];  //Tag of the 1st vertex of the triangle
      const int V1 = _VertexToInt[B->getNum()];
      const int V2 = _VertexToInt[C->getNum()];

      ///Set up an initial coordinate system per vertex:

      _pdir1[V0] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
      _pdir1[V1] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
      _pdir1[V2] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
    }
  }

  for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
  {
    _pdir1[ivertex] = crossprod(_pdir1[ivertex], _VertexNormal[ivertex]);
    _pdir1[ivertex].normalize();
    _pdir2[ivertex] = crossprod(_VertexNormal[ivertex], _pdir1[ivertex]);
  }

  // Compute curvature per face:
  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    //face is a pointer to one surface of the group "FinalEntityList"
    GFace* face = _EntityArray[i];

    //Loop over all elements of this face:
    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      MElement *E = face->getMeshElement(iElem);  //Pointer to one element
      // The NEW tag of the corresponding element
      const int EIdx = _ElementToInt[E->getNum()];

      MVertex* A = E->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = E->getVertex(1);
      MVertex* C = E->getVertex(2);

      e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
      e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
      e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());

      //SVector3 e[3] = {SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()), SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z()), SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z()) };

      //N-T-B coordinate system per face
      t = e[0];
      t.normalize();
      n = crossprod( e[0], e[1]);
      b = crossprod(n, t);
      b.normalize();

      //Estimate curvature based on variations of normals along edges:
      //intialization:
      m = SVector3(0.0, 0.0, 0.0);
      //maybe double m[3] = { 0.0, 0.0, 0.0 };
      // w *= 0.0; //Reset w to zero
      for (int i  = 0; i< 3; ++i)
      {
        for (int j = 0; j<3; ++j)
        {
          w(i,j) = 0.0;
        }
      }

      //filling:
      for (int j = 0; j< 3; ++j)
      {
        u = dot(e[j], t);
        v = dot(e[j], b);

        w(0,0) += u*u;
        w(0,1) += u*v;
        w(2,2) += v*v;

        MVertex* U = E->getVertex(PREV(j));
        MVertex* V = E->getVertex(NEXT(j));

        const int UIdx = _VertexToInt[U->getNum()];
        const int VIdx = _VertexToInt[V->getNum()];

        dn = _VertexNormal[UIdx] - _VertexNormal[VIdx];

        dnu = dot(dn, t);
        dnv = dot(dn, b);

        m[0] += dnu*u;
        m[1] += dnu*v + dnv*u;
        m[2] += dnv*v;
      }

      w(1,1) = w(0,0) + w(2,2);
      w(1,2) = w(0,1);

      //Least Squares Solution
      double diag[3];
      if (!ldltdc(w, diag))
      {
        std::cout << "ldltdc failed" << std::endl;
        continue;
      }
      ldltsl(w, diag, m, m);

      //Push it back out to the vertices
      for (int j = 0; j< 3; ++j)
      {
        const int old_vj = E->getVertex(j)->getNum();
        const int vj = _VertexToInt[old_vj];
        proj_curv(t, b, m[0], m[1], m[2], _pdir1[vj], _pdir2[vj], c1, c12, c2);
        wt = _cornerareas[EIdx][j]/_pointareas[vj];
//          wt = 1.0;
//        if (_isOnBoundary[vj])
//        {
//            wt = 2.0*wt;
//        }

        _curv1[vj]  += wt*c1;
        _curv12[vj] += wt*c12;
        _curv2[vj]  += wt*c2;
      }

    } //End of loop over the element (iElem)

  } //End of loop over "_EntityArray"


  //Compute principal directions and curvatures at each vertex
  for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)  {
    diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv12[ivertex], _curv2[ivertex],
                     _VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv2[ivertex]);
  }

  _VertexCurve.resize( _VertexToInt.size() );

  for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){

    if (isMax){
      _VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
    }
    else{
    //  _VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5; //Mean curvature
    //_VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]);
    _VertexCurve[ivertex] = _curv1[ivertex]*_curv2[ivertex]; //Gaussian
    //_VertexCurve[ivertex] = std::abs(_VertexCurve[ivertex]);
    }

  }

//   smoothCurvatureField(1);
  _alreadyComputedCurvature = true;

// Propagate the value of curvature from nodes close the edges of the geometry onto the edges
  correctOnEdges();

} //End of the "computeCurvature_Rusinkiewicz" method

 //========================================================================================================

void Curvature::computeCurvature_RBF()
{
  retrieveCompounds();
  initializeMap();
  
  //fill set of MVertex
  std::set<MVertex*> allNodes;
  for (int i = 0; i< _EntityArray.size(); ++i)  {
    GFaceCompound* face = (GFaceCompound*)_EntityArray[i];
    for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++) {
      MElement *e = face->getMeshElement(iElem);
      for(int j = 0; j < e->getNumVertices(); j++){
        allNodes.insert(e->getVertex(j));
      }
    }
  }

  //bounding box
  SBoundingBox3d bb;
  std::vector<SPoint3> vertices;
  for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
    SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
    vertices.push_back(pt);
    bb +=pt;
  }
  double sizeBox = norm(SVector3(bb.max(), bb.min()));

  //compure curvature RBF
  std::map<MVertex*, SVector3> _normals;
  std::vector<MVertex*> _ordered;
  std::map<MVertex*, double> curvRBF;
  //GLOBAL
  GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered); 
  _rbf->computeCurvature(_rbf->getXYZ(),curvRBF);
  //LOCAL FD
  //GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true); 
  //_rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);

  //fill vertex curve
  _VertexCurve.resize( _VertexToInt.size() );
  for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
    MVertex *v = *itv;
    std::map<int,int>::iterator vertexIterator;
    int V0;
    vertexIterator = _VertexToInt.find( v->getNum() );
    if ( vertexIterator != _VertexToInt.end() )  V0 = (*vertexIterator).second;
    _VertexCurve[V0] = curvRBF[v];
   }

 _alreadyComputedCurvature = true;

} //End of the "computeCurvature_RBF" method


 //========================================================================================================

void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs)
  {
    MVertex* A = triangle->getVertex(0);
    MVertex* B = triangle->getVertex(1);
    MVertex* C = triangle->getVertex(2);

    int V0 = 0;
    int V1 = 0;
    int V2 = 0;

    std::map<int,int>::iterator vertexIterator;
    vertexIterator = _VertexToInt.find( A->getNum() );
    if ( vertexIterator != _VertexToInt.end() )  V0 = (*vertexIterator).second;
    else
      std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
    
    vertexIterator = _VertexToInt.find( B->getNum() );
    if ( vertexIterator != _VertexToInt.end() )   V1 = (*vertexIterator).second;
    else
      std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
    
    vertexIterator = _VertexToInt.find( C->getNum() );
    if ( vertexIterator != _VertexToInt.end() )  V2 = (*vertexIterator).second;
    else
      std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
    
    if (isAbs){
      c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
      c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
      c2 = std::abs(_VertexCurve[V2]); //Mean curvature in vertex 2
    }
    else{
      c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
      c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
      c2 = _VertexCurve[V2]; //Mean curvature in vertex 2
    }

  }

//========================================================================================================

void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs)
{
  MVertex* A = triangle->getVertex(0);
  MVertex* B = triangle->getVertex(1);
  MVertex* C = triangle->getVertex(2);

  int V0 = 0;
  int V1 = 0;
  int V2 = 0;

  std::map<int,int>::iterator vertexIterator;
  vertexIterator = _VertexToInt.find( A->getNum() );
  if ( vertexIterator != _VertexToInt.end() )  V0 = (*vertexIterator).second;
  else
    std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;

  vertexIterator = _VertexToInt.find( B->getNum() );
  if ( vertexIterator != _VertexToInt.end() )  V1 = (*vertexIterator).second;
  else
    std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;

  vertexIterator = _VertexToInt.find( C->getNum() );
  if ( vertexIterator != _VertexToInt.end() )  V2 = (*vertexIterator).second;
  else
    std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;

  if (isAbs){
    dMax[0] = _pdir1[V0];
    dMax[1] = _pdir1[V1];
    dMax[2] = _pdir1[V2];

    dMin[0] = _pdir2[V0];
    dMin[1] = _pdir2[V1];
    dMin[2] = _pdir2[V2];

    cMax[0]  = std::abs(_curv1[V0]);
    cMax[1]  = std::abs(_curv1[V1]);
    cMax[2]  = std::abs(_curv1[V2]);

    cMin[0]  = std::abs(_curv2[V0]);
    cMin[1]  = std::abs(_curv2[V1]);
    cMin[2]  = std::abs(_curv2[V2]);

  }
  else{

    dMax[0] = _pdir1[V0];
    dMax[1] = _pdir1[V1];
    dMax[2] = _pdir1[V2];

    dMin[0] = _pdir2[V0];
    dMin[1] = _pdir2[V1];
    dMin[2] = _pdir2[V2];

    cMax[0]  = _curv1[V0];
    cMax[1]  = _curv1[V1];
    cMax[2]  = _curv1[V2];

    cMin[0]  = _curv2[V0];
    cMin[1]  = _curv2[V1];
    cMin[2]  = _curv2[V2];
  }
}



  //========================================================================================================

void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
   {
     MVertex* A = edge->getVertex(0);
     MVertex* B = edge->getVertex(1);

     int V0 = 0;
     int V1 = 0;

     std::map<int,int>::iterator vertexIterator;

     vertexIterator = _VertexToInt.find( A->getNum() );
     if ( vertexIterator != _VertexToInt.end() )  V0 = (*vertexIterator).second;
     else  std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
     
     vertexIterator = _VertexToInt.find( B->getNum() );
     if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
     else std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
     
     if (isAbs){
       c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
       c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
     }
     else{
       c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
       c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
     }

   }

//========================================================================================================

double Curvature::getAtVertex(const MVertex *v) const {
  std::map<int,int>::const_iterator it = _VertexToInt.find(v->getNum());
  if (it == _VertexToInt.end()) {
    Msg::Error("curvature has not been computed for vertex %i (%i)", v->getNum(), _VertexToInt.size());
    return 1;
  }
  return _VertexCurve[it->second];
}

//========================================================================================================

void Curvature::writeToMshFile(const std::string &filename)
{
  std::ofstream outfile;
  outfile.precision(18);
  outfile.open(filename.c_str());

  /// Write the values of curvature
  outfile << "$MeshFormat"    << std::endl;
  outfile << "2.1 0 8"        << std::endl;
  outfile << "$EndMeshFormat" << std::endl;
  outfile << "$NodeData" << std::endl;
  outfile << "1" << std::endl;                 // One string tag
  outfile << "\"Curvature\"" << std::endl;     // The name of the view
  outfile << "1"             << std::endl;     // One real tag
  outfile << "0.0"           << std::endl;     // The time value
  outfile << "3"             << std::endl;     // Three integer tags
  outfile << "0"             << std::endl;     // The time step (time steps always start at 0)
  outfile << "1"             << std::endl;     // 1-component (scalar) field
  outfile << _VertexToInt.size() << std::endl; // How many associated nodal values

  std::map<int,int>::const_iterator vertex_iterator;
  for(vertex_iterator = _VertexToInt.begin(); vertex_iterator != _VertexToInt.end(); ++vertex_iterator)
  {
    outfile << vertex_iterator->first << " " << _VertexCurve[vertex_iterator->second] << std::endl;
  }

  outfile << "$EndNodeData" << std::endl;


  /// Write the values of characteristic length

  double lc;

  outfile << "$MeshFormat"    << std::endl;
  outfile << "2.1 0 8"        << std::endl;
  outfile << "$EndMeshFormat" << std::endl;
  outfile << "$NodeData" << std::endl;
  outfile << "1" << std::endl;                 // One string tag
  outfile << "\"Characteristic mesh length\"" << std::endl;     // The name of the view
  outfile << "1"             << std::endl;     // One real tag
  outfile << "0.0"           << std::endl;     // The time value
  outfile << "3"             << std::endl;     // Three integer tags
  outfile << "0"             << std::endl;     // The time step (time steps always start at 0)
  outfile << "1"             << std::endl;     // 1-component (scalar) field
  outfile << _VertexToInt.size() << std::endl; // How many associated nodal values

  for(vertex_iterator = _VertexToInt.begin(); vertex_iterator != _VertexToInt.end(); ++vertex_iterator)
  {
    lc = 2.0*M_PI/( fabs(_VertexCurve[vertex_iterator->second]) * CTX::instance()->mesh.minCircPoints );
    lc = std::max(lc, CTX::instance()->mesh.lcMin);
    lc = std::min(lc, CTX::instance()->mesh.lcMax);
    outfile << vertex_iterator->first << " " << 1.0/(lc*lc) << std::endl;
  }

  outfile << "$EndNodeData" << std::endl;

  /// Write the values of curvature direction - principal direction 1
  outfile << "$NodeData" << std::endl;
  outfile << "1" << std::endl;                 // One string tag
  outfile << "\"Principal curvature direction 1\"" << std::endl;     // The name of the view
  outfile << "1"             << std::endl;     // One real tag
  outfile << "0.0"           << std::endl;     // The time value
  outfile << "3"             << std::endl;     // Three integer tags
  outfile << "0"             << std::endl;     // The time step (time steps always start at 0)
  outfile << "3"             << std::endl;     // 3-component (vector) field
  outfile << _VertexToInt.size() << std::endl; // How many associated nodal values

  for(vertex_iterator = _VertexToInt.begin(); vertex_iterator != _VertexToInt.end(); ++vertex_iterator)
  {
    outfile << vertex_iterator->first << " " << _pdir1[vertex_iterator->second].x() << " "
                                             << _pdir1[vertex_iterator->second].y() << " "
                                             << _pdir1[vertex_iterator->second].z() << std::endl;
  }

  outfile << "$EndNodeData" << std::endl;


  /// Write the values of curvature direction - principal direction 1
  outfile << "$NodeData" << std::endl;
  outfile << "1" << std::endl;                 // One string tag
  outfile << "\"Principal curvature direction 2\"" << std::endl;     // The name of the view
  outfile << "1"             << std::endl;     // One real tag
  outfile << "0.0"           << std::endl;     // The time value
  outfile << "3"             << std::endl;     // Three integer tags
  outfile << "0"             << std::endl;     // The time step (time steps always start at 0)
  outfile << "3"             << std::endl;     // 3-component (vector) field
  outfile << _VertexToInt.size() << std::endl; // How many associated nodal values

  for(vertex_iterator = _VertexToInt.begin(); vertex_iterator != _VertexToInt.end(); ++vertex_iterator)
  {
    outfile << vertex_iterator->first << " " << _pdir2[vertex_iterator->second].x() << " "
                                             << _pdir2[vertex_iterator->second].y() << " "
                                             << _pdir2[vertex_iterator->second].z() << std::endl;
  }

  outfile << "$EndNodeData" << std::endl;

  outfile.close();


}

//========================================================================================================

void Curvature::writeToPosFile( const std::string & filename)
{
  std::ofstream outfile;
  outfile.precision(18);
  outfile.open(filename.c_str());
  outfile << "View \"Curvature \"{" << std::endl;

  int idxelem = 0;

  for (int i = 0; i< _EntityArray.size(); ++i) {
    GFace* face = _EntityArray[i];

    for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
      MElement *e = face->getMeshElement(iElem);  
      const int E = _ElementToInt[e->getNum()]; 
      //std::cout << "We are now looking at element Nr: " << E << std::endl;

      MVertex* A = e->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int V1 = _VertexToInt[A->getNum()];                //Tag of the 1st vertex of the triangle
      const int V2 = _VertexToInt[B->getNum()];                //Tag of the 2nd vertex of the triangle
      const int V3 = _VertexToInt[C->getNum()];                //Tag of the 3rd vertex of the triangle

      //Here is printing the triplet X-Y-Z of each vertex:
      //*************************************************

      outfile << "ST("; //VT = vector triangles   //ST = scalar triangle
      outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
      outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
      outfile << C->x() << ","<< C->y() << "," << C->z();

      outfile << ")";
      outfile <<"{";

      //Here is printing the 3 components of the normal vector for each vertex:
      //**********************************************************************

//      outfile << _VertexNormal[V1].x() << ","<< _VertexNormal[V1].y() << ","<< _VertexNormal[V1].z() << ",";
//      outfile << _VertexNormal[V2].x() << ","<< _VertexNormal[V2].y() << ","<< _VertexNormal[V2].z() << ",";
//      outfile << _VertexNormal[V3].x() << ","<< _VertexNormal[V3].y() << ","<< _VertexNormal[V3].z();

      outfile << _VertexCurve[V1] << "," << _VertexCurve[V2] << "," << _VertexCurve[V3];

      outfile << "};" << std::endl;

      idxelem++;

  } //Loop over elements
 
} // Loop over ptFinalEntityList

outfile << "};" << std::endl;

outfile.close();
}

//========================================================================================================

void Curvature::writeToVtkFile( const std::string & filename)
{

  std::ofstream outfile;
  outfile.precision(18);
  outfile.open(filename.c_str());
  outfile << "# vtk DataFile Version 2.0" << std::endl;
  outfile << "Surface curvature computed by Gmsh" << std::endl;
  outfile << "ASCII" << std::endl;
  outfile << "DATASET UNSTRUCTURED_GRID" << std::endl;

  const int npoints = _VertexToInt.size();

  outfile << "POINTS " << npoints << " double" << std::endl;

  /// Build a table of coordinates
  /// Loop over all elements and look at the 'old' (not necessarily continuous) numbers of vertices
  /// Get the 'new' index of each vertex through _VertexToInt and the [x,y,z] coordinates of this vertex
  /// Store them in coordx,coordy and coordz


  std::vector<VtkPoint> coord;

  coord.resize(npoints);

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i];

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      MElement *e = face->getMeshElement(iElem);

      MVertex* A = e->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int oldIdxA = A->getNum();
      const int oldIdxB = B->getNum();
      const int oldIdxC = C->getNum();

      const int newIdxA = _VertexToInt[oldIdxA];
      const int newIdxB = _VertexToInt[oldIdxB];
      const int newIdxC = _VertexToInt[oldIdxC];

      coord[newIdxA].x = A->x();
      coord[newIdxA].y = A->y();
      coord[newIdxA].z = A->z();

      coord[newIdxB].x = B->x();
      coord[newIdxB].y = B->y();
      coord[newIdxB].z = B->z();

      coord[newIdxC].x = C->x();
      coord[newIdxC].y = C->y();
      coord[newIdxC].z = C->z();
    }
  }

  for (int v = 0; v < npoints; ++v)
  {
    outfile << coord[v].x << " " << coord[v].y << " " << coord[v].z << std::endl;
  }

  /// Empty the array 'coord' to free the memory
  /// Point coordinates will not be needed anymore
  coord.clear();

  /// Write the cell connectivity

  outfile << std::endl << "CELLS " << _ElementToInt.size() << " " << 4*_ElementToInt.size() << std::endl;

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i];

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
    {
      MElement *e = face->getMeshElement(iElem);

      MVertex* A = e->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int oldIdxA = A->getNum();
      const int oldIdxB = B->getNum();
      const int oldIdxC = C->getNum();

      const int newIdxA = _VertexToInt[oldIdxA];
      const int newIdxB = _VertexToInt[oldIdxB];
      const int newIdxC = _VertexToInt[oldIdxC];

      outfile << "3 " << newIdxA << " " << newIdxB << " " << newIdxC << std::endl;
    }
  }

  outfile << std::endl << "CELL_TYPES " << _ElementToInt.size() << std::endl;
  for(int ie = 0; ie < _ElementToInt.size(); ++ie)
  {
    outfile << "5" << std::endl; //Triangle is element type 5 in vtk

  }

  /// Write the curvature values as vtk 'point data'

  outfile << std::endl << "POINT_DATA " << npoints << std::endl;
  outfile << "SCALARS Curvature float 1" << std::endl;
  outfile << "LOOKUP_TABLE default" << std::endl;

  for (int iv = 0; iv < npoints; ++iv)
  {
    outfile << _VertexCurve[iv] << std::endl;
  }

  outfile << "SCALARS CharacteristicMeshLength float 1" << std::endl;
  outfile << "LOOKUP_TABLE default" << std::endl;

  double lc;

  for (int iv = 0; iv < npoints; ++iv)
  {
    lc = 2.0*M_PI / ( fabs(_VertexCurve[iv]) * CTX::instance()->mesh.minCircPoints );
    lc = std::max(lc, CTX::instance()->mesh.lcMin);
    lc = std::min(lc, CTX::instance()->mesh.lcMax);
    outfile << 1.0/(lc*lc) << std::endl;
  }

  outfile << "VECTORS CurvatureDir1 float" << std::endl;
  for (int iv = 0; iv < npoints; ++iv)
  {
    outfile << _pdir1[iv].x() << " " << _pdir1[iv].y() << " " << _pdir1[iv].z() << std::endl;
  }

  outfile << "VECTORS CurvatureDir2 float" << std::endl;
  for (int iv = 0; iv < npoints; ++iv)
  {
    outfile << _pdir2[iv].x() << " " << _pdir2[iv].y() << " " << _pdir2[iv].z() << std::endl;
  }

  outfile.close();

}

//========================================================================================================

void Curvature::writeDirectionsToPosFile( const std::string & filename)
{
  std::ofstream outfile;
  outfile.precision(18);
  outfile.open(filename.c_str());
  outfile << "View \"Curvature_DirMax \"{" << std::endl;

  for (int i = 0; i< _EntityArray.size(); ++i)
  {
    GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"

    for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
    {
      MElement *e = face->getMeshElement(iElem);  //Pointer to one element
      const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element

      //std::cout << "We are now looking at element Nr: " << E << std::endl;

      MVertex* A = e->getVertex(0);  //Pointers to vertices of triangle
      MVertex* B = e->getVertex(1);
      MVertex* C = e->getVertex(2);

      const int V1 = _VertexToInt[A->getNum()];                //Tag of the 1st vertex of the triangle
      const int V2 = _VertexToInt[B->getNum()];                //Tag of the 2nd vertex of the triangle
      const int V3 = _VertexToInt[C->getNum()];                //Tag of the 3rd vertex of the triangle

      //Here is printing the triplet X-Y-Z of each vertex:
      //*************************************************

      outfile << "VT("; //VT = vector triangles   //ST = scalar triangle
      outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
      outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
      outfile << C->x() << ","<< C->y() << "," << C->z();

      outfile << ")";
      outfile <<"{";

      //Here is printing the 3 components of the curvature max direction for each vertex:
      //*********************************************************************************

         outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
         outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
         outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();


      outfile << "};" << std::endl;

  } //Loop over elements

} // Loop over ptFinalEntityList

outfile << "};" << std::endl;


//----------------------------------------------------------------------------------------------

outfile << "View \"Curvature_DirMin \"{" << std::endl;

for (int i = 0; i< _EntityArray.size(); ++i)
{
  GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"

  for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
  {
    MElement *e = face->getMeshElement(iElem);  //Pointer to one element
    const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element

    //std::cout << "We are now looking at element Nr: " << E << std::endl;

    MVertex* A = e->getVertex(0);  //Pointers to vertices of triangle
    MVertex* B = e->getVertex(1);
    MVertex* C = e->getVertex(2);

    const int V1 = _VertexToInt[A->getNum()];                //Tag of the 1st vertex of the triangle
    const int V2 = _VertexToInt[B->getNum()];                //Tag of the 2nd vertex of the triangle
    const int V3 = _VertexToInt[C->getNum()];                //Tag of the 3rd vertex of the triangle

    //Here is printing the triplet X-Y-Z of each vertex:
    //*************************************************

    outfile << "VT("; //VT = vector triangles   //ST = scalar triangle
    outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
    outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
    outfile << C->x() << ","<< C->y() << "," << C->z();

    outfile << ")";
    outfile <<"{";

    //Here is printing the 3 components of the curvature min direction for each vertex:
    //*********************************************************************************

       outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
       outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
       outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();


    outfile << "};" << std::endl;

} //Loop over elements

} // Loop over ptFinalEntityList

outfile << "};" << std::endl;

outfile.close();
}
//========================================================================================================