Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2795)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2796)
@@ -109,5 +109,8 @@
 
 		//Build output
+		if (verbosity>1) printf("   Write Mesh...\n");
 		Th.WriteBamgMesh(bamgmesh,bamgopts);
+		if (verbosity>1) printf("   Write Geometry...\n");
+		Gh.WriteGeometry(bamggeom,bamgopts);
 		/*}}}*/
 	}
@@ -164,5 +167,8 @@
 
 		//Build output
+		if (verbosity>1) printf("   Write Mesh...\n");
 		Th.WriteBamgMesh(bamgmesh,bamgopts);
+		if (verbosity>1) printf("   Write Geometry...\n");
+		Th.Gh.WriteGeometry(bamggeom,bamgopts);
 
 		/*clean up*/
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2795)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2796)
@@ -990,20 +990,13 @@
   Int4 Number(const GeometricalEdge * t) const  { return t - edges;}
   Int4 Number(const Curve * c) const  { return c - curves;}
-  
+
   void UnMarkEdges() {
-    for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
-
- GeometricalEdge *  ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
+	  for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
+
+  GeometricalEdge *  ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
   GeometricalEdge *  Contening(const R2 P,  GeometricalEdge * start) const;
- friend std::ostream& operator <<(std::ostream& f, const   Geometry & Gh); 
- void Write(const char * filename);
- 
-#ifdef DEBUG
-  void inline Check();   
-#endif
-#ifdef DRAWING
-  void  Draw() const ;
-  void  InitDraw() const ;
-#endif
+  friend std::ostream& operator <<(std::ostream& f, const   Geometry & Gh); 
+  void Write(const char * filename);
+  void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
   
 }; // End Class Geometry
Index: /issm/trunk/src/c/Bamgx/MeshWrite.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshWrite.cpp	(revision 2795)
+++ /issm/trunk/src/c/Bamgx/MeshWrite.cpp	(revision 2796)
@@ -230,4 +230,163 @@
        f << *this;
     }
+}
+
+void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
+
+	int verbose;
+	int nbreq=0;
+	int nbreqv=0;
+	int nbtan=0;
+	int nbcracked=0;
+	int i,count;
+
+	verbose=bamgopts->verbose;
+
+	//Vertices
+	if(verbose>3) printf("      writing Vertices\n");
+	bamggeom->NumVertices=nbv;
+	xfree((void**)&bamggeom->Vertices);
+	if (nbv){
+		bamggeom->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
+		for (i=0;i<nbv;i++){
+			bamggeom->Vertices[i*3+0]=vertices[i].r.x;
+			bamggeom->Vertices[i*3+1]=vertices[i].r.y;
+			bamggeom->Vertices[i*3+2]=vertices[i].ref();
+
+			//update counters
+			if (vertices[i].Required()) nbreqv++;
+		}
+	}
+	else{
+		bamggeom->Vertices=NULL;
+	}
+
+	//Edges
+	if(verbose>3) printf("      writing Edges\n");
+	bamggeom->NumEdges=nbe;
+	xfree((void**)&bamggeom->Edges);
+	if (nbe){
+		bamggeom->Edges=(double*)xmalloc(3*nbe*sizeof(double));
+		for (i=0;i<nbe;i++){
+			bamggeom->Edges[i*3+0]=Number(edges[i][0])+1; //back to Matlab indexing
+			bamggeom->Edges[i*3+1]=Number(edges[i][1])+1; //back to Matlab indexing
+			//bamggeom->Edges[i*3+2]=(double)edges[i].ref(); //TEST does not compile???
+			bamggeom->Edges[i*3+2]=1; //TEST for now
+
+			//update counters
+			if (edges[i].Required()) nbreq++;
+			if (edges[i].Cracked()){
+				if (i<=Number(edges[i].link)) nbcracked++;
+			}
+			if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
+			if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
+		}
+	}
+	else{
+		bamggeom->Edges=NULL;
+	}
+
+	//CrackedEdges
+	if(verbose>3) printf("      writing CrackedEdges\n");
+	bamggeom->NumCrackedEdges=nbcracked;
+	xfree((void**)&bamggeom->CrackedEdges);
+	if (nbcracked){
+		bamggeom->CrackedEdges=(double*)xmalloc(2*nbcracked*sizeof(double));
+		count=0;
+		for (i=0;i<nbe;i++){
+			if (edges[i].Cracked()){
+				if (i<=Number(edges[i].link)){
+					bamggeom->CrackedEdges[count*2+0]=i+1;                     //back to Matlab indexing
+					bamggeom->CrackedEdges[count*2+1]=Number(edges[i].link)+1; //back to Matlab indexing
+					count=count+1;
+				}
+			}
+		}
+	}
+	else{
+		bamggeom->CrackedEdges=NULL;
+	}
+
+	//RequiredEdges
+	if(verbose>3) printf("      writing RequiredEdges\n");
+	bamggeom->NumRequiredEdges=nbreq;
+	xfree((void**)&bamggeom->RequiredEdges);
+	if (nbreq){
+		bamggeom->RequiredEdges=(double*)xmalloc(1*nbreq*sizeof(double));
+		count=0;
+		for (i=0;i<nbe;i++){
+			if (edges[i].Required()){
+				bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
+				count=count+1;
+			}
+		}
+	}
+	else{
+		bamggeom->RequiredEdges=NULL;
+	}
+
+	//No corners
+
+	//RequiredVertices
+	if(verbose>3) printf("      writing RequiredVertices\n");
+	bamggeom->NumRequiredVertices=nbreqv;
+	xfree((void**)&bamggeom->RequiredVertices);
+	if (nbreqv){
+		bamggeom->RequiredVertices=(double*)xmalloc(1*nbreqv*sizeof(double));
+		count=0;
+		for (i=0;i<nbe;i++){
+			if (vertices[i].Required()){
+				bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
+				count=count+1;
+			}
+		}
+	}
+	else{
+		bamggeom->RequiredVertices=NULL;
+	}
+
+	//SubDomains
+	if(verbose>3) printf("      writing SubDomains\n");
+	bamggeom->NumSubDomains=NbSubDomains;
+	xfree((void**)&bamggeom->SubDomains);
+	if (NbSubDomains){
+		bamggeom->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
+		for (i=0;i<NbSubDomains;i++){
+			bamggeom->SubDomains[4*i+0]=2;
+			bamggeom->SubDomains[4*i+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
+			bamggeom->SubDomains[4*i+2]=subdomains[i].sens;
+			bamggeom->SubDomains[4*i+3]=subdomains[i].ref;
+		}
+	}
+	else{
+		bamggeom->SubDomains=NULL;
+	}
+
+	//TangentAtEdges
+	if(verbose>3) printf("      writing TangentAtEdges\n");
+	bamggeom->NumTangentAtEdges=nbtan;
+	xfree((void**)&bamggeom->TangentAtEdges);
+	if (nbtan){
+		bamggeom->TangentAtEdges=(double*)xmalloc(4*nbtan*sizeof(double));
+		count=0;
+		for (i=0;i<nbe;i++){
+			if (edges[i].TgA() && edges[i][0].Corner()){
+				bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
+				bamggeom->TangentAtEdges[4*i+1]=1;
+				bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
+				bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
+			}
+			if (edges[i].TgB() && edges[i][1].Corner()){
+				bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
+				bamggeom->TangentAtEdges[4*i+1]=2;
+				bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
+				bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[2].y;
+			}
+			count=count+1;
+		}
+	}
+	else{
+		bamggeom->TangentAtEdges=NULL;
+	}
 }