From caa545a432d5fa31f86c8c50ab2e99f565927958 Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Tue, 20 Mar 2012 16:45:00 +0000
Subject: [PATCH] lloyd 2D et 3D

---
 Mesh/Levy3D.cpp         | 85 ++++++++++++++++++++---------------------
 Mesh/Levy3D.h           | 21 +++++-----
 Mesh/meshGFaceLloyd.cpp |  7 +++-
 3 files changed, 57 insertions(+), 56 deletions(-)

diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index af47895f74..310e51d1a0 100755
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -280,13 +280,13 @@ void VoronoiVertex::set_h(double new_h){
 
 Tensor::Tensor(){
   t11 = 1.0;
-  t12 = 0.0;
-  t13 = 0.0;
   t21 = 0.0;
-  t22 = 1.0;
-  t23 = 0.0;
   t31 = 0.0;
+  t12 = 0.0;
+  t22 = 1.0;
   t32 = 0.0;
+  t13 = 0.0;
+  t23 = 0.0;
   t33 = 1.0;
 }
 
@@ -296,32 +296,32 @@ void Tensor::set_t11(double new_t11){
   t11 = new_t11;
 }
 
-void Tensor::set_t12(double new_t12){
-  t12 = new_t12;
+void Tensor::set_t21(double new_t21){
+  t21 = new_t21;
 }
 
-void Tensor::set_t13(double new_t13){
-  t13 = new_t13;
+void Tensor::set_t31(double new_t31){
+  t31 = new_t31;
 }
 
-void Tensor::set_t21(double new_t21){
-  t21 = new_t21;
+void Tensor::set_t12(double new_t12){
+  t12 = new_t12;
 }
 
 void Tensor::set_t22(double new_t22){
   t22 = new_t22;
 }
 
-void Tensor::set_t23(double new_t23){
-  t23 = new_t23;
+void Tensor::set_t32(double new_t32){
+  t32 = new_t32;
 }
 
-void Tensor::set_t31(double new_t31){
-  t31 = new_t31;
+void Tensor::set_t13(double new_t13){
+  t13 = new_t13;
 }
 
-void Tensor::set_t32(double new_t32){
-  t32 = new_t32;
+void Tensor::set_t23(double new_t23){
+  t23 = new_t23;
 }
 
 void Tensor::set_t33(double new_t33){
@@ -332,32 +332,32 @@ double Tensor::get_t11(){
   return t11;
 }
 
-double Tensor::get_t12(){
-  return t12;
+double Tensor::get_t21(){
+  return t21;
 }
 
-double Tensor::get_t13(){
-  return t13;
+double Tensor::get_t31(){
+  return t31;
 }
 
-double Tensor::get_t21(){
-  return t21;
+double Tensor::get_t12(){
+  return t12;
 }
 
 double Tensor::get_t22(){
   return t22;
 }
 
-double Tensor::get_t23(){
-  return t23;
+double Tensor::get_t32(){
+  return t32;
 }
 
-double Tensor::get_t31(){
-  return t31;
+double Tensor::get_t13(){
+  return t13;
 }
 
-double Tensor::get_t32(){
-  return t32;
+double Tensor::get_t23(){
+  return t23;
 }
 
 double Tensor::get_t33(){
@@ -711,7 +711,7 @@ void LpCVT::verification(std::vector<SPoint3>& bank,std::vector<int>& movability
   std::vector<SVector3> gradients;
 
   gradients.resize(bank.size()-offset);
-  e = 0.0001;
+  e = 0.0000001;
   srand(time(NULL));
   index = rand()%(bank.size()-offset) + offset;
 
@@ -743,6 +743,7 @@ void LpCVT::verification(std::vector<SPoint3>& bank,std::vector<int>& movability
 
   printf("...%f  %f  %f\n",gradients[index-offset].x(),gradients[index-offset].y(),gradients[index-offset].z());
   printf("...%f  %f  %f\n",(right-left)/(2.0*e),(up-down)/(2.0*e),(front-back)/(2.0*e));
+  printf("...%d %d %d\n",index,bank.size(),offset);
 }
 
 void LpCVT::eval(std::vector<SPoint3>& bank,std::vector<int>& movability,int offset,std::vector<SVector3>& gradients,double& energy,int p){
@@ -858,6 +859,7 @@ void LpCVT::compute_parameters(){
 }
 
 double LpCVT::get_size(double x,double y,double z){
+  //if outside domain return 1.0 (or other value > 0.0)
   return 0.25;
 }
 
@@ -869,15 +871,15 @@ Tensor LpCVT::get_tensor(double x,double y,double z){
   t = Tensor();
 
   t.set_t11(1.0);
-  t.set_t12(0.0);
-  t.set_t13(0.0);
-
   t.set_t21(0.0);
-  t.set_t22(1.0);
-  t.set_t23(0.0);
-
   t.set_t31(0.0);
+  
+  t.set_t12(0.0);
+  t.set_t22(1.0);
   t.set_t32(0.0);
+  
+  t.set_t13(0.0);
+  t.set_t23(0.0);
   t.set_t33(1.0);
 
   return t;
@@ -1580,14 +1582,15 @@ void LpSmoother::improve_region(GRegion* gr)
 
   double initial_conditions[3*(bank.size()-offset)];
   double scales[3*(bank.size()-offset)];
-  factor = 2.0;
+  LpCVT instance;
+  factor = 0.9;
   for(i=0;i<(bank.size()-offset);i++){
     initial_conditions[i] = bank[i+offset].x();
 	initial_conditions[i+(bank.size()-offset)] = bank[i+offset].y();
 	initial_conditions[i+2*(bank.size()-offset)] = bank[i+offset].z();
-	scales[i] = factor*get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
-	scales[i+(bank.size()-offset)] = factor*get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
-	scales[i+2*(bank.size()-offset)] = factor*get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
+	scales[i] = factor*instance.get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
+	scales[i+(bank.size()-offset)] = factor*instance.get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
+	scales[i+2*(bank.size()-offset)] = factor*instance.get_size(bank[i+offset].x(),bank[i+offset].y(),bank[i+offset].z());
   }
   x.setcontent(3*(bank.size()-offset),initial_conditions);
   alglib_scales.setcontent(3*(bank.size()-offset),scales);
@@ -1626,10 +1629,6 @@ void LpSmoother::improve_region(GRegion* gr)
 #endif
 }
 
-double LpSmoother::get_size(double x,double y,double z){
-  return 0.25;
-}
-
 int LpSmoother::get_nbr_interior_vertices(){
   return interior_vertices.size();
 }
diff --git a/Mesh/Levy3D.h b/Mesh/Levy3D.h
index 1baa6ac96c..f5e0b82ffc 100755
--- a/Mesh/Levy3D.h
+++ b/Mesh/Levy3D.h
@@ -51,27 +51,27 @@ class VoronoiVertex{
 
 class Tensor{
  private:
-  double t11,t12,t13,t21,t22,t23,t31,t32,t33;
+  double t11,t21,t31,t12,t22,t32,t13,t23,t33;
  public:
   Tensor();
   ~Tensor();
   void set_t11(double);
-  void set_t12(double);
-  void set_t13(double);
   void set_t21(double);
-  void set_t22(double);
-  void set_t23(double);
   void set_t31(double);
+  void set_t12(double);
+  void set_t22(double);
   void set_t32(double);
+  void set_t13(double);
+  void set_t23(double);
   void set_t33(double);
   double get_t11();
-  double get_t12();
-  double get_t13();
   double get_t21();
-  double get_t22();
-  double get_t23();
   double get_t31();
+  double get_t12();
+  double get_t22();
   double get_t32();
+  double get_t13();
+  double get_t23();
   double get_t33();
 };
 
@@ -121,9 +121,8 @@ class LpSmoother{
   ~LpSmoother();
   void improve_model();
   void improve_region(GRegion*);
-  double get_size(double,double,double);
   static int get_nbr_interior_vertices();
   static MVertex* get_interior_vertex(int);
 };
 
-#endif
+#endif
\ No newline at end of file
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 5c3c288b2d..04316a03e3 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -433,8 +433,8 @@ void smoothing::optimize_face(GFace* gf){
   double ratio;
   double factor;
   lpcvt obj;
-  double initial_conditions[2*triangulator.numPoints];
-  double variables_scales[2*triangulator.numPoints];
+  double* initial_conditions = static_cast<double*>(malloc(2*triangulator.numPoints*sizeof(double)));
+  double* variables_scales = static_cast<double*>(malloc(2*triangulator.numPoints*sizeof(double)));
   alglib::ae_int_t maxits;
   alglib::minlbfgsstate state;
   alglib::minlbfgsreport rep;
@@ -541,6 +541,9 @@ void smoothing::optimize_face(GFace* gf){
                            gf->additionalVertices.end());
   // clear the list of additional vertices
   gf->additionalVertices.clear();
+	
+  free(initial_conditions);
+  free(variables_scales);
 }
 
 void smoothing::optimize_model(){
-- 
GitLab