From 6091f563fa44078c676d0ec24ff2ae6e4c0003ac Mon Sep 17 00:00:00 2001
From: Emilie Sauvage <emilie.sauvage@uclouvain.be>
Date: Tue, 12 Jul 2011 14:16:16 +0000
Subject: [PATCH] Take curvature into account when remeshing edges of a surface

---
 Geo/Curvature.cpp     | 82 +++++++++++++++++++++++++++++++++++++++++--
 Geo/Curvature.h       |  8 +++--
 Geo/GFaceCompound.cpp |  2 +-
 Geo/discreteEdge.cpp  | 36 +++++++++++++++++++
 Geo/discreteEdge.h    |  1 +
 5 files changed, 124 insertions(+), 5 deletions(-)

diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 62a88f310f..0622aa7742 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -7,6 +7,7 @@
 #include "MElement.h"
 #include "GEntity.h"
 #include "GFaceCompound.h"
+#include "MLine.h"
 
 #include<iostream>
 #include<fstream>
@@ -866,7 +867,7 @@ void Curvature::computeCurvature_Rusinkiewicz()
 
 //========================================================================================================
 
- void Curvature::elementNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2)
+ void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2)
  {
    MVertex* A = triangle->getVertex(0);
    MVertex* B = triangle->getVertex(1);
@@ -921,7 +922,7 @@ void Curvature::computeCurvature_Rusinkiewicz()
 
  //========================================================================================================
 
-  void Curvature::elementNodalAbsoluteValues(MTriangle* triangle, double& c0, double& c1, double& c2)
+  void Curvature::triangleNodalAbsValues(MTriangle* triangle, double& c0, double& c1, double& c2)
   {
     MVertex* A = triangle->getVertex(0);
     MVertex* B = triangle->getVertex(1);
@@ -974,6 +975,83 @@ void Curvature::computeCurvature_Rusinkiewicz()
 
   }
 
+  //========================================================================================================
+
+   void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1)
+   {
+     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;
+     }
+
+     c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+     c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
+
+   }
+
+//========================================================================================================
+
+    void Curvature::edgeNodalAbsValues(MLine* edge, double& c0, double& c1)
+    {
+      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;
+      }
+
+      c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+      c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
+
+    }
+
+
+
+
 //========================================================================================================
 
 void Curvature::writeToPosFile( const std::string & filename)
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 3134442043..6c7d30b4f5 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -187,9 +187,13 @@ public:
   /// Code taken from Rusinkiewicz' 'trimesh2' library
   void computeCurvature_Rusinkiewicz();
 
-  void elementNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2);
+  void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2);
 
-  void elementNodalAbsoluteValues(MTriangle* triangle, double& c0, double& c1, double& c2);
+  void triangleNodalAbsValues(MTriangle* triangle, double& c0, double& c1, double& c2);
+
+  void edgeNodalValues(MLine* edge, double& c0, double& c1);
+
+  void edgeNodalAbsValues(MLine* edge, double& c0, double& c1);
 
   void writeToPosFile( const std::string & filename);
 
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1bcaf02c9e..7f067a13c2 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1421,7 +1421,7 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const
     double c1;
     double c2;
     //curvature.elementNodalValues(lt->tri,c0, c1, c2);
-    curvature.elementNodalAbsoluteValues(lt->tri,c0, c1, c2);
+    curvature.triangleNodalAbsValues(lt->tri,c0, c1, c2);
 
     double cv = (1-U-V)*c0 + U*c1 + V*c2;
 //    std::cout << "(" << c0 << "," << c1 << "," << c2 << ")" << std::endl;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index bc05d3180d..7ef8239de2 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -19,6 +19,7 @@
 #include "MPyramid.h"
 #include "Geo.h"
 #include "OS.h"
+#include "Curvature.h"
 
 discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
   : GEdge(model, num, _v0, _v1)
@@ -491,6 +492,41 @@ SVector3 discreteEdge::firstDer(double par) const
   return der;
 }
 
+double discreteEdge::curvature(double par) const
+{
+//  std::cout << "Computing Curvature in Discrete Edges" << std::endl;
+  double tLoc;
+  int iEdge;
+  if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC;
+
+  double c0, c1;
+
+
+  Curvature& curvature = Curvature::getInstance();
+
+  if( !Curvature::valueAlreadyComputed() )
+  {
+    std::cout << "Need to compute discrete curvature" << std::endl;
+    std::cout << "Getting instance of curvature" << std::endl;
+
+    curvature.setGModel( model() );
+    curvature.computeCurvature_Rusinkiewicz();
+    curvature.writeToPosFile("curvature.pos");
+    curvature.writeToVtkFile("curvature.vtk");
+
+    std::cout << " ... computing curvature finished" << std::endl;
+
+  }
+
+  curvature.edgeNodalAbsValues(lines[iEdge],c0, c1);
+
+  double cv = (1-tLoc)*c0 + tLoc*c1;
+
+  return cv;
+
+}
+
+
 Range<double> discreteEdge::parBounds(int i) const
 {
   return Range<double>(0, lines.size());
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index cbedfeb592..ec54b77993 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -23,6 +23,7 @@ class discreteEdge : public GEdge {
   virtual GeomType geomType() const { return DiscreteCurve; }
   virtual GPoint point(double p) const;
   virtual SVector3 firstDer(double par) const;
+  virtual double curvature(double par) const;
   virtual Range<double> parBounds(int) const;
   void parametrize(std::map<GFace*, std::map<MVertex*, MVertex*, 
                    std::less<MVertex*> > > &face2Verts, 
-- 
GitLab