Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 5399)
+++ /issm/trunk/src/c/Makefile.am	(revision 5400)
@@ -53,5 +53,5 @@
 					./objects/Bamg/ListofIntersectionTriangles.cpp\
 					./objects/Bamg/ListofIntersectionTriangles.h\
-					./objects/Bamg/MatVVP2x2.cpp\
+					./objects/Bamg/EigenMetric.cpp\
 					./objects/Bamg/Metric.cpp\
 					./objects/Bamg/Metric.h\
@@ -605,5 +605,5 @@
 					./objects/Bamg/ListofIntersectionTriangles.cpp\
 					./objects/Bamg/ListofIntersectionTriangles.h\
-					./objects/Bamg/MatVVP2x2.cpp\
+					./objects/Bamg/EigenMetric.cpp\
 					./objects/Bamg/Metric.cpp\
 					./objects/Bamg/Metric.h\
Index: /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp	(revision 5399)
+++ /issm/trunk/src/c/modules/Bamgx/Bamgx.cpp	(revision 5400)
@@ -53,5 +53,5 @@
 		for(i=0;i<Gh.nbv;i++){
 			Metric M=Gh[i];
-			MatVVP2x2 Vp(M/coef);
+			EigenMetric Vp(M/coef);
 			Vp.Maxh(bamgopts->hmax);
 			Vp.Minh(bamgopts->hmin);
@@ -131,5 +131,5 @@
 				if (!isnan(bamgopts->hminVertices[i])){
 					Metric M=BTh.vertices[i].m;
-					MatVVP2x2 Vp(M/coef);
+					EigenMetric Vp(M/coef);
 					Vp.Minh(bamgopts->hminVertices[i]);
 					BTh.vertices[i].m=Vp;
@@ -144,5 +144,5 @@
 				if (!isnan(bamgopts->hmaxVertices[i])){
 					Metric M=BTh.vertices[i].m;
-					MatVVP2x2 Vp(M/coef);
+					EigenMetric Vp(M/coef);
 					Vp.Minh(bamgopts->hmaxVertices[i]);
 					BTh.vertices[i].m=Vp;
Index: /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/BamgVertex.cpp	(revision 5400)
@@ -83,5 +83,5 @@
 
 		//Get eigen values and vectors of Miv
-		MatVVP2x2 Vp(Miv);
+		EigenMetric Vp(Miv);
 
 		//move eigen valuse to their absolute values
Index: /issm/trunk/src/c/objects/Bamg/Curve.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Curve.cpp	(revision 5400)
@@ -11,6 +11,10 @@
 	/*Constructors/Destructors*/
 	/*FUNCTION Curve::Curve(){{{1*/
-	Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0){
-
+	Curve::Curve(){
+		FirstEdge=NULL;
+		LastEdge=NULL;
+		FirstVertexIndex=0;
+		LastVertexIndex=0;
+		next=NULL;
 	} 
 	/*}}}*/
@@ -20,6 +24,6 @@
 	void Curve::Reverse() {
 		/*reverse the direction of the curve */
-		Exchange(be,ee);
-		Exchange(kb,ke);
+		Exchange(FirstEdge,LastEdge);
+		Exchange(FirstVertexIndex,LastVertexIndex);
 	}
 	/*}}}*/
@@ -27,6 +31,6 @@
 	void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){
 		*this = rec;
-		be = GhNew.edges + Gh.GetId(be);    
-		ee = GhNew.edges + Gh.GetId(ee); 
+		FirstEdge = GhNew.edges + Gh.GetId(FirstEdge);    
+		LastEdge = GhNew.edges + Gh.GetId(LastEdge); 
 		if(next) next= GhNew.curves + Gh.GetId(next); 
 	}
Index: /issm/trunk/src/c/objects/Bamg/Curve.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Curve.h	(revision 5400)
@@ -13,6 +13,8 @@
 	class Curve {
 		public:
-			GeometricalEdge *be,*ee; // begin et end edge
-			int kb,ke;  //  begin vetex and end vertex
+			GeometricalEdge *FirstEdge; //First edge of the curve
+			GeometricalEdge *LastEdge;  //Last edge of the curve
+			int FirstVertexIndex;       //Last vertex index in the last edge
+			int LastVertexIndex;        //First Vertex index in the first edge
 			Curve* next; // next curve equi to this
 
Index: /issm/trunk/src/c/objects/Bamg/EigenMetric.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/EigenMetric.cpp	(revision 5400)
+++ /issm/trunk/src/c/objects/Bamg/EigenMetric.cpp	(revision 5400)
@@ -0,0 +1,166 @@
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <ctime>
+
+#include "Metric.h"
+
+namespace bamg {
+
+	/*Constructor*/
+	/*FUNCTION EigenMetric::EigenMetric(const Metric M){{{1*/
+	EigenMetric::EigenMetric(const Metric M){
+		/*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
+
+		/*Intermediaries*/
+		double a11=M.a11,a21=M.a21,a22=M.a22;
+		double norm1,norm2,normM;
+		double delta,b;
+
+		/*To get the eigen values, we must solve the following equation:
+		 *     | a11 - lambda    a21        |
+		 * det |                            | = 0
+		 *     | a21             a22-lambda |
+		 *
+		 * We have to solve the following polynom:
+		 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
+
+		/*Compute polynom determinant*/
+		b=-a11-a22;
+		delta=b*b - 4*(a11*a22-a21*a21);
+
+
+		/*Compute norm of M to avoid round off errors*/
+		normM=a11*a11 + a22*a22 + a21*a21;
+
+		/*1: normM too small: eigen values = 0*/
+		if(normM<1.e-30){
+			lambda1=0;
+			lambda2=0;
+			v.x=1;
+			v.y=0;
+		}
+		/*2: delta is small -> double root*/
+		else if (delta < 1.e-5*normM){
+			lambda1=-b/2;
+			lambda2=-b/2;
+			v.x=1;
+			v.y=0;
+		}
+		/*3: general case -> two roots*/
+		else{
+			delta     = sqrt(delta);
+			lambda1   = (-b-delta)/2.0;
+			lambda2   = (-b+delta)/2.0;
+
+			/*Now, one must find the eigen vectors. For that we use the following property of the inner product
+			 *    <Ax,y> = <x,tAy>
+			 * Here, M'(M-lambda*Id) is symmetrical, which gives:
+			 *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
+			 * And we have the following:
+			 *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
+			 * We have shown that
+			 *    Im(M') ⊥ Ker(M')
+			 *
+			 * To find the eigen vectors of M, we only have to find two vectors
+			 * of the image of M' and take their perpendicular as long as they are
+			 * not 0.
+			 * To do that, we take the images (1,0) and (0,1):
+			 *  x1 = (a11 - lambda)      x2 = a21
+			 *  y1 = a21                 y2 = (a22-lambda)
+			 *
+			 * We take the vector that has the larger norm and take its perpendicular.*/
+
+			double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
+			double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
+
+			if (norm2<norm1){
+				norm1=sqrt(norm1);
+				v.x = - a21/norm1;
+				v.y = (a11-lambda1)/norm1;
+			}
+			else{
+				norm2=sqrt(norm2);
+				v.x = - (a22-lambda1)/norm2;
+				v.y = a21/norm2;
+			}
+		}
+
+	}
+	/*}}}1*/
+	/*FUNCTION EigenMetric::EigenMetric(double r1,double r2,const D2 vp1){{{1*/
+	EigenMetric::EigenMetric(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){
+
+	}/*}}}*/
+
+	/*Methods*/
+	/*FUNCTION EigenMetric::Abs{{{1*/
+	void   EigenMetric::Abs(){
+		lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Aniso{{{1*/
+	double EigenMetric::Aniso() const  { 
+		return sqrt( Aniso2());
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Aniso2{{{1*/
+	double EigenMetric::Aniso2() const  { 
+		return lmax()/lmin();
+	}/*}}}*/
+	/*FUNCTION EigenMetric::BoundAniso{{{1*/
+	void   EigenMetric::BoundAniso(const double c){ 
+		BoundAniso2(1/(c*c));
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Echo {{{1*/
+	void EigenMetric::Echo(void){
+
+		printf("EigenMetric:\n");
+		printf("   lambda1: %g\n",lambda1);
+		printf("   lambda2: %g\n",lambda2);
+		printf("   v.x: %g\n",v.x);
+		printf("   v.y: %g\n",v.y);
+
+		return;
+	}
+	/*}}}*/
+	/*FUNCTION EigenMetric::hmin{{{1*/
+	double EigenMetric::hmin() const {
+		return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));
+	}/*}}}*/
+	/*FUNCTION EigenMetric::hmax{{{1*/
+	double EigenMetric::hmax() const {
+		return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Isotrope{{{1*/
+	void   EigenMetric::Isotrope() {
+		lambda1=lambda2=bamg::Max(lambda1,lambda2);
+	}/*}}}*/
+	/*FUNCTION EigenMetric::lmax{{{1*/
+	double EigenMetric::lmax() const {
+		return bamg::Max3(lambda1,lambda2,1e-30);
+	}/*}}}*/
+	/*FUNCTION EigenMetric::lmin{{{1*/
+	double EigenMetric::lmin() const {
+		return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Min{{{1*/
+	void   EigenMetric::Min(double a) { 
+		lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Max{{{1*/
+	void   EigenMetric::Max(double a) { 
+		lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Minh{{{1*/
+	void   EigenMetric::Minh(double h) {
+		Min(1.0/(h*h));
+	}/*}}}*/
+	/*FUNCTION EigenMetric::Maxh{{{1*/
+	void   EigenMetric::Maxh(double h) {
+		Max(1.0/(h*h));
+	}/*}}}*/
+	/*FUNCTION EigenMetric::pow{{{1*/
+	void   EigenMetric::pow(double p){
+		lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);
+	}/*}}}*/
+
+} 
Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp	(revision 5400)
@@ -400,4 +400,99 @@
 
 	/*Methods*/
+	/*FUNCTION Geometry::Echo {{{1*/
+	void Geometry::Echo(void){
+
+		printf("Geometry:\n");
+		printf("   nbv  (number of vertices) : %i\n",nbv);
+		printf("   nbe  (number of edges)    : %i\n",nbe);
+		printf("   nbsubdomains: %i\n",nbsubdomains);
+		printf("   nbcurves: %i\n",nbcurves);
+		printf("   vertices: %p\n",vertices);
+		printf("   edges: %p\n",edges);
+		printf("   quadtree: %p\n",quadtree);
+		printf("   subdomains: %p\n",subdomains);
+		printf("   curves: %p\n",curves);
+		printf("   pmin (x,y): (%g %g)\n",pmin.x,pmin.y);
+		printf("   pmax (x,y): (%g %g)\n",pmax.x,pmax.y);
+		printf("   coefIcoor: %g\n",coefIcoor);
+		printf("   MaxCornerAngle: %g\n",MaxCornerAngle);
+
+		return;
+	}
+	/*}}}*/
+	/*FUNCTION Geometry::Init{{{1*/
+	void Geometry::Init(void){
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
+
+		NbRef=0;
+		nbv=0;
+		nbe=0;
+		quadtree=NULL;
+		curves=NULL;
+		edges=NULL;
+		vertices=NULL;
+		nbsubdomains=0;
+		nbcurves=0;
+		subdomains=NULL;
+		MaxCornerAngle = 10*Pi/180; //default is 10 degres
+	}
+	/*}}}1*/
+	/*FUNCTION Geometry::MinimalHmin{{{1*/
+	double Geometry::MinimalHmin() {
+		/* coeffIcoor = (2^30-1)/D
+		 * We cannot go beyond hmin = D/2^30 because of the quadtree
+		 * hmin is therefore approximately 2/coeffIcoor */
+		return 2.0/coefIcoor;
+	}/*}}}*/
+	/*FUNCTION Geometry::MaximalHmax{{{1*/
+	double Geometry::MaximalHmax() {
+		return Max(pmax.x-pmin.x,pmax.y-pmin.y);
+	}/*}}}*/
+	/*FUNCTION Geometry::GetId(const GeometricalVertex &t){{{1*/
+	long Geometry::GetId(const GeometricalVertex & t) const  {
+		return &t - vertices;
+	}/*}}}*/
+	/*FUNCTION Geometry::GetId(const GeometricalVertex * t){{{1*/
+	long Geometry::GetId(const GeometricalVertex * t) const  {
+		return t - vertices;
+	}/*}}}*/
+	/*FUNCTION Geometry::GetId(const GeometricalEdge & t){{{1*/
+	long Geometry::GetId(const GeometricalEdge & t) const  {
+		return &t - edges;
+	}/*}}}*/
+	/*FUNCTION Geometry::GetId(const GeometricalEdge * t){{{1*/
+	long Geometry::GetId(const GeometricalEdge * t) const  {
+		return t - edges;
+	}/*}}}*/
+	/*FUNCTION Geometry::GetId(const Curve * c){{{1*/
+	long Geometry::GetId(const Curve * c) const  {
+		return c - curves;
+	}/*}}}*/
+	/*FUNCTION Geometry::Containing{{{1*/
+	GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
+
+		GeometricalEdge* on =start,* pon=0;
+		// walk with the cos on geometry
+		int counter=0;
+		while(pon != on){  
+			counter++;
+			ISSMASSERT(counter<100);
+			pon = on;
+			R2 A= (*on)[0];
+			R2 B= (*on)[1];
+			R2 AB = B-A;
+			R2 AP = P-A;
+			R2 BP = P-B;
+			if ( (AB,AP) < 0) 
+			 on = on->Adj[0];
+			else if ( (AB,BP)  > 0) 
+			 on = on->Adj[1];
+			else
+			 return on;
+		}
+		return on;
+	}
+	/*}}}1*/
 	/*FUNCTION Geometry::PostRead{{{1*/
 	void Geometry::PostRead(){
@@ -420,9 +515,9 @@
 
 			/*build integer coordinates (non unique)
-			these coordinates are used by the quadtree to group
-			the vertices by groups of 5:
-			All the coordinates are transformed to ]0,1[^2
-			then, the integer coordinates are computed using 
-			the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
+			  these coordinates are used by the quadtree to group
+			  the vertices by groups of 5:
+			  All the coordinates are transformed to ]0,1[^2
+			  then, the integer coordinates are computed using 
+			  the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
 			vertices[i].i=R2ToI2(vertices[i].r); 
 
@@ -669,6 +764,6 @@
 							GeometricalVertex *a=(*e)(k0); // begin 
 							if(curves){
-								curves[nbcurves].be=e;
-								curves[nbcurves].kb=k0;
+								curves[nbcurves].FirstEdge=e;
+								curves[nbcurves].FirstVertexIndex=k0;
 							}
 							int nee=0;
@@ -679,13 +774,12 @@
 								nb_marked_edges++;
 								e->CurveNumber=nbcurves;
-								if(curves){
-									curves[nbcurves].ee=e;
-									curves[nbcurves].ke=k1;
-								}
-
 								GeometricalVertex *b=(*e)(k1);
 
 								//break if we have reached the other end of the curve
 								if (a==b || b->Required()){
+									if(curves){
+										curves[nbcurves].LastEdge=e;
+										curves[nbcurves].LastVertexIndex=k1;
+									}
 									break;
 								}
@@ -706,8 +800,4 @@
 			if(step==0) curves=new Curve[nbcurves];
 		} 
-		for(int i=0;i<nbcurves ;i++){
-			GeometricalEdge * be=curves[i].be, *eqbe=be;
-			//GeometricalEdge * ee=curves[i].ee, *eqee=be;
-		}
 
 		/*clean up*/
@@ -716,99 +806,4 @@
 		delete [] eangle;
 
-	}
-	/*}}}1*/
-	/*FUNCTION Geometry::Echo {{{1*/
-	void Geometry::Echo(void){
-
-		printf("Geometry:\n");
-		printf("   nbv  (number of vertices) : %i\n",nbv);
-		printf("   nbe  (number of edges)    : %i\n",nbe);
-		printf("   nbsubdomains: %i\n",nbsubdomains);
-		printf("   nbcurves: %i\n",nbcurves);
-		printf("   vertices: %p\n",vertices);
-		printf("   edges: %p\n",edges);
-		printf("   quadtree: %p\n",quadtree);
-		printf("   subdomains: %p\n",subdomains);
-		printf("   curves: %p\n",curves);
-		printf("   pmin (x,y): (%g %g)\n",pmin.x,pmin.y);
-		printf("   pmax (x,y): (%g %g)\n",pmax.x,pmax.y);
-		printf("   coefIcoor: %g\n",coefIcoor);
-		printf("   MaxCornerAngle: %g\n",MaxCornerAngle);
-
-		return;
-	}
-	/*}}}*/
-	/*FUNCTION Geometry::Init{{{1*/
-	void Geometry::Init(void){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
-
-		NbRef=0;
-		nbv=0;
-		nbe=0;
-		quadtree=NULL;
-		curves=NULL;
-		edges=NULL;
-		vertices=NULL;
-		nbsubdomains=0;
-		nbcurves=0;
-		subdomains=NULL;
-		MaxCornerAngle = 10*Pi/180; //default is 10 degres
-	}
-	/*}}}1*/
-	/*FUNCTION Geometry::MinimalHmin{{{1*/
-	double Geometry::MinimalHmin() {
-		/* coeffIcoor = (2^30-1)/D
-		 * We cannot go beyond hmin = D/2^30 because of the quadtree
-		 * hmin is therefore approximately 2/coeffIcoor */
-		return 2.0/coefIcoor;
-	}/*}}}*/
-	/*FUNCTION Geometry::MaximalHmax{{{1*/
-	double Geometry::MaximalHmax() {
-		return Max(pmax.x-pmin.x,pmax.y-pmin.y);
-	}/*}}}*/
-	/*FUNCTION Geometry::GetId(const GeometricalVertex &t){{{1*/
-	long Geometry::GetId(const GeometricalVertex & t) const  {
-		return &t - vertices;
-	}/*}}}*/
-	/*FUNCTION Geometry::GetId(const GeometricalVertex * t){{{1*/
-	long Geometry::GetId(const GeometricalVertex * t) const  {
-		return t - vertices;
-	}/*}}}*/
-	/*FUNCTION Geometry::GetId(const GeometricalEdge & t){{{1*/
-	long Geometry::GetId(const GeometricalEdge & t) const  {
-		return &t - edges;
-	}/*}}}*/
-	/*FUNCTION Geometry::GetId(const GeometricalEdge * t){{{1*/
-	long Geometry::GetId(const GeometricalEdge * t) const  {
-		return t - edges;
-	}/*}}}*/
-	/*FUNCTION Geometry::GetId(const Curve * c){{{1*/
-	long Geometry::GetId(const Curve * c) const  {
-		return c - curves;
-	}/*}}}*/
-	/*FUNCTION Geometry::Containing{{{1*/
-	GeometricalEdge* Geometry::Containing(const R2 P,  GeometricalEdge * start) const {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
-
-		GeometricalEdge* on =start,* pon=0;
-		// walk with the cos on geometry
-		int counter=0;
-		while(pon != on){  
-			counter++;
-			ISSMASSERT(counter<100);
-			pon = on;
-			R2 A= (*on)[0];
-			R2 B= (*on)[1];
-			R2 AB = B-A;
-			R2 AP = P-A;
-			R2 BP = P-B;
-			if ( (AB,AP) < 0) 
-			 on = on->Adj[0];
-			else if ( (AB,BP)  > 0) 
-			 on = on->Adj[1];
-			else
-			 return on;
-		}
-		return on;
 	}
 	/*}}}1*/
Index: sm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp	(revision 5399)
+++ 	(revision )
@@ -1,166 +1,0 @@
-#include <cstdio>
-#include <cstring>
-#include <cmath>
-#include <ctime>
-
-#include "Metric.h"
-
-namespace bamg {
-
-	/*Constructor*/
-	/*FUNCTION MatVVP2x2::MatVVP2x2(const Metric M){{{1*/
-	MatVVP2x2::MatVVP2x2(const Metric M){
-		/*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
-
-		/*Intermediaries*/
-		double a11=M.a11,a21=M.a21,a22=M.a22;
-		double norm1,norm2,normM;
-		double delta,b;
-
-		/*To get the eigen values, we must solve the following equation:
-		 *     | a11 - lambda    a21        |
-		 * det |                            | = 0
-		 *     | a21             a22-lambda |
-		 *
-		 * We have to solve the following polynom:
-		 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
-
-		/*Compute polynom determinant*/
-		b=-a11-a22;
-		delta=b*b - 4*(a11*a22-a21*a21);
-
-
-		/*Compute norm of M to avoid round off errors*/
-		normM=a11*a11 + a22*a22 + a21*a21;
-
-		/*1: normM too small: eigen values = 0*/
-		if(normM<1.e-30){
-			lambda1=0;
-			lambda2=0;
-			v.x=1;
-			v.y=0;
-		}
-		/*2: delta is small -> double root*/
-		else if (delta < 1.e-5*normM){
-			lambda1=-b/2;
-			lambda2=-b/2;
-			v.x=1;
-			v.y=0;
-		}
-		/*3: general case -> two roots*/
-		else{
-			delta     = sqrt(delta);
-			lambda1   = (-b-delta)/2.0;
-			lambda2   = (-b+delta)/2.0;
-
-			/*Now, one must find the eigen vectors. For that we use the following property of the inner product
-			 *    <Ax,y> = <x,tAy>
-			 * Here, M'(M-lambda*Id) is symmetrical, which gives:
-			 *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
-			 * And we have the following:
-			 *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
-			 * We have shown that
-			 *    Im(M') ⊥ Ker(M')
-			 *
-			 * To find the eigen vectors of M, we only have to find two vectors
-			 * of the image of M' and take their perpendicular as long as they are
-			 * not 0.
-			 * To do that, we take the images (1,0) and (0,1):
-			 *  x1 = (a11 - lambda)      x2 = a21
-			 *  y1 = a21                 y2 = (a22-lambda)
-			 *
-			 * We take the vector that has the larger norm and take its perpendicular.*/
-
-			double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
-			double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
-
-			if (norm2<norm1){
-				norm1=sqrt(norm1);
-				v.x = - a21/norm1;
-				v.y = (a11-lambda1)/norm1;
-			}
-			else{
-				norm2=sqrt(norm2);
-				v.x = - (a22-lambda1)/norm2;
-				v.y = a21/norm2;
-			}
-		}
-
-	}
-	/*}}}1*/
-	/*FUNCTION MatVVP2x2::MatVVP2x2(double r1,double r2,const D2 vp1){{{1*/
-	MatVVP2x2::MatVVP2x2(double r1,double r2,const D2 vp1): lambda1(r1),lambda2(r2),v(vp1){
-
-	}/*}}}*/
-
-	/*Methods*/
-	/*FUNCTION MatVVP2x2::Abs{{{1*/
-	void   MatVVP2x2::Abs(){
-		lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Aniso{{{1*/
-	double MatVVP2x2::Aniso() const  { 
-		return sqrt( Aniso2());
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Aniso2{{{1*/
-	double MatVVP2x2::Aniso2() const  { 
-		return lmax()/lmin();
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::BoundAniso{{{1*/
-	void   MatVVP2x2::BoundAniso(const double c){ 
-		BoundAniso2(1/(c*c));
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Echo {{{1*/
-	void MatVVP2x2::Echo(void){
-
-		printf("MatVVP2x2:\n");
-		printf("   lambda1: %g\n",lambda1);
-		printf("   lambda2: %g\n",lambda2);
-		printf("   v.x: %g\n",v.x);
-		printf("   v.y: %g\n",v.y);
-
-		return;
-	}
-	/*}}}*/
-	/*FUNCTION MatVVP2x2::hmin{{{1*/
-	double MatVVP2x2::hmin() const {
-		return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::hmax{{{1*/
-	double MatVVP2x2::hmax() const {
-		return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Isotrope{{{1*/
-	void   MatVVP2x2::Isotrope() {
-		lambda1=lambda2=bamg::Max(lambda1,lambda2);
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::lmax{{{1*/
-	double MatVVP2x2::lmax() const {
-		return bamg::Max3(lambda1,lambda2,1e-30);
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::lmin{{{1*/
-	double MatVVP2x2::lmin() const {
-		return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Min{{{1*/
-	void   MatVVP2x2::Min(double a) { 
-		lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Max{{{1*/
-	void   MatVVP2x2::Max(double a) { 
-		lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Minh{{{1*/
-	void   MatVVP2x2::Minh(double h) {
-		Min(1.0/(h*h));
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::Maxh{{{1*/
-	void   MatVVP2x2::Maxh(double h) {
-		Max(1.0/(h*h));
-	}/*}}}*/
-	/*FUNCTION MatVVP2x2::pow{{{1*/
-	void   MatVVP2x2::pow(double p){
-		lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);
-	}/*}}}*/
-
-} 
Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp	(revision 5400)
@@ -957,5 +957,5 @@
 					c=bamgopts->metric[i*3+2];
 					Metric M(a,b,c);
-					MatVVP2x2 Vp(M/coef);
+					EigenMetric Vp(M/coef);
 
 					Vp.Maxh(hmax);
@@ -1028,5 +1028,5 @@
 					ISSMERROR("ht<=0 || hn<=0");
 				}
-				MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
+				EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg);
 				Metric MVp(Vp);
 				edges[i][j].m.IntersectWith(MVp);
@@ -1212,5 +1212,5 @@
 		//loop over all vertices
 		for (int i=0;i<nbv;i++){
-			MatVVP2x2 Vp(vertices[i]);
+			EigenMetric Vp(vertices[i]);
 			double lmax=Vp.lmax();
 			Vp*=Min(lminaniso,lmax)/lmax;
@@ -2109,5 +2109,5 @@
 						// Compute the matrix with abs(eigen value)
 						Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
-						MatVVP2x2 Vp(M);
+						EigenMetric Vp(M);
 						Vp.Abs();
 						M = Vp;
@@ -3935,5 +3935,5 @@
 		 D2xD2 B1B1 = B1K.t()*B1K;
 		 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
-		 MatVVP2x2 VMK(MK);
+		 EigenMetric VMK(MK);
 		 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
 		 double betaK=0;
@@ -3950,5 +3950,5 @@
 			 D2xD2 BMB = BK.t()*M*BK;
 			 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
-			 MatVVP2x2 VM1(M1);
+			 EigenMetric VM1(M1);
 			 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
 			 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
@@ -5278,16 +5278,16 @@
 					
 					//printf("Dealing with curve number %i\n",nc);
-					//printf("edge on geometry is same as GhCurve? %s\n",(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee)?"yes":"no");
-					//if(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee){
-					//	printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb])?"yes":"no");
-					//	printf("Do we have the right extremity? curve last  vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke])?"yes":"no");
+					//printf("edge on geometry is same as GhCurve? %s\n",(ei.GeometricalEdgeHook==Gh.curves[nc].FirstEdge || ei.GeometricalEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no");
+					//if(ei.GeometricalEdgeHook==Gh.curves[nc].FirstEdge || ei.GeometricalEdgeHook==Gh.curves[nc].LastEdge){
+					//	printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no");
+					//	printf("Do we have the right extremity? curve last  vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no");
 					//}
 					//BUG FIX from original bamg
 					/*Check that we are on the same edge and right vertex (0 or 1) */
-					if(ei.GeometricalEdgeHook==Gh.curves[nc].be  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){
+					if(ei.GeometricalEdgeHook==Gh.curves[nc].FirstEdge  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
 						bcurve[nc]=iedge*2+je;
 						bfind++;	
 					}
-					else if ((ei.GeometricalEdgeHook==Gh.curves[nc].ee  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke]) && bcurve[nc]==-1){
+					else if ((ei.GeometricalEdgeHook==Gh.curves[nc].LastEdge  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
 						bcurve[nc]=iedge*2+je;
 						bfind++;	
Index: /issm/trunk/src/c/objects/Bamg/Metric.cpp
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Metric.cpp	(revision 5400)
@@ -34,5 +34,5 @@
 					a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22);
 
-		MatVVP2x2 vab(mab);
+		EigenMetric vab(mab);
 
 		R2 v1(vab.v.x,vab.v.y);
@@ -49,8 +49,8 @@
 	/*FUNCTION Metric::Metric(double  a,const  Metric ma, double  b,const  Metric mb){{{1*/
 	Metric::Metric(double  a,const  Metric ma, double  b,const  Metric mb) { 
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MatVVP2x2)*/
+		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/EigenMetric)*/
 
 		Metric mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
-		MatVVP2x2 vab(mab);
+		EigenMetric vab(mab);
 
 		R2 v1(vab.v.x,vab.v.y);
Index: /issm/trunk/src/c/objects/Bamg/Metric.h
===================================================================
--- /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5399)
+++ /issm/trunk/src/c/objects/Bamg/Metric.h	(revision 5400)
@@ -12,5 +12,5 @@
 
 	class Metric;
-	class MatVVP2x2;
+	class EigenMetric;
 
 	class Metric{
@@ -22,9 +22,9 @@
 
 			//friends
-			friend class MatVVP2x2;
+			friend class EigenMetric;
 
 			//functions
 			Metric();
-			Metric(const MatVVP2x2);
+			Metric(const EigenMetric);
 			Metric(double a);
 			Metric(double a,double b,double c);
@@ -51,5 +51,5 @@
 	};
 
-	class MatVVP2x2{
+	class EigenMetric{
 		public:
 
@@ -62,6 +62,6 @@
 
 			//functions
-			MatVVP2x2(const Metric );
-			MatVVP2x2(double r1,double r2,const D2 vp1);
+			EigenMetric(const Metric );
+			EigenMetric(double r1,double r2,const D2 vp1);
 			void   Echo();
 			void   Abs();
@@ -100,5 +100,5 @@
 
 	//inlines
-	inline void  MatVVP2x2::BoundAniso2(const double coef){
+	inline void  EigenMetric::BoundAniso2(const double coef){
 		if (coef<=1.00000000001){
 			if (lambda1 < lambda2)
@@ -114,5 +114,5 @@
 		}
 	}
-	inline Metric::Metric(const MatVVP2x2 M) {
+	inline Metric::Metric(const EigenMetric M) {
 		double v00=M.v.x*M.v.x;
 		double v11=M.v.y*M.v.y;