Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 3267)
@@ -113,5 +113,5 @@
 		if (verbosity>0) printf("Anisotropic mesh adaptation\n");
 		if (verbosity>1) printf("   Processing initial mesh and geometry...\n");
-		Triangles BTh(bamgmesh_in,bamgopts); 
+		Triangles BTh(bamggeom_in,bamgmesh_in,bamgopts); 
 		hmin = Max(hmin,BTh.MinimalHmin());
 		hmax = Min(hmax,BTh.MaximalHmax());
Index: /issm/trunk/src/c/Bamgx/objects/BamgObjects.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/BamgObjects.h	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/BamgObjects.h	(revision 3267)
@@ -2,5 +2,5 @@
 #define _BAMGOBEJECTS_H_
 
-/*Bamg objects (Must be compiled in this order!!*/
+/*Bamg objects (Must be compiled in this order!!)*/
 #include "./Metric.h"
 #include "./DoubleAndInt.h"
@@ -22,4 +22,6 @@
 #include "./Triangles.h"
 #include "./Geometry.h"
+#include "./QuadTree.h"
+#include "./SetOfE4.h"
 
 #endif
Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp	(revision 3267)
@@ -220,6 +220,6 @@
 		register long n0;
 		register QuadTreeBox* b;
-		IntQuad  h=MaxISize,h0;
-		IntQuad  hb=MaxISize;
+		long     h=MaxISize,h0;
+		long     hb=MaxISize;
 		Icoor1   i0=0,j0=0;
 		//iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0)
@@ -356,6 +356,6 @@
 		int l; // level
 		QuadTreeBox * b;
-		IntQuad  h=MaxISize,h0;
-		IntQuad hb =  MaxISize;
+		long     h =MaxISize,h0;
+		long     hb=MaxISize;
 		Icoor1  i0=0,j0=0;
 		Icoor1  iplus( i<MaxISize?(i<0?0:i):MaxISize-1);
Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/QuadTree.h	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/QuadTree.h	(revision 3267)
@@ -10,8 +10,6 @@
 namespace bamg {
 
-	typedef long IntQuad;
-
-	const int MaxDeep = 30;
-	const IntQuad MaxISize = ( 1L << MaxDeep); 
+	const int  MaxDeep = 30;
+	const long MaxISize = ( 1L << MaxDeep); 
 
 	class Triangles;
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp	(revision 3267)
@@ -3,9 +3,8 @@
 #include <cmath>
 #include <ctime>
+#include <assert.h>
 
 #include "BamgObjects.h"
 #include "../shared/shared.h"
-#include "QuadTree.h"
-#include "SetOfE4.h"
 
 #undef __FUNCT__ 
@@ -18,10 +17,30 @@
 
 	/*Constructors/Destructors*/
-	/*FUNCTION Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
-	Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
+	/*FUNCTION Triangles::Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
+	Triangles::Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
+
+		/*Initialize fields*/
 		PreInit(0);
+
+		/*Read background mesh*/
 		ReadMesh(bamgmesh,bamgopts);
+
+		/*Read Geometry*/
+		if(bamggeom->NumEdges==0) {
+			/*Recreate geometry if needed*/
+			printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
+			BuildGeometryFromMesh(bamgopts);
+			Gh.AfterRead();
+		}
+		else{
+			Gh.ReadGeometry(bamggeom,bamgopts);
+			Gh.AfterRead();
+		}
+
+		/*Set integer coordinates*/
 		SetIntCoor();
-		GenerateMeshProperties();
+
+		/*Fill holes and generate mesh properties*/
+		ReconstructExistingMesh();
 	}
 	/*}}}1*/
@@ -32,5 +51,5 @@
 		ReadMesh(index,x,y,nods,nels);
 		SetIntCoor();
-		GenerateMeshProperties();
+		ReconstructExistingMesh();
 	}
 	/*}}}1*/
@@ -127,5 +146,5 @@
 		  Gh.AfterRead(); 
 		  SetIntCoor();
-		  GenerateMeshProperties();
+		  ReconstructExistingMesh();
 
 		  if (!NbSubDomains){
@@ -523,11 +542,4 @@
 		}
 
-		/*Recreate geometry if needed*/
-		if(1!=0) {
-			printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
-			/*Recreate geometry: */
-			BuildGeometryFromMesh(bamgopts);
-			Gh.AfterRead();
-		}
 	}
 	/*}}}1*/
@@ -2604,7 +2616,15 @@
 	}
 	/*}}}1*/
-	/*FUNCTION Triangles::GenerateMeshProperties{{{1*/
-	void Triangles::GenerateMeshProperties() {
+	/*FUNCTION Triangles::ReconstructExistingMesh{{{1*/
+	void Triangles::ReconstructExistingMesh() {
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
+
+		/*This routine reconstruct an existing mesh to make it CONVEX:
+		 * -all the holes are filled
+		 * -concave boundaries are filled
+		 * A convex mesh is required for a lot of operations. This is why every mesh
+		 * goes through this process.
+		 * This routine also generates mesh properties such as adjencies,...
+		 */
 
 		int verbosity=0;
@@ -2614,10 +2634,8 @@
 		// find extrema coordinates of vertices pmin,pmax
 		long i;
-		if(verbosity>2) printf("      Filling holes in mesh of %i vertices\n",nbv); 
+		if(verbosity>2) printf("      Completing Triangles structure for a mesh of %i vertices\n",nbv); 
 
 		//initialize ordre
-		if (!ordre){
-			throw ErrorException(__FUNCT__,exprintf("!ordre"));
-		}
+		assert(ordre);
 		for (i=0;i<nbv;i++) ordre[i]=0;
 
@@ -2651,7 +2669,6 @@
 				}
 				else if(st[k]>=0) {
-					if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){
-						throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))"));
-					}
+					assert(!triangles[i].TriangleAdj(j));
+					assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
 
 					triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
@@ -2681,5 +2698,5 @@
 		for (i=0;i<edge4->nb();i++){
 			if (st[i] >=0){ // edge alone 
-				if (i < nbe) {
+				if (i<nbe){
 					long i0=edge4->i(i);
 					ordre[i0] = vertices+i0;
@@ -2689,12 +2706,16 @@
 				else {
 					k=k+1;
-					if (k <20) {
-						throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));
+					if (k<10) {
+						//print only 10 edges
+						printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
 					}
-				}
-			}
-		}
-		if(k != 0) {
-			throw ErrorException(__FUNCT__,exprintf("%i boundary edges  are not defined as edges",k));
+					else if (k==10){
+						printf("Other lost boundary edges not shown...\n");
+					}
+				}
+			}
+		}
+		if(k) {
+			throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k));
 		}
 
@@ -2705,13 +2726,13 @@
 			vertices[i].vint=0;
 			if (ordre[i]){ 
-				ordre[nbvb++] = ordre[i];
-			}
-		}
-
-		Triangle* savetriangles= triangles;
+				ordre[nbvb++]=ordre[i];
+			}
+		}
+
+		Triangle* savetriangles=triangles;
 		long savenbt=nbt;
 		long savenbtx=nbtx;
-		SubDomain * savesubdomains = subdomains;
-		subdomains = 0;
+		SubDomain* savesubdomains=subdomains;
+		subdomains=0;
 
 		long  Nbtriafillhole = 2*nbvb;
@@ -2723,17 +2744,16 @@
 
 		for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 
-		 if  ( ++i >= nbvb) {
-			 throw ErrorException(__FUNCT__,exprintf("GenerateMeshProperties: All the vertices are aligned"));
+		 if  (++i>=nbvb) {
+			 throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));
 		 }
-		Exchange( ordre[2], ordre[i]);
+		Exchange(ordre[2], ordre[i]);
 
 		Vertex *  v0=ordre[0], *v1=ordre[1];
 
-
-		triangles[0](0) = 0; // sommet pour infini 
+		triangles[0](0) = NULL; // Infinite vertex
 		triangles[0](1) = v0;
 		triangles[0](2) = v1;
 
-		triangles[1](0) = 0;// sommet pour infini 
+		triangles[1](0) = NULL;// Infinite vertex
 		triangles[1](2) = v0;
 		triangles[1](1) = v1;
@@ -2745,6 +2765,6 @@
 		triangles[0].SetAdj2(e2, &triangles[1] ,e1);
 
-		triangles[0].det = -1;  // faux triangles
-		triangles[1].det = -1;  // faux triangles
+		triangles[0].det = -1;  // boundary triangles
+		triangles[1].det = -1;  // boundary triangles
 
 		triangles[0].SetTriangleContainingTheVertex();
@@ -2881,11 +2901,24 @@
 		SetVertexFieldOn();
 
-		for (i=0;i<nbe;i++)
-		 if(edges[i].onGeometry) 
-		  for(int j=0;j<2;j++)
-			if (!edges[i].adj[j])
-			 if(!edges[i][j].onGeometry->IsRequiredVertex()) {
-				 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));
-			 }
+		/*Check requirements consistency*/
+		for (i=0;i<nbe;i++){
+			/*If the current mesh edge is on Geometry*/
+			if(edges[i].onGeometry){
+				for(int j=0;j<2;j++){
+					if (!edges[i].adj[j]){
+						if(!edges[i][j].onGeometry->IsRequiredVertex()) {
+							printf("Error: adj and vertex requires edges [%i %i]=%i on = %i\n",i,j,Number(edges[i][j]),Gh.Number(edges[i].onGeometry));
+							if (edges[i][j].onGeometry->OnGeomVertex())
+							 printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv));
+							else if (edges[i][j].onGeometry->OnGeomEdge())
+							 printf(" Edges %i \n",Gh.Number(edges[i][j].onGeometry->ge));
+							else
+							 printf(" = %p\n",edges[i][j].onGeometry);
+							throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)"));
+						}
+					}
+				}
+			}
+		}
 	}
 	/*}}}1*/
@@ -5271,5 +5304,5 @@
 		//  ReMakeTriangleContainingTheVertex();
 
-		GenerateMeshProperties();
+		ReconstructExistingMesh();
 
 		if (verbosity>2){
Index: /issm/trunk/src/c/Bamgx/objects/Triangles.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/Triangles.h	(revision 3267)
@@ -68,5 +68,5 @@
 
 			//Constructors/Destructors
-			Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+			Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
 			Triangles(double* index,double* x,double* y,int nods,int nels);
 			Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,long nbvxx=0 ); //copy operator
@@ -145,5 +145,5 @@
 			int  UnCrack();
 			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
-			void GenerateMeshProperties() ;
+			void ReconstructExistingMesh();
 			int  CrackMesh();
 
Index: /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3266)
+++ /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h	(revision 3267)
@@ -24,6 +24,6 @@
 			double abscisse;  
 			union{ 
-				GeometricalVertex * gv; // if abscisse <0; 
-				GeometricalEdge * ge;  // if abscisse in [0..1]
+				GeometricalVertex* gv; // if abscisse <0; 
+				GeometricalEdge*   ge; // if abscisse in [0..1]
 			};
 
@@ -46,5 +46,4 @@
 
 			//Inline methods
-			//inline void Set(const Triangles &,long,Triangles &);
 			void Set(const VertexOnGeom&,const Triangles &,Triangles &);