Index: /issm/trunk-jpl/src/c/bamg/Geometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Geometry.cpp	(revision 22379)
+++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp	(revision 22380)
@@ -438,5 +438,5 @@
 		return c - curves;
 	}/*}}}*/
-	void Geometry::PostRead(){/*{{{*/
+	void Geometry::PostRead(bool checkcurve){/*{{{*/
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
 
@@ -606,17 +606,44 @@
 			}
 
-			//all vertices provided in geometry are corners (ord = number of edges holding i)
-			vertices[i].SetCorner() ; 
-			if(ord==2){ 
-				long  n1 = head_v[i];
-				long  n2 = next_p[n1];
-				long  i1 = n1/2, i2 = n2/2; // edge number
-				long  j1 = n1%2, j2 = n2%2; // vertex in the edge
-				/* if the edge type/referencenumber a changing then is SetRequired();*/
-				if (edges[i1].type != edges[i2].type || edges[i1].Required()){
-					vertices[i].SetRequired();
-				}
-				if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
-					vertices[i].SetRequired();
+			/*Do we want to check for curve? Default is no, but if we are reconstructing, it's better to turn it on with a small threshold*/
+			if(checkcurve){
+				/*angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) */
+				if(ord==2){
+					IssmDouble MaxCornerAngle = 1*Pi/180; /*default is 1 degree*/
+					long  n1 = head_v[i];
+					long  n2 = next_p[n1];
+					long  i1 = n1/2, i2 = n2/2; // edge number
+					long  j1 = n1%2, j2 = n2%2; // vertex in the edge
+					float angle1=  j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
+					float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
+					float da12 = Abs(angle2-angle1);
+					if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
+						vertices[i].SetCorner() ;
+					}
+					/* if the edge type/referencenumber a changing then is SetRequired();*/
+					if(edges[i1].type != edges[i2].type || edges[i1].Required()){
+						vertices[i].SetRequired();
+					}
+					if(edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
+						vertices[i].SetRequired();
+					}
+				}
+				if(ord!=2) vertices[i].SetCorner();
+			}
+			else{
+				/*all vertices provided in geometry are corners (ord = number of edges holding i)*/
+				vertices[i].SetCorner() ; 
+				if(ord==2){ 
+					long  n1 = head_v[i];
+					long  n2 = next_p[n1];
+					long  i1 = n1/2, i2 = n2/2; // edge number
+					long  j1 = n1%2, j2 = n2%2; // vertex in the edge
+					/* if the edge type/referencenumber a changing then is SetRequired();*/
+					if (edges[i1].type != edges[i2].type || edges[i1].Required()){
+						vertices[i].SetRequired();
+					}
+					if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
+						vertices[i].SetRequired();
+					}
 				}
 			}
Index: /issm/trunk-jpl/src/c/bamg/Geometry.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Geometry.h	(revision 22379)
+++ /issm/trunk-jpl/src/c/bamg/Geometry.h	(revision 22380)
@@ -52,5 +52,5 @@
 			void             ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
 			void             Init(void);
-			void             PostRead();
+			void             PostRead(bool checkcurve=false);
 			long             GetId(const GeomVertex &t) const;
 			long             GetId(const GeomVertex *t) const;
Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 22379)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 22380)
@@ -30,5 +30,5 @@
 			_printf_("WARNING: mesh present but no geometry found. Reconstructing...\n");
 			BuildGeometryFromMesh(bamgopts);
-			Gh.PostRead();
+			Gh.PostRead(true);
 		}
 
@@ -1215,16 +1215,16 @@
 		/*Construction of the edges*/
 
-		//initialize st and edge4
+		/*initialize st and edge4*/
 		SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
 		long*        st   = new long[nbt*3];
 
-		//initialize st as -1 (chaining algorithm)
+		/*initialize st as -1 (chaining algorithm)*/
 		for (i=0;i<nbt*3;i++) st[i]=-1;
 
-		//build edge4 (chain)
+		/*build edge4 (chain)*/
 		for (i=0;i<nbe;i++){
 			edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1]));
 		}
-		//check that there is no double edge
+		/*check that there is no double edge*/
 		if (nbe !=  edge4->nb()){ 
 			delete [] st;
@@ -1234,5 +1234,5 @@
 		long nbeold = nbe;
 
-		//Go through the triangles and ass the edges in edge4 if they are not there yet
+		//Go through the triangles and add the edges in edge4 if they are not there yet
 		for (i=0;i<nbt;i++){
 			//3 edges per triangle
@@ -1251,16 +1251,16 @@
 						_error_("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)");
 					}
-					//OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list
+					/*OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list*/
 					triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
 					if (invisible)  triangles[i].SetHidden(j);
-					// if k < nbe mark the edge as on Boundary (Locked)
+					/* if k < nbe mark the edge as on Boundary (Locked)*/
 					if (k<nbe) {
 						triangles[i].SetLocked(j);
 					}
-					//set st[k] as negative so that the edge should not be called again
+					/*set st[k] as negative so that the edge should not be called again*/
 					st[k]=-2-st[k]; 
 				}
-				//else (see 3 lines above), the edge has been called more than twice: return error
 				else {
+					/*else (see 3 lines above), the edge has been called more than twice: return error*/
 					_printf_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")\n");
 					_printf_("Edge " << j << " of triangle " << i << "\n");
@@ -1286,5 +1286,5 @@
 		}
 
-		// check consistency of edge[].adj and geometrical required  vertices
+		/*check consistency of edge[].adj and geometrical required  vertices*/
 		k=0; kk=0;
 		for (i=0;i<nbedges;i++){
@@ -1304,5 +1304,4 @@
 
 		/*Constructions of edges*/
-
 		k += kk;
 		kk=0;