Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2805)
+++ /issm/trunk/src/c/Bamgx/Bamgx.cpp	(revision 2806)
@@ -109,5 +109,5 @@
 		//Build output
 		if (verbosity>1) printf("   Write Mesh...\n");
-		Th.WriteBamgMesh(bamgmesh,bamgopts);
+		Th.WriteMesh(bamgmesh,bamgopts);
 		if (verbosity>1) printf("   Write Geometry...\n");
 		Gh.WriteGeometry(bamggeom,bamgopts);
@@ -167,5 +167,5 @@
 		//Build output
 		if (verbosity>1) printf("   Write Mesh...\n");
-		Th.WriteBamgMesh(bamgmesh,bamgopts);
+		Th.WriteMesh(bamgmesh,bamgopts);
 		if (verbosity>1) printf("   Write Geometry...\n");
 		Th.Gh.WriteGeometry(bamggeom,bamgopts);
Index: /issm/trunk/src/c/Bamgx/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Geometry.cpp	(revision 2806)
+++ /issm/trunk/src/c/Bamgx/Geometry.cpp	(revision 2806)
@@ -0,0 +1,435 @@
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <ctime>
+#include "Mesh2.h"
+#include "QuadTree.h"
+#include "SetOfE4.h"
+
+
+namespace bamg {
+
+	static const  Direction NoDirOfSearch=Direction();
+
+	/*Constructors/Destructors*/
+
+	/*IO*/
+	/*FUNCTION Geometry::WriteGeometry{{{1*/
+	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;
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Geometry::ReadGeometry{{{1*/
+	void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
+
+		int verbose;
+		nbiv=nbv=nbvx=0;
+		nbe=nbt=nbtx=0;
+		NbOfCurves=0;
+
+		Real8 Hmin = HUGE_VAL;// the infinie value 
+		Int4 hvertices=0;
+		int i,j,k,n;
+
+		//initialize some variables
+		int Version=1,dim=2;
+		nbv=bamggeom->NumVertices;
+		nbe=bamggeom->NumEdges;
+		nbvx = nbv;
+		nbiv = nbv;
+		verbose=bamgopts->verbose;
+
+		//some checks
+		if (bamggeom->NumVertices<=0 || bamggeom->Vertices==NULL){
+			throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any vertex"));
+		}
+		if (bamggeom->NumEdges<=0 || bamggeom->Edges==NULL){
+			throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any edge"));
+		}
+
+		//Vertices
+		if (bamggeom->Vertices){
+			if(verbose>3) printf("      processing Vertices\n");
+			vertices = new GeometricalVertex[nbvx];
+			for (i=0;i<nbv;i++) {
+				vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
+				vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
+				vertices[i].ReferenceNumber=(Int4)bamggeom->Vertices[i*3+2];
+				vertices[i].DirOfSearch=NoDirOfSearch;
+				vertices[i].color =0;
+				vertices[i].Set();
+			}
+			//find domain extrema for pmin,pmax
+			pmin =  vertices[0].r;
+			pmax =  vertices[0].r;
+			for (i=0;i<nbv;i++) {
+				pmin.x = Min(pmin.x,vertices[i].r.x);
+				pmin.y = Min(pmin.y,vertices[i].r.y);
+				pmax.x = Max(pmax.x,vertices[i].r.x);
+				pmax.y = Max(pmax.y,vertices[i].r.y);
+			}
+			R2 DD05 = (pmax-pmin)*0.05;
+			pmin -=  DD05;
+			pmax +=  DD05;
+			coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
+			if(coefIcoor <=0){
+				throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive"));
+			}
+		}
+		else{
+			throw ErrorException(__FUNCT__,exprintf("No Vertex provided"));
+		}
+
+		//Edges
+		if (bamggeom->Edges){
+			int i1,i2;
+			R2 zero2(0,0);
+			Real4 *len =0;
+
+			if(verbose>3) printf("      processing Edges\n");
+			edges = new GeometricalEdge[nbe];
+
+			if (!hvertices) {
+				len = new Real4[nbv];
+				for(i=0;i<nbv;i++)
+				 len[i]=0;
+			}
+
+			for (i=0;i<nbe;i++){
+				i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
+				i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
+				edges[i].ref=(Int4)bamggeom->Edges[i*3+2];
+				edges[i].v[0]=  vertices + i1;
+				edges[i].v[1]=  vertices + i2;
+				R2 x12 = vertices[i2].r-vertices[i1].r;
+				Real8 l12=sqrt((x12,x12));
+				edges[i].tg[0]=zero2;
+				edges[i].tg[1]=zero2;
+				edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;
+				edges[i].Adj[0] = edges[i].Adj[1] = 0;
+				edges[i].flag = 0;
+				if (!hvertices) {
+					vertices[i1].color++;
+					vertices[i2].color++;
+					len[i1] += l12;
+					len[i2] += l12;
+				}
+
+				Hmin = Min(Hmin,l12);
+			}
+
+			// definition  the default of the given mesh size 
+			if (!hvertices){
+				for (i=0;i<nbv;i++) 
+				 if (vertices[i].color > 0) 
+				  vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
+				 else 
+				  vertices[i].m=  Metric(Hmin);
+				delete [] len;
+			}
+		}
+		else{
+			throw ErrorException(__FUNCT__,exprintf("No edges provided"));
+		}
+
+		//hVertices
+		if(bamggeom->hVertices){
+			if(verbose>3) printf("      processing hVertices\n");
+			for (i=0;i< nbv;i++){
+				vertices[i].m=Metric((Real4)bamggeom->hVertices[i]);
+			}
+		}
+		else{
+			if(verbose>3) printf("      no hVertices found\n");
+		}
+
+		//MetricVertices
+		if(bamggeom->MetricVertices){
+			if(verbose>3) printf("      processing MetricVertices\n");
+			hvertices=1;
+			for (i=0;i< nbv;i++) {
+				vertices[i].m = Metric((Real4)bamggeom->MetricVertices[i*3+0],(Real4)bamggeom->MetricVertices[i*3+1],(Real4)bamggeom->MetricVertices[i*3+2]);
+			}
+		}
+		else{
+			if(verbose>3) printf("      no MetricVertices found\n");
+		}
+
+		//h1h2VpVertices
+		if(bamggeom->h1h2VpVertices){
+			if(verbose>3) printf("      processing h1h2VpVertices\n");
+			Real4 h1,h2,v1,v2;
+			hvertices =1;
+			for (i=0;i< nbv;i++) {
+				h1=(Real4)bamggeom->MetricVertices[i*4+0];
+				h2=(Real4)bamggeom->MetricVertices[i*4+1];
+				v1=(Real4)bamggeom->MetricVertices[i*4+2];
+				v2=(Real4)bamggeom->MetricVertices[i*4+3];
+				vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));
+			}
+		}
+		else{
+			if(verbose>3) printf("      no h1h2VpVertices found\n");
+		}
+
+		//MaximalAngleOfCorner
+		if (bamgopts->MaximalAngleOfCorner){
+			if(verbose>3) printf("      processing MaximalAngleOfCorner\n");
+			MaximalAngleOfCorner=bamgopts->MaximalAngleOfCorner*Pi/180;
+		}
+		else{
+			if(verbose>3) printf("      no MaximalAngleOfCorner found\n");
+		}
+
+		//TangentAtEdges
+		if (bamggeom->TangentAtEdges){
+			if(verbose>3) printf("      processing TangentAtEdges");
+			int n,i,j,k;
+			R2 tg;
+
+			n=bamggeom->NumTangentAtEdges;
+			for (k=0;k<n;k++) {
+				i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
+				j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
+				tg.x=bamggeom->TangentAtEdges[k*4+2];
+				tg.y=bamggeom->TangentAtEdges[k*4+3];
+				if (j!=0 && j!=1){
+					throw ErrorException(__FUNCT__,exprintf("TangentAtEdges second index should be 1 or 2 only"));
+				}
+				edges[i].tg[j] = tg;
+			}
+		}
+		else{
+			if(verbose>3) printf("      no TangentAtEdges found\n");
+		}
+
+		//Corners
+		if(bamggeom->Corners){
+			if(verbose>3) printf("      processing Corners");
+			n=bamggeom->NumCorners;
+			for (i=0;i<n;i++) {     
+				j=(int)bamggeom->Corners[i]-1; //for C indexing
+				if (j>nbv-1 || j<0){
+					throw ErrorException(__FUNCT__,exprintf("Bad corner definition: should in [0 %i]",nbv));
+				}
+				vertices[j].SetCorner();
+				vertices[j].SetRequired();  }
+		}
+		else{
+			if(verbose>3) printf("      no Corners found\n");
+		}
+
+		//RequiredVertices
+		if(bamggeom->RequiredVertices){
+			if(verbose>3) printf("      processing RequiredVertices");
+			n=bamggeom->NumRequiredVertices;
+			for (i=0;i<n;i++) {     
+				j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
+				if (j>nbv-1 || j<0){
+					throw ErrorException(__FUNCT__,exprintf("Bad RequiredVerticess  definition: should in [0 %i]",nbv));
+				}
+				vertices[j].SetRequired();  }
+		}
+		else{
+			if(verbose>3) printf("      no RequiredVertices found\n");
+		}
+
+		//RequiredEdges
+		if(bamggeom->RequiredEdges){
+			if(verbose>3) printf("      processing RequiredEdges");
+			n=bamggeom->NumRequiredEdges;
+			for (i=0;i<n;i++) {     
+				j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
+				if (j>nbe-1 || j<0){
+					throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv));
+				}
+				edges[j].SetRequired();  }
+		}
+		else{
+			if(verbose>3) printf("      no RequiredEdges found\n");
+		}
+
+		//SubDomain
+		if(bamggeom->SubDomains){
+			Int4 i0,i1,i2,i3;
+			if(verbose>3) printf("      processing SubDomains\n");
+			NbSubDomains=bamggeom->NumSubDomains;
+			subdomains = new GeometricalSubDomain[NbSubDomains];
+			for (i=0;i<NbSubDomains;i++) {
+				i0=(int)bamggeom->SubDomains[i*4+0];
+				i1=(int)bamggeom->SubDomains[i*4+1];
+				i2=(int)bamggeom->SubDomains[i*4+2];
+				i3=(int)bamggeom->SubDomains[i*4+3];
+				if (i0!=2) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 2 (for Edges)"));
+				if (i1>nbe || i1<=0) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: second number should in [1 %i] (edge number)",nbe));
+				subdomains[i].edge=edges + (i1-1);
+				subdomains[i].sens = (int) i2;
+				subdomains[i].ref = i3;
+			}
+		}
+		else{
+			if(verbose>3) printf("      no SubDomains found\n");
+		}
+	}
+	/*}}}1*/
+
+} 
Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2805)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2806)
@@ -838,5 +838,5 @@
 
   void ReadFromMatlabMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
-  void WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+  void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
 
   void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
Index: sm/trunk/src/c/Bamgx/MeshRead.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshRead.cpp	(revision 2805)
+++ 	(revision )
@@ -1,527 +1,0 @@
-/*This file is largely inspired by Freefem++
- 
- Freefem++ is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 2.1 of the License, or
- (at your option) any later version.
- 
- Freefem++  is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU Lesser General Public License for more details.
- 
- You should have received a copy of the GNU Lesser General Public License
- along with Freefem++; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- */
-
-#include "../shared/shared.h"
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-#include <time.h>
-
-#include "Mesh2.h"
-#include "QuadTree.h"
-#include "SetOfE4.h"
-#ifdef __MWERKS__
-#pragma optimization_level 2
-//#pragma inline_depth 1
-#endif
-
-
-namespace bamg {
-
-static const  Direction NoDirOfSearch=Direction();
-
-void Triangles::ReadFromMatlabMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){
-
-	int verbose;
-	Real8 Hmin = HUGE_VAL;// the infinie value 
-	Int4 i1,i2,i3,iref;
-	Int4 i,j;
-	Int4 hvertices =0;
-	Int4 ifgeom=0;
-	Metric M1(1);
-	int Version=1,dim=2;
-
-	verbose=bamgopts->verbose;
-
-	nbv=bamgmesh->NumVertices;
-	nbvx=nbv;
-	nbt=bamgmesh->NumTriangles;
-	nbiv=nbvx;
-
-	//Vertices
-	if(bamgmesh->Vertices){
-		if(verbose>3) printf("      processing Vertices\n");
-
-		vertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
-		ordre=(Vertex**)xmalloc(nbv*sizeof(Vertex*));
-
-		for (i=0;i<nbv;i++){
-			vertices[i].r.x=bamgmesh->Vertices[i*3+0];
-			vertices[i].r.y=bamgmesh->Vertices[i*3+1];
-			vertices[i].ReferenceNumber=1;
-			vertices[i].DirOfSearch =NoDirOfSearch;
-			vertices[i].m=M1;
-			vertices[i].color=(Int4)bamgmesh->Vertices[i*3+2];
-		}
-		nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
-	}
-	else{
-		if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Vertices found in the initial mesh"));
-	}
-
-	//Triangles
-	if(bamgmesh->Triangles){
-		if(verbose>3) printf("      processing Triangles\n");
-		triangles =new Triangle[nbtx]; //we cannot allocate only nbt triangles since 
-												 //other triangles will be added for each edge
-		for (i=0;i<nbt;i++){
-			Triangle & t = triangles[i];
-			i1=(Int4)bamgmesh->Triangles[i*4+0]-1; //for C indexing
-			i2=(Int4)bamgmesh->Triangles[i*4+1]-1; //for C indexing
-			i3=(Int4)bamgmesh->Triangles[i*4+2]-1; //for C indexing
-			t=Triangle(this,i1,i2,i3);
-			t.color=(Int4)bamgmesh->Triangles[i*4+3];
-		}
-	}
-	else{
-		if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Triangles found in the initial mesh"));
-	}
-
-	//Quadrilaterals
-	if(bamgmesh->Quadrilaterals){
-		if(verbose>3) printf("      processing Quadrilaterals\n");
-		Int4 i1,i2,i3,i4,iref;
-		triangles =new Triangle[nbt];
-		for (i=0;i<bamgmesh->NumQuadrilaterals;i++){
-			Triangle & t1 = triangles[2*i];
-			Triangle & t2 = triangles[2*i+1];
-			i1=(Int4)bamgmesh->Quadrilaterals[i*4+0]-1; //for C indexing
-			i2=(Int4)bamgmesh->Quadrilaterals[i*4+1]-1; //for C indexing
-			i3=(Int4)bamgmesh->Quadrilaterals[i*4+2]-1; //for C indexing
-			i4=(Int4)bamgmesh->Quadrilaterals[i*4+3]-1; //for C indexing
-			t1=Triangle(this,i1,i2,i3);
-			t2=Triangle(this,i3,i4,i1);
-			t1.color=1; //reference = 1 for all triangles since it has not been specified
-			t2.color=1; //reference = 1 for all triangles since it has not been specified
-			t1.SetHidden(OppositeEdge[1]); // two times  because the adj was not created 
-			t2.SetHidden(OppositeEdge[1]); 
-		}
-	}
-	else{
-		if(verbose>3) printf("      no Quadrilaterals found\n");
-	}
-
-	//hVertices
-	if(bamgmesh->hVertices){
-		if(verbose>3) printf("      processing hVertices\n");
-		hvertices=1;
-		for (i=0;i< nbv;i++){
-			vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
-		}
-	}
-	else{
-		if(verbose>3) printf("      no hVertices found\n");
-	}
-
-	//VerVerticesOnGeometricEdge
-	if(bamgmesh->VerticesOnGeometricEdge){
-		if(verbose>3) printf("      processing VerticesOnGeometricEdge\n");
-		NbVerticesOnGeomEdge=bamgmesh->NumVerticesOnGeometricEdge;
-		VerticesOnGeomEdge= new  VertexOnGeom[NbVerticesOnGeomEdge] ;
-		for (i=0;i<NbVerticesOnGeomEdge;i++){
-			Int4  i1,i2;
-			Real8 s;
-			i1=(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing
-			i2=(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing
-			s =(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+2];
-			VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
-		}
-	}
-	else{
-		if(verbose>3) printf("      no VertexOnGeometricEdge found\n");
-	}
-
-	//Edges
-	if (bamgmesh->Edges){
-		int i1,i2;
-		R2 zero2(0,0);
-		Real4 *len =0;
-
-		if(verbose>3) printf("      processing Edges\n");
-		nbe=bamgmesh->NumEdges;
-		edges = new Edge[nbe];
-
-		if (!hvertices) {
-			len = new Real4[nbv];
-			for(i=0;i<nbv;i++)
-			 len[i]=0;
-		}
-
-		for (i=0;i<nbe;i++){
-			i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
-			i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
-			edges[i].v[0]= vertices +i1;
-			edges[i].v[1]= vertices +i2;
-			edges[i].adj[0]=0;
-			edges[i].adj[1]=0;
-			R2 x12 = vertices[i2].r-vertices[i1].r;
-			Real8 l12=sqrt( (x12,x12));        
-
-			if (!hvertices) {
-				vertices[i1].color++;
-				vertices[i2].color++;
-				len[i1]+= l12;
-				len[i2] += l12;
-			}
-			Hmin = Min(Hmin,l12);
-		}
-
-		// definition  the default of the given mesh size 
-		if (!hvertices){
-			for (i=0;i<nbv;i++) 
-			 if (vertices[i].color > 0) 
-			  vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
-			 else 
-			  vertices[i].m=  Metric(Hmin);
-			delete [] len;
-		}
-
-		// construction of edges[].adj 
-		for (i=0;i<nbv;i++){ 
-		 vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
-		}
-		for (i=0;i<nbe;i++){
-			for (j=0;j<2;j++) { 
-				Vertex *v=edges[i].v[j];
-				Int4 i0=v->color,j0;
-				if(i0==-1){
-					v->color=i*2+j;
-				}
-				else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
-					j0 =  i0%2;
-					i0 =  i0/2;
-					assert( v ==  edges[i0 ].v[j0]);
-					edges[i ].adj[ j ] =edges +i0;
-					edges[i0].adj[ j0] =edges +i ;
-					assert(edges[i0].v[j0] == v);
-					v->color = -3;
-				}
-			}
-		}
-	}
-	else{
-		if(verbose>3) printf("      no Edges found\n");
-	}
-
-	//EdgeOnGeometricEdge
-	if(bamgmesh->EdgesOnGeometricEdge){
-		if(verbose>3) printf("      processing EdgesOnGeometricEdge\n");
-		int i1,i2,i,j;
-		i2=bamgmesh->NumEdgesOnGeometricEdge;
-		for (i1=0;i1<i2;i1++) {
-			i=(int)bamgmesh->EdgesOnGeometricEdge[i*2+0];
-			j=(int)bamgmesh->EdgesOnGeometricEdge[i*2+1];
-			if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe)) {
-				throw ErrorException(__FUNCT__,exprintf("EdgesOnGeometricEdge error: We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe)"));
-			}
-			edges[i-1].on = Gh.edges + j-1;
-		}
-	}
-	else{
-		if(verbose>3) printf("      no EdgesOnGeometricEdge found\n");
-	}
-
-	//SubDomain
-	if(bamgmesh->SubDomains){
-		Int4 i3,head,sens;
-		if(verbose>3) printf("      processing SubDomains\n");
-		NbSubDomains=bamgmesh->NumSubDomains;
-		subdomains = new SubDomain [ NbSubDomains ];
-		for (i=0;i<NbSubDomains;i++) {
-			i3  =(int)bamgmesh->SubDomains[i*3+0];
-			head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
-			sens=(int)bamgmesh->SubDomains[i*3+2];
-			if (i3!=23) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 3"));
-			if (head<0 || head>=nbt) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: head should in [1 %i] (triangle number)",nbt));
-			subdomains[i].head = triangles+head;
-		}
-	}
-	else{
-		if(verbose>3) printf("      no SubDomains found\n");
-	}
-
-	/*Recreate geometry if needed*/
-	if(1!=0) {
-		printf("Warning: mesh present but no geometry found. Reconstructing...\n");
-		/*Recreate geometry: */
-		ConsGeometry(bamgopts->MaximalAngleOfCorner*Pi/180);
-		Gh.AfterRead();
-	}
-}
-
-Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
-
-  PreInit(0,"none");
-  OnDisk = 1;
-
-  ReadFromMatlabMesh(bamgmesh,bamgopts);
-  SetIntCoor();
-  FillHoleInMesh();
-}
-
-void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
-
-	int verbose;
-	nbiv=nbv=nbvx=0;
-	nbe=nbt=nbtx=0;
-	NbOfCurves=0;
-
-	Real8 Hmin = HUGE_VAL;// the infinie value 
-	Int4 hvertices=0;
-	int i,j,k,n;
-
-	//initialize some variables
-	int Version=1,dim=2;
-	nbv=bamggeom->NumVertices;
-	nbe=bamggeom->NumEdges;
-	nbvx = nbv;
-	nbiv = nbv;
-	verbose=bamgopts->verbose;
-
-	//some checks
-	if (bamggeom->NumVertices<=0 || bamggeom->Vertices==NULL){
-		throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any vertex"));
-	}
-	if (bamggeom->NumEdges<=0 || bamggeom->Edges==NULL){
-		throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any edge"));
-	}
-
-	//Vertices
-	if (bamggeom->Vertices){
-		if(verbose>3) printf("      processing Vertices\n");
-		vertices = new GeometricalVertex[nbvx];
-		for (i=0;i<nbv;i++) {
-			vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
-			vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
-			vertices[i].ReferenceNumber=(Int4)bamggeom->Vertices[i*3+2];
-			vertices[i].DirOfSearch=NoDirOfSearch;
-			vertices[i].color =0;
-			vertices[i].Set();
-		}
-		//find domain extrema for pmin,pmax
-		pmin =  vertices[0].r;
-		pmax =  vertices[0].r;
-		for (i=0;i<nbv;i++) {
-			pmin.x = Min(pmin.x,vertices[i].r.x);
-			pmin.y = Min(pmin.y,vertices[i].r.y);
-			pmax.x = Max(pmax.x,vertices[i].r.x);
-			pmax.y = Max(pmax.y,vertices[i].r.y);
-		}
-		R2 DD05 = (pmax-pmin)*0.05;
-		pmin -=  DD05;
-		pmax +=  DD05;
-		coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
-		if(coefIcoor <=0){
-			throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive"));
-		}
-	}
-	else{
-		throw ErrorException(__FUNCT__,exprintf("No Vertex provided"));
-	}
-
-	//Edges
-	if (bamggeom->Edges){
-		int i1,i2;
-		R2 zero2(0,0);
-		Real4 *len =0;
-
-		if(verbose>3) printf("      processing Edges\n");
-		edges = new GeometricalEdge[nbe];
-
-		if (!hvertices) {
-			len = new Real4[nbv];
-			for(i=0;i<nbv;i++)
-			 len[i]=0;
-		}
-
-		for (i=0;i<nbe;i++){
-			i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
-			i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
-			edges[i].ref=(Int4)bamggeom->Edges[i*3+2];
-			edges[i].v[0]=  vertices + i1;
-			edges[i].v[1]=  vertices + i2;
-			R2 x12 = vertices[i2].r-vertices[i1].r;
-			Real8 l12=sqrt((x12,x12));
-			edges[i].tg[0]=zero2;
-			edges[i].tg[1]=zero2;
-			edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;
-			edges[i].Adj[0] = edges[i].Adj[1] = 0;
-			edges[i].flag = 0;
-			if (!hvertices) {
-				vertices[i1].color++;
-				vertices[i2].color++;
-				len[i1] += l12;
-				len[i2] += l12;
-			}
-
-			Hmin = Min(Hmin,l12);
-		}
-
-		// definition  the default of the given mesh size 
-		if (!hvertices){
-			for (i=0;i<nbv;i++) 
-			 if (vertices[i].color > 0) 
-			  vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
-			 else 
-			  vertices[i].m=  Metric(Hmin);
-			delete [] len;
-		}
-	}
-	else{
-		throw ErrorException(__FUNCT__,exprintf("No edges provided"));
-	}
-
-	//hVertices
-	if(bamggeom->hVertices){
-		if(verbose>3) printf("      processing hVertices\n");
-		for (i=0;i< nbv;i++){
-			vertices[i].m=Metric((Real4)bamggeom->hVertices[i]);
-		}
-	}
-	else{
-		if(verbose>3) printf("      no hVertices found\n");
-	}
-
-	//MetricVertices
-	if(bamggeom->MetricVertices){
-		if(verbose>3) printf("      processing MetricVertices\n");
-		hvertices=1;
-		for (i=0;i< nbv;i++) {
-			vertices[i].m = Metric((Real4)bamggeom->MetricVertices[i*3+0],(Real4)bamggeom->MetricVertices[i*3+1],(Real4)bamggeom->MetricVertices[i*3+2]);
-		}
-	}
-	else{
-		if(verbose>3) printf("      no MetricVertices found\n");
-	}
-
-	//h1h2VpVertices
-	if(bamggeom->h1h2VpVertices){
-		if(verbose>3) printf("      processing h1h2VpVertices\n");
-		Real4 h1,h2,v1,v2;
-		hvertices =1;
-		for (i=0;i< nbv;i++) {
-			h1=(Real4)bamggeom->MetricVertices[i*4+0];
-			h2=(Real4)bamggeom->MetricVertices[i*4+1];
-			v1=(Real4)bamggeom->MetricVertices[i*4+2];
-			v2=(Real4)bamggeom->MetricVertices[i*4+3];
-			vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));
-		}
-	}
-	else{
-		if(verbose>3) printf("      no h1h2VpVertices found\n");
-	}
-
-	//MaximalAngleOfCorner
-	if (bamgopts->MaximalAngleOfCorner){
-		if(verbose>3) printf("      processing MaximalAngleOfCorner\n");
-		MaximalAngleOfCorner=bamgopts->MaximalAngleOfCorner*Pi/180;
-	}
-	else{
-		if(verbose>3) printf("      no MaximalAngleOfCorner found\n");
-	}
-
-	//TangentAtEdges
-	if (bamggeom->TangentAtEdges){
-		if(verbose>3) printf("      processing TangentAtEdges");
-		int n,i,j,k;
-		R2 tg;
-
-		n=bamggeom->NumTangentAtEdges;
-		for (k=0;k<n;k++) {
-			i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
-			j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
-			tg.x=bamggeom->TangentAtEdges[k*4+2];
-			tg.y=bamggeom->TangentAtEdges[k*4+3];
-			if (j!=0 && j!=1){
-				throw ErrorException(__FUNCT__,exprintf("TangentAtEdges second index should be 1 or 2 only"));
-			}
-			edges[i].tg[j] = tg;
-		}
-	}
-	else{
-		if(verbose>3) printf("      no TangentAtEdges found\n");
-	}
-
-	//Corners
-	if(bamggeom->Corners){
-		if(verbose>3) printf("      processing Corners");
-		n=bamggeom->NumCorners;
-		for (i=0;i<n;i++) {     
-			j=(int)bamggeom->Corners[i]-1; //for C indexing
-			if (j>nbv-1 || j<0){
-				throw ErrorException(__FUNCT__,exprintf("Bad corner definition: should in [0 %i]",nbv));
-			}
-			vertices[j].SetCorner();
-			vertices[j].SetRequired();  }
-	}
-	else{
-		if(verbose>3) printf("      no Corners found\n");
-	}
-
-	//RequiredVertices
-	if(bamggeom->RequiredVertices){
-		if(verbose>3) printf("      processing RequiredVertices");
-		n=bamggeom->NumRequiredVertices;
-		for (i=0;i<n;i++) {     
-			j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
-			if (j>nbv-1 || j<0){
-				throw ErrorException(__FUNCT__,exprintf("Bad RequiredVerticess  definition: should in [0 %i]",nbv));
-			}
-			vertices[j].SetRequired();  }
-	}
-	else{
-		if(verbose>3) printf("      no RequiredVertices found\n");
-	}
-
-	//RequiredEdges
-	if(bamggeom->RequiredEdges){
-		if(verbose>3) printf("      processing RequiredEdges");
-		n=bamggeom->NumRequiredEdges;
-		for (i=0;i<n;i++) {     
-			j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
-			if (j>nbe-1 || j<0){
-				throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv));
-			}
-			edges[j].SetRequired();  }
-	}
-	else{
-		if(verbose>3) printf("      no RequiredEdges found\n");
-	}
-
-	//SubDomain
-	if(bamggeom->SubDomains){
-		Int4 i0,i1,i2,i3;
-		if(verbose>3) printf("      processing SubDomains\n");
-		NbSubDomains=bamggeom->NumSubDomains;
-		subdomains = new GeometricalSubDomain[NbSubDomains];
-		for (i=0;i<NbSubDomains;i++) {
-			i0=(int)bamggeom->SubDomains[i*4+0];
-			i1=(int)bamggeom->SubDomains[i*4+1];
-			i2=(int)bamggeom->SubDomains[i*4+2];
-			i3=(int)bamggeom->SubDomains[i*4+3];
-			if (i0!=2) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 2 (for Edges)"));
-			if (i1>nbe || i1<=0) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: second number should in [1 %i] (edge number)",nbe));
-			subdomains[i].edge=edges + (i1-1);
-			subdomains[i].sens = (int) i2;
-			subdomains[i].ref = i3;
-		}
-	}
-	else{
-		if(verbose>3) printf("      no SubDomains found\n");
-	}
-}
-}  // end of namespace bamg 
Index: sm/trunk/src/c/Bamgx/MeshWrite.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshWrite.cpp	(revision 2805)
+++ 	(revision )
@@ -1,375 +1,0 @@
-#include "../shared/shared.h"
-#include "../include/macros.h"
-#include "../toolkits/toolkits.h"
-
-#include <cstdio>
-#include <cstring>
-#include <cmath>
-#include <ctime>
-#include "Mesh2.h"
-#include "QuadTree.h"
-#include "SetOfE4.h"
-
-
-namespace bamg {
-
-void Triangles::WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
-
-	int i,j;
-	int verbose;
-
-	verbose=bamgopts->verbose;
-	Int4 *reft = new Int4[nbt];
-	Int4 nbInT = ConsRefTriangle(reft);
-
-	//Vertices
-	if(verbose>3) printf("      writing Vertices\n");
-	bamgmesh->NumVertices=nbv;
-	xfree((void**)&bamgmesh->Vertices);
-	if (nbv){
-		bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
-		for (i=0;i<nbv;i++){
-			bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
-			bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
-			bamgmesh->Vertices[i*3+2]=vertices[i].ref();
-		}
-	}
-	else{
-		bamgmesh->Vertices=NULL;
-	}
-
-	//Edges
-	if(verbose>3) printf("      writing Edges\n");
-	bamgmesh->NumEdges=nbe;
-	xfree((void**)&bamgmesh->Edges);
-	if (nbe){
-		bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double));
-		for (i=0;i<nbe;i++){
-			bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing
-			bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
-			bamgmesh->Edges[i*3+2]=edges[i].ref;
-		}
-	}
-	else{
-		bamgmesh->Edges=NULL;
-	}
-
-	//CrackedEdges
-	if(verbose>3) printf("      writing CrackedEdges\n");
-	bamgmesh->NumCrackedEdges=NbCrackedEdges;
-	xfree((void**)&bamgmesh->CrackedEdges);
-	if (NbCrackedEdges){
-		bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double));
-		for (i=0;i<NbCrackedEdges;i++){
-			bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing
-			bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing
-		}
-	}
-	else{
-		bamgmesh->CrackedEdges=NULL;
-	}
-
-	//Triangles
-	if(verbose>3) printf("      writing Triangles\n");
-	Int4 k=nbInT-NbOfQuad*2;
-	Int4 num=0;
-	bamgmesh->NumTriangles=k;
-	xfree((void**)&bamgmesh->Triangles);
-	if (k){
-		bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double));
-		for (i=0;i<nbt;i++){
-			Triangle &t=triangles[i];
-			if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
-				bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing
-				bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing
-				bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing
-				bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
-				num=num+1;
-			}
-		}
-	}
-	else{
-		bamgmesh->Triangles=NULL;
-	}
-
-	//Quadrilaterals
-	if(verbose>3) printf("      writing Quadrilaterals\n");
-	bamgmesh->NumQuadrilaterals=NbOfQuad;
-	xfree((void**)&bamgmesh->Quadrilaterals);
-	if (NbOfQuad){
-		bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));
-		for (i=0;i<nbt;i++){
-			Triangle &t =triangles[i];
-			Triangle* ta;
-			Vertex *v0,*v1,*v2,*v3;
-			if (reft[i]<0) continue;
-			if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 
-				bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing
-				bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing
-				bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing
-				bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing
-				bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;
-			}
-		}
-	}
-	else{
-		bamgmesh->Quadrilaterals=NULL;
-	}
-
-	//SubDomains
-	if(verbose>3) printf("      writing SubDomains\n");
-	bamgmesh->NumSubDomains=NbSubDomains;
-	xfree((void**)&bamgmesh->SubDomains);
-	if (NbSubDomains){
-		bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
-		for (i=0;i<NbSubDomains;i++){
-			bamgmesh->SubDomains[i*4+0]=3;
-			bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];
-			bamgmesh->SubDomains[i*4+2]=1;
-			bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;
-		}
-	}
-	else{
-		bamgmesh->SubDomains=NULL;
-	}
-
-	//SubDomainsFromGeom
-	if(verbose>3) printf("      writing SubDomainsFromGeom\n");
-	bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains;
-	xfree((void**)&bamgmesh->SubDomainsFromGeom);
-	if (Gh.NbSubDomains){
-		bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));
-		for (i=0;i<Gh.NbSubDomains;i++){
-			bamgmesh->SubDomainsFromGeom[i*4+0]=2;
-			bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
-			bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;
-			bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;
-		}
-	}
-	else{
-		bamgmesh->SubDomainsFromGeom=NULL;
-	}
-
-	//VerticesOnGeomVertex
-	if(verbose>3) printf("      writing VerticesOnGeometricVertex\n");
-	bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex;
-	xfree((void**)&bamgmesh->VerticesOnGeometricVertex);
-	if (NbVerticesOnGeomVertex){
-		bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double));
-		for (i=0;i<NbVerticesOnGeomVertex;i++){
-			VertexOnGeom &v=VerticesOnGeomVertex[i];
-			if (!v.OnGeomVertex()){
-				throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));
-			}
-			bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing
-			bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing
-		}
-	}
-	else{
-		bamgmesh->VerticesOnGeometricVertex=NULL;
-	}
-
-	//VertexOnGeometricEdge
-	if(verbose>3) printf("      writing VerticesOnGeometricEdge\n");
-	bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge;
-	xfree((void**)&bamgmesh->VerticesOnGeometricEdge);
-	if (NbVerticesOnGeomEdge){
-		bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double));
-		for (i=0;i<NbVerticesOnGeomEdge;i++){
-			const VertexOnGeom &v=VerticesOnGeomEdge[i];
-			if (!v.OnGeomEdge()){
-				throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not"));
-			}
-			bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing
-			bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
-			bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;
-		}
-	}
-	else{
-		bamgmesh->VerticesOnGeometricEdge=NULL;
-	}
-
-	//EdgesOnGeometricEdge
-	if(verbose>3) printf("      writing EdgesOnGeometricEdge\n");
-	k=0;
-	for (i=0;i<nbe;i++){
-		if (edges[i].on) k=k+1;
-	}
-	bamgmesh->NumEdgesOnGeometricEdge=k;
-	xfree((void**)&bamgmesh->EdgesOnGeometricEdge);
-	if (k){
-		bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double));
-		int count=0;
-		for (i=0;i<nbe;i++){
-			if (edges[i].on){
-				bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing
-				bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing
-				count=count+1;
-			}
-		}
-	}
-	else{
-		bamgmesh->EdgesOnGeometricEdge=NULL;
-	}
-}
-
-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;
-	}
-}
-
-} // end of namespace bamg 
Index: /issm/trunk/src/c/Bamgx/Triangles.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Triangles.cpp	(revision 2806)
+++ /issm/trunk/src/c/Bamgx/Triangles.cpp	(revision 2806)
@@ -0,0 +1,463 @@
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <ctime>
+#include "Mesh2.h"
+#include "QuadTree.h"
+#include "SetOfE4.h"
+
+namespace bamg {
+
+	static const  Direction NoDirOfSearch=Direction();
+
+	/*Constructors/Destructors*/
+	/*FUNCTION Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
+	Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
+
+		PreInit(0,"none");
+		OnDisk = 1;
+
+		ReadFromMatlabMesh(bamgmesh,bamgopts);
+		SetIntCoor();
+		FillHoleInMesh();
+	}
+	/*}}}1*/
+
+	/*IO*/
+	/*FUNCTION Triangles::WriteMesh {{{1*/
+	void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
+
+		int i,j;
+		int verbose;
+
+		verbose=bamgopts->verbose;
+		Int4 *reft = new Int4[nbt];
+		Int4 nbInT = ConsRefTriangle(reft);
+
+		//Vertices
+		if(verbose>3) printf("      writing Vertices\n");
+		bamgmesh->NumVertices=nbv;
+		xfree((void**)&bamgmesh->Vertices);
+		if (nbv){
+			bamgmesh->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
+			for (i=0;i<nbv;i++){
+				bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
+				bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
+				bamgmesh->Vertices[i*3+2]=vertices[i].ref();
+			}
+		}
+		else{
+			bamgmesh->Vertices=NULL;
+		}
+
+		//Edges
+		if(verbose>3) printf("      writing Edges\n");
+		bamgmesh->NumEdges=nbe;
+		xfree((void**)&bamgmesh->Edges);
+		if (nbe){
+			bamgmesh->Edges=(double*)xmalloc(3*nbe*sizeof(double));
+			for (i=0;i<nbe;i++){
+				bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing
+				bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
+				bamgmesh->Edges[i*3+2]=edges[i].ref;
+			}
+		}
+		else{
+			bamgmesh->Edges=NULL;
+		}
+
+		//CrackedEdges
+		if(verbose>3) printf("      writing CrackedEdges\n");
+		bamgmesh->NumCrackedEdges=NbCrackedEdges;
+		xfree((void**)&bamgmesh->CrackedEdges);
+		if (NbCrackedEdges){
+			bamgmesh->CrackedEdges=(double*)xmalloc(2*NbCrackedEdges*sizeof(double));
+			for (i=0;i<NbCrackedEdges;i++){
+				bamgmesh->CrackedEdges[i*2+0]=Number(CrackedEdges[i].a.edge)+1; //back to M indexing
+				bamgmesh->CrackedEdges[i*2+1]=Number(CrackedEdges[i].b.edge)+1; //back to M indexing
+			}
+		}
+		else{
+			bamgmesh->CrackedEdges=NULL;
+		}
+
+		//Triangles
+		if(verbose>3) printf("      writing Triangles\n");
+		Int4 k=nbInT-NbOfQuad*2;
+		Int4 num=0;
+		bamgmesh->NumTriangles=k;
+		xfree((void**)&bamgmesh->Triangles);
+		if (k){
+			bamgmesh->Triangles=(double*)xmalloc(4*k*sizeof(double));
+			for (i=0;i<nbt;i++){
+				Triangle &t=triangles[i];
+				if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
+					bamgmesh->Triangles[num*4+0]=Number(t[0])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing
+					bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
+					num=num+1;
+				}
+			}
+		}
+		else{
+			bamgmesh->Triangles=NULL;
+		}
+
+		//Quadrilaterals
+		if(verbose>3) printf("      writing Quadrilaterals\n");
+		bamgmesh->NumQuadrilaterals=NbOfQuad;
+		xfree((void**)&bamgmesh->Quadrilaterals);
+		if (NbOfQuad){
+			bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));
+			for (i=0;i<nbt;i++){
+				Triangle &t =triangles[i];
+				Triangle* ta;
+				Vertex *v0,*v1,*v2,*v3;
+				if (reft[i]<0) continue;
+				if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { 
+					bamgmesh->Quadrilaterals[i*5+0]=Number(v0)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+1]=Number(v1)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing
+					bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;
+				}
+			}
+		}
+		else{
+			bamgmesh->Quadrilaterals=NULL;
+		}
+
+		//SubDomains
+		if(verbose>3) printf("      writing SubDomains\n");
+		bamgmesh->NumSubDomains=NbSubDomains;
+		xfree((void**)&bamgmesh->SubDomains);
+		if (NbSubDomains){
+			bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
+			for (i=0;i<NbSubDomains;i++){
+				bamgmesh->SubDomains[i*4+0]=3;
+				bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];
+				bamgmesh->SubDomains[i*4+2]=1;
+				bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;
+			}
+		}
+		else{
+			bamgmesh->SubDomains=NULL;
+		}
+
+		//SubDomainsFromGeom
+		if(verbose>3) printf("      writing SubDomainsFromGeom\n");
+		bamgmesh->NumSubDomainsFromGeom=Gh.NbSubDomains;
+		xfree((void**)&bamgmesh->SubDomainsFromGeom);
+		if (Gh.NbSubDomains){
+			bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));
+			for (i=0;i<Gh.NbSubDomains;i++){
+				bamgmesh->SubDomainsFromGeom[i*4+0]=2;
+				bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
+				bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;
+				bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;
+			}
+		}
+		else{
+			bamgmesh->SubDomainsFromGeom=NULL;
+		}
+
+		//VerticesOnGeomVertex
+		if(verbose>3) printf("      writing VerticesOnGeometricVertex\n");
+		bamgmesh->NumVerticesOnGeometricVertex=NbVerticesOnGeomVertex;
+		xfree((void**)&bamgmesh->VerticesOnGeometricVertex);
+		if (NbVerticesOnGeomVertex){
+			bamgmesh->VerticesOnGeometricVertex=(double*)xmalloc(2*NbVerticesOnGeomVertex*sizeof(double));
+			for (i=0;i<NbVerticesOnGeomVertex;i++){
+				VertexOnGeom &v=VerticesOnGeomVertex[i];
+				if (!v.OnGeomVertex()){
+					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricVertex is actually not"));
+				}
+				bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((Vertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number(( GeometricalVertex*)v)+1; //back to Matlab indexing
+			}
+		}
+		else{
+			bamgmesh->VerticesOnGeometricVertex=NULL;
+		}
+
+		//VertexOnGeometricEdge
+		if(verbose>3) printf("      writing VerticesOnGeometricEdge\n");
+		bamgmesh->NumVerticesOnGeometricEdge=NbVerticesOnGeomEdge;
+		xfree((void**)&bamgmesh->VerticesOnGeometricEdge);
+		if (NbVerticesOnGeomEdge){
+			bamgmesh->VerticesOnGeometricEdge=(double*)xmalloc(3*NbVerticesOnGeomEdge*sizeof(double));
+			for (i=0;i<NbVerticesOnGeomEdge;i++){
+				const VertexOnGeom &v=VerticesOnGeomEdge[i];
+				if (!v.OnGeomEdge()){
+					throw ErrorException(__FUNCT__,exprintf("A vertices supposed to be OnGeometricEdge is actually not"));
+				}
+				bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
+				bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;
+			}
+		}
+		else{
+			bamgmesh->VerticesOnGeometricEdge=NULL;
+		}
+
+		//EdgesOnGeometricEdge
+		if(verbose>3) printf("      writing EdgesOnGeometricEdge\n");
+		k=0;
+		for (i=0;i<nbe;i++){
+			if (edges[i].on) k=k+1;
+		}
+		bamgmesh->NumEdgesOnGeometricEdge=k;
+		xfree((void**)&bamgmesh->EdgesOnGeometricEdge);
+		if (k){
+			bamgmesh->EdgesOnGeometricEdge=(double*)xmalloc(2*(int)k*sizeof(double));
+			int count=0;
+			for (i=0;i<nbe;i++){
+				if (edges[i].on){
+					bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing
+					bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].on)+1; //back to Matlab indexing
+					count=count+1;
+				}
+			}
+		}
+		else{
+			bamgmesh->EdgesOnGeometricEdge=NULL;
+		}
+	}
+	/*}}}1*/
+	/*FUNCTION Triangles::ReadFromMatlabMesh{{{1*/
+	void Triangles::ReadFromMatlabMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){
+
+		int verbose;
+		Real8 Hmin = HUGE_VAL;// the infinie value 
+		Int4 i1,i2,i3,iref;
+		Int4 i,j;
+		Int4 hvertices =0;
+		Int4 ifgeom=0;
+		Metric M1(1);
+		int Version=1,dim=2;
+
+		verbose=bamgopts->verbose;
+
+		nbv=bamgmesh->NumVertices;
+		nbvx=nbv;
+		nbt=bamgmesh->NumTriangles;
+		nbiv=nbvx;
+
+		//Vertices
+		if(bamgmesh->Vertices){
+			if(verbose>3) printf("      processing Vertices\n");
+
+			vertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
+			ordre=(Vertex**)xmalloc(nbv*sizeof(Vertex*));
+
+			for (i=0;i<nbv;i++){
+				vertices[i].r.x=bamgmesh->Vertices[i*3+0];
+				vertices[i].r.y=bamgmesh->Vertices[i*3+1];
+				vertices[i].ReferenceNumber=1;
+				vertices[i].DirOfSearch =NoDirOfSearch;
+				vertices[i].m=M1;
+				vertices[i].color=(Int4)bamgmesh->Vertices[i*3+2];
+			}
+			nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
+		}
+		else{
+			if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Vertices found in the initial mesh"));
+		}
+
+		//Triangles
+		if(bamgmesh->Triangles){
+			if(verbose>3) printf("      processing Triangles\n");
+			triangles =new Triangle[nbtx]; //we cannot allocate only nbt triangles since 
+			//other triangles will be added for each edge
+			for (i=0;i<nbt;i++){
+				Triangle & t = triangles[i];
+				i1=(Int4)bamgmesh->Triangles[i*4+0]-1; //for C indexing
+				i2=(Int4)bamgmesh->Triangles[i*4+1]-1; //for C indexing
+				i3=(Int4)bamgmesh->Triangles[i*4+2]-1; //for C indexing
+				t=Triangle(this,i1,i2,i3);
+				t.color=(Int4)bamgmesh->Triangles[i*4+3];
+			}
+		}
+		else{
+			if(verbose>3) throw ErrorException(__FUNCT__,exprintf("no Triangles found in the initial mesh"));
+		}
+
+		//Quadrilaterals
+		if(bamgmesh->Quadrilaterals){
+			if(verbose>3) printf("      processing Quadrilaterals\n");
+			Int4 i1,i2,i3,i4,iref;
+			triangles =new Triangle[nbt];
+			for (i=0;i<bamgmesh->NumQuadrilaterals;i++){
+				Triangle & t1 = triangles[2*i];
+				Triangle & t2 = triangles[2*i+1];
+				i1=(Int4)bamgmesh->Quadrilaterals[i*4+0]-1; //for C indexing
+				i2=(Int4)bamgmesh->Quadrilaterals[i*4+1]-1; //for C indexing
+				i3=(Int4)bamgmesh->Quadrilaterals[i*4+2]-1; //for C indexing
+				i4=(Int4)bamgmesh->Quadrilaterals[i*4+3]-1; //for C indexing
+				t1=Triangle(this,i1,i2,i3);
+				t2=Triangle(this,i3,i4,i1);
+				t1.color=1; //reference = 1 for all triangles since it has not been specified
+				t2.color=1; //reference = 1 for all triangles since it has not been specified
+				t1.SetHidden(OppositeEdge[1]); // two times  because the adj was not created 
+				t2.SetHidden(OppositeEdge[1]); 
+			}
+		}
+		else{
+			if(verbose>3) printf("      no Quadrilaterals found\n");
+		}
+
+		//hVertices
+		if(bamgmesh->hVertices){
+			if(verbose>3) printf("      processing hVertices\n");
+			hvertices=1;
+			for (i=0;i< nbv;i++){
+				vertices[i].m=Metric((Real4)bamgmesh->hVertices[i]);
+			}
+		}
+		else{
+			if(verbose>3) printf("      no hVertices found\n");
+		}
+
+		//VerVerticesOnGeometricEdge
+		if(bamgmesh->VerticesOnGeometricEdge){
+			if(verbose>3) printf("      processing VerticesOnGeometricEdge\n");
+			NbVerticesOnGeomEdge=bamgmesh->NumVerticesOnGeometricEdge;
+			VerticesOnGeomEdge= new  VertexOnGeom[NbVerticesOnGeomEdge] ;
+			for (i=0;i<NbVerticesOnGeomEdge;i++){
+				Int4  i1,i2;
+				Real8 s;
+				i1=(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing
+				i2=(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing
+				s =(Int4)bamgmesh->VerticesOnGeometricEdge[i*3+2];
+				VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
+			}
+		}
+		else{
+			if(verbose>3) printf("      no VertexOnGeometricEdge found\n");
+		}
+
+		//Edges
+		if (bamgmesh->Edges){
+			int i1,i2;
+			R2 zero2(0,0);
+			Real4 *len =0;
+
+			if(verbose>3) printf("      processing Edges\n");
+			nbe=bamgmesh->NumEdges;
+			edges = new Edge[nbe];
+
+			if (!hvertices) {
+				len = new Real4[nbv];
+				for(i=0;i<nbv;i++)
+				 len[i]=0;
+			}
+
+			for (i=0;i<nbe;i++){
+				i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
+				i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
+				edges[i].v[0]= vertices +i1;
+				edges[i].v[1]= vertices +i2;
+				edges[i].adj[0]=0;
+				edges[i].adj[1]=0;
+				R2 x12 = vertices[i2].r-vertices[i1].r;
+				Real8 l12=sqrt( (x12,x12));        
+
+				if (!hvertices) {
+					vertices[i1].color++;
+					vertices[i2].color++;
+					len[i1]+= l12;
+					len[i2] += l12;
+				}
+				Hmin = Min(Hmin,l12);
+			}
+
+			// definition  the default of the given mesh size 
+			if (!hvertices){
+				for (i=0;i<nbv;i++) 
+				 if (vertices[i].color > 0) 
+				  vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
+				 else 
+				  vertices[i].m=  Metric(Hmin);
+				delete [] len;
+			}
+
+			// construction of edges[].adj 
+			for (i=0;i<nbv;i++){ 
+				vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
+			}
+			for (i=0;i<nbe;i++){
+				for (j=0;j<2;j++) { 
+					Vertex *v=edges[i].v[j];
+					Int4 i0=v->color,j0;
+					if(i0==-1){
+						v->color=i*2+j;
+					}
+					else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
+						j0 =  i0%2;
+						i0 =  i0/2;
+						assert( v ==  edges[i0 ].v[j0]);
+						edges[i ].adj[ j ] =edges +i0;
+						edges[i0].adj[ j0] =edges +i ;
+						assert(edges[i0].v[j0] == v);
+						v->color = -3;
+					}
+				}
+			}
+		}
+		else{
+			if(verbose>3) printf("      no Edges found\n");
+		}
+
+		//EdgeOnGeometricEdge
+		if(bamgmesh->EdgesOnGeometricEdge){
+			if(verbose>3) printf("      processing EdgesOnGeometricEdge\n");
+			int i1,i2,i,j;
+			i2=bamgmesh->NumEdgesOnGeometricEdge;
+			for (i1=0;i1<i2;i1++) {
+				i=(int)bamgmesh->EdgesOnGeometricEdge[i*2+0];
+				j=(int)bamgmesh->EdgesOnGeometricEdge[i*2+1];
+				if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe)) {
+					throw ErrorException(__FUNCT__,exprintf("EdgesOnGeometricEdge error: We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe)"));
+				}
+				edges[i-1].on = Gh.edges + j-1;
+			}
+		}
+		else{
+			if(verbose>3) printf("      no EdgesOnGeometricEdge found\n");
+		}
+
+		//SubDomain
+		if(bamgmesh->SubDomains){
+			Int4 i3,head,sens;
+			if(verbose>3) printf("      processing SubDomains\n");
+			NbSubDomains=bamgmesh->NumSubDomains;
+			subdomains = new SubDomain [ NbSubDomains ];
+			for (i=0;i<NbSubDomains;i++) {
+				i3  =(int)bamgmesh->SubDomains[i*3+0];
+				head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
+				sens=(int)bamgmesh->SubDomains[i*3+2];
+				if (i3!=23) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 3"));
+				if (head<0 || head>=nbt) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: head should in [1 %i] (triangle number)",nbt));
+				subdomains[i].head = triangles+head;
+			}
+		}
+		else{
+			if(verbose>3) printf("      no SubDomains found\n");
+		}
+
+		/*Recreate geometry if needed*/
+		if(1!=0) {
+			printf("Warning: mesh present but no geometry found. Reconstructing...\n");
+			/*Recreate geometry: */
+			ConsGeometry(bamgopts->MaximalAngleOfCorner*Pi/180);
+			Gh.AfterRead();
+		}
+	}
+	/*}}}1*/
+
+} // end of namespace bamg 
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 2805)
+++ /issm/trunk/src/c/Makefile.am	(revision 2806)
@@ -324,7 +324,7 @@
 					./Bamgx/MeshGeom.cpp \
 					./Bamgx/MeshQuad.cpp \
-					./Bamgx/MeshRead.cpp \
 					./Bamgx/meshtype.h \
-					./Bamgx/MeshWrite.cpp	\
+					./Bamgx/Triangles.cpp	\
+					./Bamgx/Geometry.cpp	\
 					./Bamgx/Metric.cpp \
 					./Bamgx/Metric.h \
@@ -664,7 +664,7 @@
 					./Bamgx/MeshGeom.cpp \
 					./Bamgx/MeshQuad.cpp \
-					./Bamgx/MeshRead.cpp \
 					./Bamgx/meshtype.h \
-					./Bamgx/MeshWrite.cpp	\
+					./Bamgx/Triangles.cpp	\
+					./Bamgx/Geometry.cpp	\
 					./Bamgx/Metric.cpp \
 					./Bamgx/Metric.h \