diff --git a/Cassiopee/CPlot/apps/tkCADFix.py b/Cassiopee/CPlot/apps/tkCADFix.py
index c83251e11..1219576e3 100644
--- a/Cassiopee/CPlot/apps/tkCADFix.py
+++ b/Cassiopee/CPlot/apps/tkCADFix.py
@@ -26,15 +26,15 @@ def readCAD(event=None):
         faces[2] = []
     hook = OCC.readCAD(fileName, fileFmt)
     # Previous hmax, hausd?
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     if hmax is None or (hmax < 0 and hausd < 0):
         (hmax,hmin,hausd) = OCC.occ.analyseEdges(hook)
     OCC._setCADcontainer(CTK.t, fileName, fileFmt, hmax, hausd)
     CTK.CADHOOK = hook
     # remesh and redisplay
     CTK.setCursor(2, WIDGETS['frame'])
-    OCC._meshAllEdges(hook, CTK.t, hmax=hmax, hausd=hausd)
-    OCC._meshAllFacesTri(hook, CTK.t, hmax=hmax, hausd=hausd)
+    OCC._meshAllEdges(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
+    OCC._meshAllFacesTri(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
     CTK.setCursor(0, WIDGETS['frame'])
     CTK.display(CTK.t)
     CTK.TXT.insert('START', 'CAD loaded from %s.\n'%fileName)
@@ -55,7 +55,7 @@ def sewCAD(event=None):
     if CTK.CADHOOK is None: return
     tol = CTK.varsFromWidget(VARS[2].get(), 1)[0]
     hook = CTK.CADHOOK
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     faces = []
     nzs = CPlot.getSelectedZones()
     for nz in nzs:
@@ -79,8 +79,8 @@ def sewCAD(event=None):
     edges[2] = []
     faces = Internal.getNodeFromName1(CTK.t, 'FACES')
     faces[2] = []
-    OCC._meshAllEdges(hook, CTK.t, hmax=hmax, hausd=hausd) # loose manual remeshing
-    OCC._meshAllFacesTri(hook, CTK.t, hmax=hmax, hausd=hausd)
+    OCC._meshAllEdges(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd) # loose manual remeshing
+    OCC._meshAllFacesTri(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
 
     CTK.setCursor(0, WIDGETS['frame'])
     CTK.setCursor(0, WIDGETS['sewingButton'])
@@ -99,7 +99,7 @@ def filletCAD(event=None):
     if CTK.CADHOOK is None: return
     radius = CTK.varsFromWidget(VARS[3].get(), 1)[0]
     hook = CTK.CADHOOK
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     # Get selected edges
     nzs = CPlot.getSelectedZones()
     edges = []
@@ -128,8 +128,8 @@ def filletCAD(event=None):
     edges[2] = []
     faces = Internal.getNodeFromName1(CTK.t, 'FACES')
     faces[2] = []
-    OCC._meshAllEdges(hook, CTK.t, hmax=hmax, hausd=hausd) # loose manual remeshing...
-    OCC._meshAllFacesTri(hook, CTK.t, hmax=hmax, hausd=hausd)
+    OCC._meshAllEdges(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd) # loose manual remeshing...
+    OCC._meshAllFacesTri(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
     CTK.setCursor(0, WIDGETS['frame'])
     CTK.setCursor(0, WIDGETS['filletButton'])
 
@@ -143,7 +143,7 @@ def removeFaces(event=None):
     import OCC.PyTree as OCC
     if CTK.CADHOOK is None: return
     hook = CTK.CADHOOK
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     # Get selected faces
     nzs = CPlot.getSelectedZones()
     faces = []
@@ -198,7 +198,7 @@ def fillHole(event=None):
     import OCC.PyTree as OCC
     if CTK.CADHOOK is None: return
     hook = CTK.CADHOOK
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     # Get selected edges
     nzs = CPlot.getSelectedZones()
     edges = []
@@ -232,8 +232,8 @@ def fillHole(event=None):
     edges[2] = []
     faces = Internal.getNodeFromName1(CTK.t, 'FACES')
     faces[2] = []
-    OCC._meshAllEdges(hook, CTK.t, hmax=hmax, hausd=hausd) # loose manual remeshing...
-    OCC._meshAllFacesTri(hook, CTK.t, hmax=hmax, hausd=hausd)
+    OCC._meshAllEdges(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd) # loose manual remeshing...
+    OCC._meshAllFacesTri(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
     CTK.setCursor(0, WIDGETS['frame'])
     CTK.setCursor(0, WIDGETS['fillHoleButton'])
 
@@ -267,7 +267,7 @@ def setTrimFace1():
             no = Internal.getNodeFromName1(CAD, 'no')
             no = Internal.getValue(no)
             selected += str(no)+' '
-    print(selected, flush=True)
+    #print(selected, flush=True)
     VARS[5].set(selected)
 
 #==============================================================================
@@ -275,7 +275,7 @@ def trimFaces(event=None):
     import OCC.PyTree as OCC
     if CTK.CADHOOK is None: return
     hook = CTK.CADHOOK
-    [hmax, hausd] = OCC.getCADcontainer(CTK.t)
+    [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t)
     # Get selected faces
     nzs = CPlot.getSelectedZones()
     faces1 = []
@@ -316,8 +316,8 @@ def trimFaces(event=None):
     edges[2] = []
     faces = Internal.getNodeFromName1(CTK.t, 'FACES')
     faces[2] = []
-    OCC._meshAllEdges(hook, CTK.t, hmax=hmax, hausd=hausd) # loose manual remeshing...
-    OCC._meshAllFacesTri(hook, CTK.t, hmax=hmax, hausd=hausd)
+    OCC._meshAllEdges(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd) # loose manual remeshing...
+    OCC._meshAllFacesTri(hook, CTK.t, hmin=hmin, hmax=hmax, hausd=hausd)
 
     CTK.setCursor(0, WIDGETS['frame'])
     CTK.setCursor(0, WIDGETS['trimFacesButton'])
diff --git a/Cassiopee/Converter/Converter/PyTree.py b/Cassiopee/Converter/Converter/PyTree.py
index 84652de79..5a6f2abe1 100755
--- a/Cassiopee/Converter/Converter/PyTree.py
+++ b/Cassiopee/Converter/Converter/PyTree.py
@@ -1063,7 +1063,7 @@ def convertFile2PyTree(fileName, format=None, nptsCurve=20, nptsLine=2,
                 # auto setting
             (hmin,hmax,hausd) = OCC.occ.analyseEdges(hook)
         CTK.CADHOOK = hook
-        t = OCC.meshAll(hook, hmax, hmax, -1.) # constant hmax
+        t = OCC.meshAll(hook, hmax, hmax, hausd) # constant hmax
         _upgradeTree(t)
         return t
 
diff --git a/Cassiopee/KCore/KCore/Nuga/include/DelaunayMath.h b/Cassiopee/KCore/KCore/Nuga/include/DelaunayMath.h
index 48c70e00b..90b4a44c1 100755
--- a/Cassiopee/KCore/KCore/Nuga/include/DelaunayMath.h
+++ b/Cassiopee/KCore/KCore/Nuga/include/DelaunayMath.h
@@ -31,10 +31,10 @@ namespace K_LINEAR
   {
   public:
     static void eigen_vectors 
-      (E_Float a00, E_Float a11, E_Float a10, E_Float& lambda0, E_Float& lambda1, E_Float* v0, E_Float* v1) ;
+      (E_Float a00, E_Float a11, E_Float a10, E_Float& lambda0, E_Float& lambda1, E_Float* v0, E_Float* v1);
 
     static void eigen_values 
-      (E_Float a00, E_Float a11, E_Float a10, E_Float& lambda0, E_Float& lambda1) ;
+      (E_Float a00, E_Float a11, E_Float a10, E_Float& lambda0, E_Float& lambda1);
 
     static void simultaneous_reduction
       (const K_FLD::FloatArray& M1, const K_FLD::FloatArray& M2,
diff --git a/Cassiopee/KCore/KCore/Nuga/include/Metric.h b/Cassiopee/KCore/KCore/Nuga/include/Metric.h
index 8dab36258..acd0c3947 100755
--- a/Cassiopee/KCore/KCore/Nuga/include/Metric.h
+++ b/Cassiopee/KCore/KCore/Nuga/include/Metric.h
@@ -64,7 +64,7 @@ namespace DELAUNAY{
 
     VarMetric(K_FLD::FloatArray& pos, E_Float hmin, E_Float hmax, eInterpType interp_type = LINEAR);
 
-    virtual ~VarMetric(void){if (_interpol){delete _interpol; _interpol = 0;} }
+    virtual ~VarMetric(void){if (_interpol) {delete _interpol; _interpol=0;} }
 
     VarMetric& operator=(const  VarMetric& rhs) { _hmin = _hmax = -1; _field = rhs._field; _pos = rhs._pos; _interpol = nullptr; return *this; }
 
@@ -810,10 +810,10 @@ namespace DELAUNAY{
     assert (isValidMetric(mj));
 #endif
 
-    E_Float v[2]; E_Float vi[2]; E_Float vj[2];
-    //E_Float* v = _pta2;
-    //E_Float* vi = _ptb2;
-    //E_Float* vj = _ptc2;
+    E_Float v[3]; E_Float vi[3]; E_Float vj[3];
+    //E_Float* v = _pta3;
+    //E_Float* vi = _ptb3;
+    //E_Float* vj = _ptc3;
 
     NUGA::diff<3> (_pos->col(Nj), _pos->col(Ni), v);
 
diff --git a/Cassiopee/OCC/OCC/Atomic/analyse.cpp b/Cassiopee/OCC/OCC/Atomic/analyse.cpp
index c1056ea6a..6e503296b 100644
--- a/Cassiopee/OCC/OCC/Atomic/analyse.cpp
+++ b/Cassiopee/OCC/OCC/Atomic/analyse.cpp
@@ -68,12 +68,12 @@ PyObject* K_OCC::analyseEdges(PyObject* self, PyObject* args)
   }
   E_Float emean = ltot / edges.Extent();
 
-  // calcul hmax, hausd : 30 pts par edges moyen
-  E_Float hmax = emean / 30.;
+  // calcul hmax, hausd : 20 pts par edges moyen
+  E_Float hmax = emean / 20.;
   // calcul du nbre de points sur la longueur totale
   E_Int Np = ltot / hmax;
   if (Np > 20000) hmax = ltot / 20000.;
-  E_Float hmin = emin / 30.;
+  E_Float hmin = emin / 20.;
   E_Float hausd = hmax / 10.;
   printf("INFO: suggested hmin=%g hmax=%g hausd=%g\n", hmin, hmax, hausd);
   return Py_BuildValue("ddd", hmin, hmax, hausd);
diff --git a/Cassiopee/OCC/OCC/Atomic/evalFace.cpp b/Cassiopee/OCC/OCC/Atomic/evalFace.cpp
index b3a97b9c2..edeef8fa7 100644
--- a/Cassiopee/OCC/OCC/Atomic/evalFace.cpp
+++ b/Cassiopee/OCC/OCC/Atomic/evalFace.cpp
@@ -30,17 +30,17 @@
 void evalFace__(E_Int npts, E_Float* u, E_Float* v, const TopoDS_Face& F,
                 E_Float* x, E_Float* y, E_Float* z)
 {
-    Handle(Geom_Surface) face = BRep_Tool::Surface(F);
+  Handle(Geom_Surface) face = BRep_Tool::Surface(F);
 #pragma omp parallel
-{
+  {
     gp_Pnt Pt;
 #pragma omp for 
     for (E_Int i = 0; i < npts; i++)
     {
-        face->D0(u[i], v[i], Pt);
-        x[i] = Pt.X(); y[i] = Pt.Y(); z[i] = Pt.Z();
+      face->D0(u[i], v[i], Pt);
+      x[i] = Pt.X(); y[i] = Pt.Y(); z[i] = Pt.Z();
     }
-}
+  }
 }
 
 // evalFace
diff --git a/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp b/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp
index 6bcfb4623..7455f56ba 100644
--- a/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp
+++ b/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp
@@ -19,6 +19,7 @@
 
 #include "occ.h"
 
+#include "Nuga/include/DelaunayMath.h"
 #include "TopoDS.hxx"
 #include "BRep_Tool.hxx"
 #include "TopoDS_Wire.hxx"
@@ -35,6 +36,12 @@
 #include "BRepBuilderAPI_MakeFace.hxx"
 #include "BRepGProp.hxx"
 #include "GProp_GProps.hxx"
+#include "TopTools_IndexedDataMapOfShapeListOfShape.hxx"
+#include "TopTools_ListIteratorOfListOfShape.hxx"
+#include "Geom_Surface.hxx"
+#include "TopExp.hxx"
+#include "GeomLProp_SLProps.hxx"
+#include <Precision.hxx>
 
 // ultimate (best) functions
 
@@ -217,7 +224,7 @@ void geom3(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u
 void geom4(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& ue)
 {
   E_Float r = (u1-u0+h1-h0)/(u1-u0);
-  printf("r=%f\n", r);
+  //printf("r=%f\n", r);
   E_Float a = log(r);
   if (a > 1.e-12) // r!=1
   {
@@ -238,7 +245,7 @@ void geom4(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u
   }
   //for (E_Int i = 1; i < N; i++) ue[i] = ue[i-1] + h0*pow(r, i-1);  
   //ue[N-1] = u1; // force
-  for (E_Int i = 0; i < N; i++) printf("%d %f\n", i, ue[i]);
+  //for (E_Int i = 0; i < N; i++) printf("%d %f\n", i, ue[i]);
   printf("h0/r=%f real=%f\n", h0/r, ue[1]-ue[0]);
   printf("h1=%f real=%f\n", h1, ue[N-1]-ue[N-2]);
 }
@@ -410,9 +417,16 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
   GeomAdaptor_Curve geomAdap(C0.Curve());
   Standard_Real u0 = geomAdap.FirstParameter();
   Standard_Real u1 = geomAdap.LastParameter();
+  if (BRep_Tool::Degenerated(E))
+  { 
+    nbPoints = 2;
+    ue = new E_Float [nbPoints];
+    for (E_Int i = 0; i < nbPoints; i++) ue[i] = u0;
+    return 1; 
+  }
   E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1);
   L = K_FUNC::E_max(L, 1.e-12);
-
+  
   // Compute local h
   Handle(Geom_Curve) curve;
   Standard_Real first, last;
@@ -421,10 +435,8 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
   gp_Vec dE, d2E;
   E_Float U, dE1, dE2, dE3, d2E1, d2E2, d2E3;
   E_Float t1, t2, t3, a, b, rho;
-  // Compute the first derivative
-  //curve->D1(U, point, dE);
 
-  E_Int npts = 3;
+  E_Int npts = 2;
   std::vector<E_Float> Us(npts);
   for (E_Int i = 0; i < npts; i++) Us[i] = u0+(i/(npts-1.))*(u1-u0);
   std::vector<E_Float> h(npts);
@@ -454,10 +466,10 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
     else rho = 1.e12;
 
     // Compute local h
-    h[i] = std::sqrt(8.*hausd*rho);
+    h[i] = std::sqrt(1.*hausd*rho);
     h[i] = K_FUNC::E_min(h[i], hmax);
     h[i] = K_FUNC::E_max(h[i], hmin);
-    printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(8.*hausd*rho), h[i]);
+    //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(8.*hausd*rho), h[i]);
     h[i] = h[i]*(u1-u0)/L;
   }
 
@@ -487,11 +499,17 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
     N += Np[i]-1;
   }
   ue[N+1] = u1;
-
-  //printf("final: %g %g - ", u0, u1);
+  
+  //printf("final: %g %g \n", u0, u1);
   //for (E_Int i = 0; i < Ntot; i++) printf("%g ", ue[i]);
   //printf("\n");
 
+  // verifie que la suite est bien croissante
+  for (E_Int i = 0; i < Ntot-1; i++) 
+  {
+    if (ue[i+1] < ue[i]) printf("warning: %g ", ue[i]);
+  }
+
   for (E_Int i = 0; i < npts; i++)
   {
     delete [] uel[i];
@@ -499,6 +517,256 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
   return 1;
 }
 
+// ===============================================================================
+// Return the nbPoints and ue for meshing E with hmin/hmax/hausd with max of faces
+// ===============================================================================
+E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue,
+  TopoDS_Shape* shape)
+{
+  TopTools_IndexedDataMapOfShapeListOfShape edgeFaceMap;
+  TopExp::MapShapesAndAncestors(*shape, TopAbs_EDGE, TopAbs_FACE, edgeFaceMap);
+
+  BRepAdaptor_Curve C0(E);
+  GeomAdaptor_Curve geomAdap(C0.Curve());
+  Standard_Real u0 = geomAdap.FirstParameter();
+  Standard_Real u1 = geomAdap.LastParameter();
+  //printf("===============edge =========================\n");
+  //printf("u0=%g u1=%g\n", u0, u1);
+  if (BRep_Tool::Degenerated(E))
+  { 
+    nbPoints = 2;
+    ue = new E_Float [nbPoints];
+    for (E_Int i = 0; i < nbPoints; i++) ue[i] = u0;
+    return 1; 
+  }
+  E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1);
+  L = K_FUNC::E_max(L, 1.e-12);
+
+  gp_Pnt2d Puv;
+  E_Float u,v;
+  E_Float ues; 
+  Standard_Real aFirst = C0.FirstParameter(), aEnd=C0.LastParameter();
+  //printf("afirst=%g aend=%g\n", aFirst, aEnd);
+  Standard_Real pFirst = aFirst, pEnd=aEnd;
+  
+  gp_Pnt point;
+  gp_Vec DU1, DV1, DU2, DV2, DUV;
+  E_Float nx, ny, nz, n2;
+  E_Float G11, G12, G21, G22, detG;
+  E_Float M11, M12, M21, M22, detM;
+  E_Float rho, hh, K1, K2;
+  E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS;
+
+  E_Int npts = 2;
+  std::vector<E_Float> Us(npts);
+  for (E_Int i = 0; i < npts; i++) Us[i] = i/(npts-1.);
+  
+  //Us[0] = Us[0]+0.001*(Us[1]-Us[0]);
+  //Us[npts-1] = Us[npts-1]-0.001*(Us[npts-1]-Us[npts-2]);
+  
+  std::vector<E_Float> h(npts);
+  for (E_Int i = 0; i < npts; i++) h[i] = K_CONST::E_MAX_FLOAT;
+  
+  // Find faces of E
+  const TopTools_ListOfShape& connectedFaces = edgeFaceMap.FindFromKey(E);
+  for (TopTools_ListIteratorOfListOfShape it(connectedFaces); it.More(); it.Next()) 
+  {
+    const TopoDS_Face& F = TopoDS::Face(it.Value());
+    Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
+    
+    Standard_Real BU1, BU2, BV1, BV2; surface->Bounds(BU1, BU2, BV1, BV2);
+    //printf("face bound = %g %g and %g %g\n", BU1, BU2, BV1, BV2);
+
+    // p curve on F
+    Handle(Geom2d_Curve) pCurve = BRep_Tool::CurveOnSurface(E, F, pFirst, pEnd);
+    //printf("pfirst=%g pend=%g\n", pFirst, pEnd);
+    //printf("pCurve = %p\n", (void*)pCurve);
+    // estimate at some points
+    for (E_Int i = 0; i < npts; i++)
+    {
+      ues = Us[i]*(pEnd-pFirst)+pFirst;
+
+// must be MINE (0) or OCC (1)
+#define CURVATURE 1
+#if CURVATURE == 0
+      pCurve->D0(ues, Puv);
+      u = Puv.X(); v = Puv.Y(); // u,v edge on surface
+      if (u <= BU1) u = BU1+1.e-1;
+      if (u >= BU2) u = BU2-1.e-1;
+      if (v <= BV1) v = BV1+1.e-1;
+      if (v >= BV2) v = BV2-1.e-1;
+      
+      //printf("surf: u=%g, v=%g\n", u, v);
+      // estimer les derivees sur la surface
+      surface->D2(u, v, point, DU1, DV1, DU2, DV2, DUV);
+      //printf("point=%g %g %g\n", point.X(), point.Y(), point.Z());
+      //printf("dU1=%g %g %g\n", DU1.X(), DU1.Y(), DU1.Z());
+      //printf("dV1=%g %g %g\n", DV1.X(), DV1.Y(), DV1.Z());
+      
+      // get normal
+      nx = DU1.Y()*DV1.Z()-DU1.Z()*DV1.Y();
+      ny = DU1.Z()*DV1.X()-DU1.X()*DV1.Z();
+      nz = DU1.X()*DV1.Y()-DU1.Y()*DV1.X();
+      n2 = nx*nx+ny*ny+nz*nz;
+      n2 = std::sqrt(n2);
+      if (n2 < 1.e-24) n2 = 1.e12;
+      else n2 = 1./n2;
+      nx = nx*n2;
+      ny = ny*n2;
+      nz = nz*n2;
+      //printf("normale= %g %g %g norm=%g\n", nx, ny, nz, nx*nx+ny*ny+nz*nz);
+
+      // Metric
+      M11 = DU1.X()*DU1.X()+DU1.Y()*DU1.Y()+DU1.Z()*DU1.Z();
+      M22 = DV1.X()*DV1.X()+DV1.Y()*DV1.Y()+DV1.Z()*DV1.Z();
+      M12 = M21 = DU1.X()*DV1.X()+DU1.Y()*DV1.Y()+DU1.Z()*DV1.Z();
+      //printf("matrice M:\n");
+      //printf("  %g %g\n", M11, M12);
+      //printf("  %g %g\n", M21, M22);
+
+      detM = M11*M22-M12*M21;
+      if (std::abs(detM) < 1.e-12) det = 1.e12; // invalid metric
+      else det = 1./detM;
+      Mi11 = det*M22;
+      Mi21 = -det*M21;
+      Mi12 = -det*M12;
+      Mi22 = det*M11;
+      detM = std::abs(detM);
+      //printf("detM=%g\n", detM);
+
+      // matrice de courbure
+      G11 = DU2.X()*nx+DU2.Y()*ny+DU2.Z()*nz;
+      G22 = DV2.X()*nx+DV2.Y()*ny+DV2.Z()*nz;
+      G12 = G21 = DUV.X()*nx+DUV.Y()*ny+DUV.Z()*nz;
+      detG = G11*G22-G12*G21;
+      detG = std::abs(detG);
+      
+      S11 = G11*Mi11 + G12*Mi21;
+      S12 = G11*Mi12 + G12*Mi22;
+      S21 = G21*Mi11 + G22*Mi21;
+      S22 = G21*Mi12 + G22*Mi22;
+      
+      detS = S11*S22-S12*S21;
+      //printf("dU2=%g %g %g\n", DU2.X(), DU2.Y(), DU2.Z());
+      //printf("dV2=%g %g %g\n", DV2.X(), DV2.Y(), DV2.Z());
+      //printf("dUV=%g %g %g\n", DUV.X(), DUV.Y(), DUV.Z());
+      //printf("detG=%g\n", detG);
+
+      // courbure de gauss - rayon K1*K2
+      //if (detS > 1.e-12) rho = 1./detS;
+      //else rho = 1.e12;
+      //printf("gauss=%g ou bien %g\n", detG/detM, detS);
+
+      // K1 et K2
+      //printf("matrice G:\n");
+      //printf("  %g %g\n", G11, G12);
+      //printf("  %g %g\n", G21, G22);
+      
+      K_LINEAR::DelaunayMath sl;
+      sl.eigen_values(S11, S22, S12, K2, K1);
+      //printf("K1=%g K2=%g K1K2=%g\n", K1, K2, K1*K2);
+      //printf("mean=%g\n", 0.5*(S11+S22));
+
+      rho = K_FUNC::E_max(std::abs(K1), std::abs(K2));
+      if (rho > 1.e-12) rho = 1./rho;
+      else rho = 1.e12;
+
+#else
+      // with open cascade
+      GeomLProp_SLProps props(surface, u, v, 2, Precision::Confusion());
+
+      if (props.IsCurvatureDefined()) 
+      {
+        Standard_Real gaussianCurvature = props.GaussianCurvature();
+        Standard_Real meanCurvature = props.MeanCurvature();
+
+        //std::cout << "OCC: gauss: " << gaussianCurvature << std::endl;
+        //std::cout << "OCC: mean: " << meanCurvature << std::endl;
+
+        // Principal curvatures
+        K1 = meanCurvature + sqrt(meanCurvature * meanCurvature - gaussianCurvature);
+        K2 = meanCurvature - sqrt(meanCurvature * meanCurvature - gaussianCurvature);
+        //std::cout << "OCC: k1: " << K1 << " k2: " <<K2 << std::endl;
+        rho = K_FUNC::E_max(std::abs(K1), std::abs(K2));
+        if (rho > 1.e-12) rho = 1./rho;
+        else rho = 1.e12;
+      } 
+      else 
+      {
+        rho = 1.e12;
+        //std::cout << "OCC: Curvature is not defined at the given parameters." << std::endl;
+      }
+#endif
+
+      hh = std::sqrt(8.*hausd*rho);
+      hh = K_FUNC::E_min(hh, hmax);
+      hh = K_FUNC::E_max(hh, hmin);
+      //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(1.*hausd*rho), hh);
+      h[i] = K_FUNC::E_min(h[i], hh);
+    }
+  }
+
+  //printf("hi: ");
+  //for (E_Int i = 0; i < npts; i++) 
+  //{
+  //  printf("%g ", h[i]);
+  //}
+  //printf("\n");
+
+  for (E_Int i = 0; i < npts; i++)
+  {
+    // passage en parametres
+    h[i] = h[i]*(u1-u0)/L;
+  }
+
+  E_Int N;
+  std::vector<E_Float*> uel(npts);
+  std::vector<E_Int> Np(npts);
+
+  E_Int Ntot = 0;
+  for (E_Int i = 0; i < npts-1; i++)
+  {
+    // geom2
+    geom2(Us[i]*(aEnd-aFirst)+aFirst, Us[i+1]*(aEnd-aFirst)+aFirst, h[i], h[i+1], N, uel[i]);
+    Np[i] = N;
+    Ntot += N-1;
+  }
+  Ntot += 1;
+  geom1(Us[npts-2]*(aEnd-aFirst)+aFirst, Us[npts-1]*(aEnd-aFirst)+aFirst, h[npts-2], h[npts-1], N, uel[npts-1]);
+
+  // reassemble
+  ue = new E_Float [Ntot];
+  nbPoints = Ntot;
+  N = 0;
+  for (E_Int i = 0; i < npts; i++)
+  {
+    E_Float* uep = uel[i];
+    for (E_Int j = 0; j < Np[i]-1; j++) ue[j+N] = uep[j];
+    N += Np[i]-1;
+  }
+  ue[N+1] = u1;
+
+  // verifie que la suite est bien croissante
+  for (E_Int i = 0; i < Ntot-1; i++) 
+  {
+    if (ue[i+1] < ue[i]) printf("warning: %g ", ue[i]);
+  }
+
+  //printf("final: ");
+  //for (E_Int i = 0; i < Ntot; i++) 
+  //{
+  //  printf("%g ", ue[i]);
+  //}
+  //printf("\n");
+
+  for (E_Int i = 0; i < npts; i++)
+  {
+    delete [] uel[i];
+  }
+
+  return 1;
+
+}
 // ============================================================================
 // Return ue for meshing with given param in [0,1]
 // ============================================================================
@@ -627,7 +895,8 @@ E_Int __meshEdgeByFace(const TopoDS_Edge& E, const TopoDS_Face& F,
 PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
 {
   PyObject* hook; E_Int i;
-  E_Float hmin; E_Float hmax; E_Float hausd; E_Int N; PyObject* externalEdge;
+  E_Float hmin; E_Float hmax; E_Float hausd; 
+  E_Int N; PyObject* externalEdge;
   if (!PYPARSETUPLE_(args, O_ I_ RRR_ I_ O_, &hook, &i, &hmin, &hmax, &hausd, &N, &externalEdge)) return NULL;
     
   void** packet = NULL;
@@ -641,6 +910,10 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
   TopTools_IndexedMapOfShape& edges = *(TopTools_IndexedMapOfShape*)packet[2];
   const TopoDS_Edge& E = TopoDS::Edge(edges(i));
   
+  // use neighbouring faces to evaluate curvature of edge
+  bool useFaces = true;
+  TopoDS_Shape* shape = (TopoDS_Shape*)packet[0];
+
   E_Int nbPoints = 0; // nbre of points of discretized edge
   E_Float* ue; // edge param
 
@@ -655,7 +928,9 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
     RELEASESHAREDS(o, f);
     return o;
   }
-  else if (hmax > 0 && hausd < 0 && externalEdge == Py_None) // pure hmax
+  else if ( ((hausd > 0 && std::abs(hmax-hmin) < 1.e-12 && hmax > 0) ||
+            (hausd < 0 && hmax > 0)) &&
+            externalEdge == Py_None ) // pure hmax
   {
     __getParamHmax(E, hmax, nbPoints, ue);
     PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1);
@@ -675,10 +950,9 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
     RELEASESHAREDS(o, f);
     return o;
   }
-  else if (hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd
+  else if (not useFaces && hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd
   {
-    // pour l'instant on retourne hmax comme pour les mailleurs precedents
-    //__getParamHmax(E, hmax, nbPoints, ue);
+    // courbure uniquement sur edge
     __getParamHminHmaxHausd(E, hmin, hmax, hausd, nbPoints, ue);
     PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1);
     FldArrayF* f; K_ARRAY::getFromArray2(o, f);
@@ -687,6 +961,18 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
     RELEASESHAREDS(o, f);
     return o;
   }
+  else if (useFaces && hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd
+  {
+    // courbure sur le max des faces
+    __getParamHminHmaxHausd2(E, hmin, hmax, hausd, nbPoints, ue, shape);
+    PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1);
+    FldArrayF* f; K_ARRAY::getFromArray2(o, f);
+    __meshEdge(E, nbPoints, ue, *f, false);
+    delete [] ue;
+    RELEASESHAREDS(o, f);
+    return o;
+  }
+
   else if (externalEdge != Py_None) // external parametrization
   {
     E_Int ni, nj, nk;
diff --git a/Cassiopee/OCC/OCC/Atomic/trimesh.cpp b/Cassiopee/OCC/OCC/Atomic/trimesh.cpp
index 98914ba07..670d11739 100644
--- a/Cassiopee/OCC/OCC/Atomic/trimesh.cpp
+++ b/Cassiopee/OCC/OCC/Atomic/trimesh.cpp
@@ -104,7 +104,8 @@ PyObject* K_OCC::trimesh(PyObject* self, PyObject* args)
 
   E_Int aniso = false;
 
-  if (hausd < 0 && hmax > 0 && hmin > 0)
+  if ( (hausd < 0 && hmax > 0) ||
+       (hausd > 0 && std::abs(hmax-hmin) < 1.e-12) ) // iso hmax
   {
     // mode pure hmax
     E_Float dx = (hmax-hmin)/hmax;
@@ -119,7 +120,7 @@ PyObject* K_OCC::trimesh(PyObject* self, PyObject* args)
     if (dx > 0.2) mode.growth_ratio = 1.1;
     //printf("trimesh uniform hmin=" SF_F_ " hmax=" SF_F_ " grading=" SF_F_ "\n", mode.hmin, mode.hmax, mode.growth_ratio);      
   }
-  else if (hausd > 0 && hmax > 0 && hmin > 0 && aniso == true)
+  else if (hausd > 0 && hmax > 0 && hmin >= 0 && aniso == true) // aniso mix
   {
     // mode pure hausd
     mode.metric_mode = mode.ANISO; //ISO_RHO impose la courbure minimum dans les deux directions
@@ -130,7 +131,7 @@ PyObject* K_OCC::trimesh(PyObject* self, PyObject* args)
     mode.nb_smooth_iter = 2; // iter de lissage pour assurer le grading
     mode.symmetrize = true;
   }
-  else if (hausd > 0 && hmax > 0 && hmin > 0 && aniso == false)
+  else if (hausd > 0 && hmax > 0 && hmin >= 0 && aniso == false) // iso rho mix
   {
     // mode mix hmin/hmax/hausd
     mode.metric_mode = mode.ISO_RHO; //ISO_RHO impose la courbure minimum dans les deux directions;
@@ -149,7 +150,10 @@ PyObject* K_OCC::trimesh(PyObject* self, PyObject* args)
     return NULL;
   }
 
-  printf("trimesh hmin=%g hmax=%g hausd=%g sym=%d grading=%g  smooth=" SF_D_ "\n", 
+  if (mode.metric_mode == mode.ISO_CST) printf("trimesh mode=ISO_CST ");
+  else if (mode.metric_mode == mode.ISO_RHO) printf("trimesh mode=ISO_RHO ");
+  else printf("trimesh mode=ANISO ");
+  printf("hmin=%g hmax=%g hausd=%g sym=%d grading=%g  smooth=" SF_D_ "\n", 
           mode.hmin, mode.hmax, mode.chordal_error, 
           mode.symmetrize, mode.growth_ratio, mode.nb_smooth_iter);
 
diff --git a/Cassiopee/OCC/OCC/PyTree.py b/Cassiopee/OCC/OCC/PyTree.py
index 2adfb119f..ebb2ecb2a 100644
--- a/Cassiopee/OCC/OCC/PyTree.py
+++ b/Cassiopee/OCC/OCC/PyTree.py
@@ -652,7 +652,7 @@ def _remeshTreeFromEdges(hook, t, edges):
     for edge in edges:
         edgeno = getNo(edge)
         aedge = C.getFields([Internal.__GridCoordinates__, Internal.__FlowSolutionNodes__], edge, api=2)[0]
-        e = occ.meshOneEdge(hook, edgeno, -1, -1, -1, aedge)
+        e = occ.meshOneEdge(hook, edgeno, -1, -1, -1, -1, aedge)
         dedges[edgeno-1] = e
         cad = Internal.getNodeFromName1(edge, 'CAD')
         render = Internal.getNodeFromName1(edge, '.RenderInfo')
@@ -778,7 +778,7 @@ def _remeshAllEdgesOdd(hook, t):
             D._getCurvilinearAbscissa(edge)
             edgeno = getNo(edge)
             aedge = C.getFields([Internal.__GridCoordinates__, Internal.__FlowSolutionNodes__], edge, api=2)[0]
-            e = occ.meshOneEdge(hook, edgeno, -1, -1, -1, aedge)
+            e = occ.meshOneEdge(hook, edgeno, -1, -1, -1, -1, aedge)
             cad = Internal.getNodeFromName1(edge, 'CAD')
             render = Internal.getNodeFromName1(edge, '.RenderInfo')
             z = Internal.createZoneNode('edge%03d'%(edgeno), e, [],
@@ -795,7 +795,7 @@ def getCADcontainer(t):
     CAD = Internal.getNodeFromName1(t, 'CAD')
     if CAD is None: return [hmin, hmax, hausd]
     hmin = Internal.getNodeFromName1(CAD, 'hmin')
-    if hmin is not None: hmax = Internal.getValue(hmin)
+    if hmin is not None: hmin = Internal.getValue(hmin)
     hmax = Internal.getNodeFromName1(CAD, 'hmax')
     if hmax is not None: hmax = Internal.getValue(hmax)
     hausd = Internal.getNodeFromName1(CAD, 'hausd')
diff --git a/Cassiopee/XCore/XCore/AdaptMesh/AdaptMesh_AdaptGeom.cpp b/Cassiopee/XCore/XCore/AdaptMesh/AdaptMesh_AdaptGeom.cpp
index 08eb133db..8cc650141 100644
--- a/Cassiopee/XCore/XCore/AdaptMesh/AdaptMesh_AdaptGeom.cpp
+++ b/Cassiopee/XCore/XCore/AdaptMesh/AdaptMesh_AdaptGeom.cpp
@@ -792,7 +792,7 @@ PyObject *K_XCORE::AdaptMesh_AdaptGeom(PyObject *self, PyObject *args)
 
         // We need the skin connectivity graph
         S.make_tri_graph();
-        printf("Skin: %lu faces\n", S.tri_graph.nf);
+        printf("Skin: %zu faces\n", S.tri_graph.nf);
 
         // BVH the skin
         skin_fc = (Vec3f *)XCALLOC(S.tri_graph.nf, sizeof(Vec3f));