From e82c24751280f2a9e50621ef20fe5ed9ee056b3b Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Fri, 10 Jun 2011 14:15:13 +0000
Subject: [PATCH] lpcvt

---
 Geo/GFaceCompound.h     |  2 ++
 Mesh/meshGFaceLloyd.cpp | 19 +++++++++++++++++--
 Mesh/meshGFaceLloyd.h   |  2 ++
 3 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index e3da867669..efc4338452 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -24,6 +24,8 @@
 #endif
 #define AR_MAX 5 //maximal geometrical aspect ratio
 
+class rbf;
+
 /*
 A GFaceCompound is a model face that is the compound of model faces.
 
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index c9b7d225d0..b06f3a2b1e 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -402,7 +402,7 @@ void smoothing::optimize_model(){
   for(it=model->firstFace();it!=model->lastFace();it++)
   {
     gf = *it;
-	if(gf->getNumMeshElements()>0){
+	if(gf->getNumMeshElements()>0 && gf->geomType()==GEntity::CompoundSurface){
 	  optimize_face(gf);
 	  recombineIntoQuads(gf,1,1);
 	}
@@ -645,6 +645,9 @@ SPoint2 lpcvt::seed(DocRecord& triangulator,GFace* gf){
   int index2;
   double x,y;
   SPoint2 x0,x1,x2;
+	
+  work_around = gf;	
+	
   for(i=0;i<triangulator.numPoints;i++){
     if(interior(triangulator,gf,i)){
 	  num = triangulator._adjacencies[i].t_length;
@@ -1126,7 +1129,7 @@ double lpcvt::get_rho(SPoint2 point,int p){
 	
   x = point.x();
   y = point.y();
-  h = backgroundMesh::current()->operator()(x,y,0.0); //0.1;
+  h = backgroundMesh::current()->operator()(x,y,0.0)*ratio(point); //0.1;
   if(h>=0.0){
     rho = pow_int(1.0/h,p+1); //integer only, theoric value : p+2
   }
@@ -1220,6 +1223,18 @@ double lpcvt::drho_dy(SPoint2 point,int p){
   return val;
 }
 
+double lpcvt::ratio(SPoint2 point){
+  double val;
+  double uv[2];
+  double metric[3];
+	
+  uv[0] = point.x();
+  uv[1] = point.y();
+  buildMetric(work_around,uv,metric);
+  val = 1.0/pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
+  return val;
+}
+
 void lpcvt::write(DocRecord& triangulator,GFace* gf,int p){
   int i;
   double energy;
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index 58942d703c..045e6c1a7d 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/Mesh/meshGFaceLloyd.h
@@ -46,6 +46,7 @@ class lpcvt{
   fullMatrix<double> gauss_weights;
   std::vector<metric> metrics;
   int gauss_num;
+  GFace* work_around;
  public :
   lpcvt();
   ~lpcvt();
@@ -81,6 +82,7 @@ class lpcvt{
   double get_rho(SPoint2,int);
   double drho_dx(SPoint2,int);
   double drho_dy(SPoint2,int);
+  double ratio(SPoint2);
   void write(DocRecord&,GFace*,int);
   void eval(DocRecord&,std::vector<SVector3>&,double&,int);
   void swap();
-- 
GitLab