Index: /issm/trunk/src/c/Bamgx/Mesh2.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2794)
+++ /issm/trunk/src/c/Bamgx/Mesh2.h	(revision 2795)
@@ -889,21 +889,7 @@
   Vertex * NearestVertex(Icoor1 i,Icoor1 j) ;
   Triangle * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
-  void Write(const char * filename,const TypeFileMesh type = AutoMesh);
-  void Write_am_fmt(std::ostream &) const ;
-  void Write_am(std::ostream &) const ;
-  void Write_ftq(std::ostream &) const ;
-  void Write_nopo(std::ostream &) const ;
-  void Write_msh(std::ostream &) const ;
-  void Write_amdba(std::ostream &) const ;
-
-  void Read(MeshIstream &,int version,Real8 cutoffradian);
+
   void ReadFromMatlabMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
   void WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
-  void Read_am_fmt(MeshIstream &);
-  void Read_amdba(MeshIstream &);
-  void Read_am(MeshIstream &);
-  void Read_nopo(MeshIstream &);
-  void Read_ftq(MeshIstream &);
-  void Read_msh(MeshIstream &);
 
   void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
Index: /issm/trunk/src/c/Bamgx/MeshRead.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshRead.cpp	(revision 2794)
+++ /issm/trunk/src/c/Bamgx/MeshRead.cpp	(revision 2795)
@@ -268,876 +268,4 @@
 }
 
-void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian) {
-
-
-	long int verbosity=0;
-
-  Real8 hmin = HUGE_VAL;// the infinie value 
-  //  Real8 MaximalAngleOfCorner = 10*Pi/180;// 
-  Int4 i;
-  Int4 dim=0;
-  Int4 hvertices =0;
-  Int4 ifgeom=0;
-  Metric M1(1);
-  if (verbosity>1)
-    cout << "  -- ReadMesh " << f_in.CurrentFile  << " Version = " << Version << endl;
-  int field=0;
-  int showfield=0;
-  while (f_in.cm()) 
-    {  
-      field=0;
-      char  fieldname[256];
-      if(f_in.eof()) break;
-      f_in.cm() >> fieldname ;;
-      if(f_in.eof()) break;
-      f_in.err() ;
-      if(verbosity>9)
-	cout <<  "    " << fieldname << endl;
-      if (!strcmp(fieldname,"MeshVersionFormatted") )
-         f_in >>   Version ;
-      else if(!strcmp(fieldname,"End"))
-         break;
-      else if (!strcmp(fieldname,"Dimension"))
-         {
-           f_in >>   dim ;
-           assert(dim ==2);
-         }
-      else if  (!strcmp(fieldname,"Geometry"))
-	{ 
-	  char * fgeom ;
-	  f_in >> fgeom;
-	  //	  if (cutoffradian>=0) => bug if change edit the file geometry 
-	  //  Gh.MaximalAngleOfCorner = cutoffradian;
-	  if (strlen(fgeom))
-	      Gh.ReadGeometry(fgeom);
-	  else 
-	    { 
-	      // include geometry 
-	      f_in.cm();
-	      Gh.ReadGeometry(f_in,fgeom);
-	    }
-
-	  Gh.AfterRead();
-	  ifgeom=1;
-	  delete [] fgeom;
-	}
-      else if (!strcmp(fieldname,"Identifier"))
-	{
-	  if (identity) delete [] identity;
-	  f_in >> identity;
-	}
-      else if (!strcmp(fieldname,"hVertices"))
-         { hvertices =1;
-	   Real4 h;
-           for (i=0;i< nbv;i++) {
-            f_in >>  h ;
-	    vertices[i].m = Metric(h);}
-         }
-      else if (!strcmp(fieldname,"Vertices"))
-         { 
-           assert(dim ==2);
-           f_in >>   nbv ;
-	   if(verbosity>3)
-	     cout << "   Nb of Vertices = " << nbv << endl;
-	   nbvx=nbv;
-	   vertices=new Vertex[nbvx];
-	   assert(vertices);
-	   ordre=new (Vertex* [nbvx]);
-	   assert(ordre);
-	   
-	   nbiv = nbv;
-           for (i=0;i<nbv;i++) {
-             f_in >> 
-	       	 	 vertices[i].r.x >>  
-	           vertices[i].r.y >> 
-             vertices[i].ReferenceNumber  ;
-             vertices[i].DirOfSearch =NoDirOfSearch;
-	           vertices[i].m=M1;
-	     vertices[i].color =0;}
-	   nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-	   triangles =new Triangle[nbtx];
-	   assert(triangles);
-	   nbt =0;
-         }
-      else if (!strcmp(fieldname,"Triangles"))
-	{ 
-	  if (dim !=2)
-	    cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
-	  if(! vertices || !triangles  || !nbv )
-	    cerr << "ReadMesh:Triangles before Vertices" <<endl,
-	      f_in.ShowIoErr(0);
-	  int NbOfTria;
-	  f_in >>  NbOfTria ;
-	  if(verbosity>3)
-	    cout << "   NbOfTria = " << NbOfTria << endl;
-	  if (nbt+NbOfTria >= nbtx)
-	    cerr << "ReadMesh: We must have 2*NbOfQuad + NbOfTria  = "
-		 << nbt + NbOfTria<<"  < 2*nbv-2 ="  << nbtx << endl,
-	      f_in.ShowIoErr(0); 
-	  //	  begintria = nbt;
-	  for (i=0;i<NbOfTria;i++)  
-	    {
-	      Int4 i1,i2,i3,iref;
-	      Triangle & t = triangles[nbt++];
-	      f_in >>  i1 >>  i2 >>   i3 >>   iref ;
-	      t = Triangle(this,i1-1,i2-1,i3-1);
-	      t.color=iref;
-	    }
-	  // endtria=nbt;
-	}
-      else if (!strcmp(fieldname,"Quadrilaterals"))
-        { 
-
-	  if (dim !=2)
-	    cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
-	  if(! vertices || !triangles  || !nbv )
-	    cerr << "ReadMesh:Quadrilaterals before Vertices" <<endl,
-	      f_in.ShowIoErr(0);
-	  f_in >>   NbOfQuad ;
-	  if(verbosity>3)
-	    cout << "   NbOfQuad= " << NbOfQuad << endl;
-	  if (nbt+2*NbOfQuad >= nbtx)
-	    cerr << "ReadMesh: We must have 2*NbOfQuad + NbOfTria  = "
-		 << nbt + 2*NbOfQuad <<"  < 2*nbv-2 ="  << nbtx << endl,
-	      f_in.ShowIoErr(0);
-	  //  beginquad=nbt;
-	  for (i=0;i<NbOfQuad;i++)  
-	    { 
-	      Int4 i1,i2,i3,i4,iref;
-	      Triangle & t1 = triangles[nbt++];
-	      Triangle & t2 = triangles[nbt++];
-	      f_in >>  i1 >>  i2 >>   i3 >> i4 >>    iref ;
-	      t1 = Triangle(this,i1-1,i2-1,i3-1);
-	      t1.color=iref;
-	      t2 = Triangle(this,i3-1,i4-1,i1-1);
-	      t2.color=iref;
-	      t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created 
-	      t2.SetHidden(OppositeEdge[1]); 
-	    }
-	  // endquad=nbt;
-	}
-      else if (!strcmp(fieldname,"VertexOnGeometricVertex"))
-	{ 
-           f_in  >> NbVerticesOnGeomVertex ;
- 	   if(verbosity>5)
-	     cout << "   NbVerticesOnGeomVertex = "   << NbVerticesOnGeomVertex << endl
-		  << " Gh.vertices " << Gh.vertices << endl;
-	   if( NbVerticesOnGeomVertex) 
-	     {
-	       VerticesOnGeomVertex= new  VertexOnGeom[NbVerticesOnGeomVertex] ;
-	       if(verbosity>7)
-		 cout << "   VerticesOnGeomVertex = " << VerticesOnGeomVertex << endl
-		      << "   Gh.vertices " << Gh.vertices << endl;
-	       assert(VerticesOnGeomVertex);
-	       for (Int4 i0=0;i0<NbVerticesOnGeomVertex;i0++)  
-		 { 
-		   Int4  i1,i2;
-		   //VertexOnGeom & v =VerticesOnGeomVertex[i0];
-		   f_in >>  i1 >> i2 ;
-		   VerticesOnGeomVertex[i0]=VertexOnGeom(vertices[i1-1],Gh.vertices[i2-1]);
-		 }
-	     }
-	 }
-      else if (!strcmp(fieldname,"VertexOnGeometricEdge"))
-         { 
-           f_in >>  NbVerticesOnGeomEdge ;
- 	   if(verbosity>3)
-	     cout << "   NbVerticesOnGeomEdge = " << NbVerticesOnGeomEdge << endl;
-	   if(NbVerticesOnGeomEdge) 
-	     {
-	       VerticesOnGeomEdge= new  VertexOnGeom[NbVerticesOnGeomEdge] ;
-	       assert(VerticesOnGeomEdge);
-	       for (Int4 i0=0;i0<NbVerticesOnGeomEdge;i0++)  
-		 { 
-		   Int4  i1,i2;
-		   Real8 s;
-		   //VertexOnGeom & v =VerticesOnGeomVertex[i0];
-		   f_in >>  i1 >> i2  >> s;
-		   VerticesOnGeomEdge[i0]=VertexOnGeom(vertices[i1-1],Gh.edges[i2-1],s);
-		 }
-	     }
-         }
-      else if (!strcmp(fieldname,"Edges"))
-	{ 
-          Int4 i,j, i1,i2;
-          f_in >>  nbe ;
-          edges = new Edge[nbe];
-	  if(verbosity>5)
-	    cout << "     Record Edges: Nb of Edge " << nbe << " edges " <<  edges << endl;
-          assert(edges);
-	  Real4 *len =0;
-          if (!hvertices) 
-	    {
-	      len = new Real4[nbv];
-	      for(i=0;i<nbv;i++)
-		len[i]=0;
-	    }
-          for (i=0;i<nbe;i++) 
-	    {
-	      f_in  >> i1  >> i2  >> edges[i].ref ;
-                 
-	      assert(i1 >0 && i2 >0);
-	      assert(i1 <= nbv && i2 <= nbv);
-	      i1--;
-	      i2--;
-	      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;
-	    }
-	  if(verbosity>5)
-	    cout << "     hmin " << hmin << endl;
-	  // 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);
-		    //	    if(verbosity>8)
-		    //  cout << " edges adj " << i0 << " "<< j0 << " <-->  "  << i << " " << j << endl;
-		      v->color = -3;}
-		}
-
-	}
-	
-/*   ne peut pas marche car il n'y a pas de flag dans les aretes du maillages
-       else if (!strcmp(fieldname,"RequiredEdges"))
-        { 
-          int i,j,n;
-          f_in  >> n ;
-
-          for (i=0;i<n;i++) {     
-            f_in  >>  j ;
-            assert( j <= nbe );
-            assert( j > 0 );
-            j--;
-            edges[j].SetRequired();  }
-      }
-*/	
-      else if (!strcmp(fieldname,"EdgeOnGeometricEdge"))
-	{ 
-	  assert(edges);
-	  int i1,i2,i,j;
-	  f_in >> i2 ;
-	  if(verbosity>3)
-	    cout << "     Record EdgeOnGeometricEdge: Nb " << i2 <<endl;
-	  for (i1=0;i1<i2;i1++)
-	    {
-	      f_in  >> i >> j ;
-	      if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe))
-		{
-		  cerr << "      Record EdgeOnGeometricEdge i=" << i << " j = " << j;
-		  cerr << " nbe = " << nbe << " Gh.nbe = " <<  Gh.nbe << endl;
-		  cerr << " We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe) ";
-		  cerr << " Fatal error in file " << name << " line " << f_in.LineNumber << endl;
-		  MeshError(1);
-		}
-		
-		
-	      edges[i-1].on = Gh.edges + j-1;
-	    }
-	  //	  cout << "end EdgeOnGeometricEdge" << endl;
-	}
-      else if (!strcmp(fieldname,"SubDomain") || !strcmp(fieldname,"SubDomainFromMesh") )
-	{ 
-	  
-	  f_in >> NbSubDomains ;
-	  subdomains = new SubDomain [ NbSubDomains ];
-	    for (i=0;i< NbSubDomains;i++)
-	      { Int4 i3,head,sens;
-	      f_in  >>  i3 >>	head >> sens  >> subdomains[i].ref ;
-		assert (i3==3);
-		head --;
-		assert(head < nbt && head >=0);
-		subdomains[i].head = triangles+head;		
-	      }
-	}
-      else
-	{ // unkown field
-	  field = ++showfield;
-	  if(showfield==1) // just to show one time 
-	    if (verbosity>5)
-	      cout << "     Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;
-	}
-      showfield=field; // just to show one time 
-    }
-    
-    if (ifgeom==0)
-      {
-	if (verbosity)
-	  cout  << " ## Warning: Mesh without geometry we construct a geometry (theta =" 
-		<< cutoffradian*180/Pi << " degres )" << endl;
-	ConsGeometry(cutoffradian);	
-	Gh.AfterRead();
-      }	  
-}
-
-
-void Triangles::Read_am_fmt(MeshIstream & f_in)
-{
-	long int verbosity=0;
-
-  Metric M1(1);
-
-  if (verbosity>1)
-    cout << "  -- ReadMesh .am_fmt file " <<  f_in.CurrentFile  << endl;
-   	
-     Int4 i;
-     f_in.cm() >> nbv >> nbt ;
-     if (verbosity>3)
-       cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
-     f_in.eol() ;// 
-     nbvx = nbv;
-     nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-     triangles =new Triangle[nbtx];
-     assert(triangles);
-     vertices=new Vertex[nbvx];
-     ordre=new (Vertex* [nbvx]);
-     
-     for (     i=0;i<nbt;i++)
-       {
-	 Int4 i1,i2,i3;
-	 f_in >>  i1 >>  i2 >>   i3 ;
-	 triangles[i]  = Triangle(this,i1-1,i2-1,i3-1);	 
-       }
-      f_in.eol() ;// 
-      
-      for ( i=0;i<nbv;i++)
-	f_in >> vertices[i].r.x >>   vertices[i].r.y,
-	  vertices[i].m = M1,vertices[i].DirOfSearch =NoDirOfSearch;
-
-      f_in.eol() ;// 
-      
-      for ( i=0;i<nbt;i++)  
-	f_in >> triangles[i].color;
-      f_in.eol() ;// 
-      
-      for ( i=0;i<nbv;i++)  
-	     f_in >> vertices[i].ReferenceNumber;
-      
-      
-}
-
-////////////////////////
-
-void  Triangles::Read_am(MeshIstream &ff)
-{
-	long int verbosity=0;
-
-  if (verbosity>1)
-    cout << "  -- ReadMesh .am_fmt file " <<  ff.CurrentFile  << endl;
-    Metric M1(1);	
-  
-  IFortranUnFormattedFile f_in(ff);
-
-  Int4  l=f_in.Record();
-  assert(l==2*sizeof(Int4));
-  f_in >> nbv >> nbt ;
-  l=f_in.Record();
-  assert((size_t) l==nbt*sizeof(long)*4 + nbv*(2*sizeof(float)+sizeof(long)));
-  if (verbosity>3)
-    cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
-  
-  nbvx = nbv;
-  nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-  triangles =new Triangle[nbtx];
-  assert(triangles);
-  vertices=new Vertex[nbvx];
-  ordre=new (Vertex* [nbvx]);
-  
-
-  Int4 i;
-  for (     i=0;i<nbt;i++) {
-    long i1,i2,i3;
-    f_in >>  i1 >>  i2 >>   i3 ;
-    triangles[i]  = Triangle(this,i1-1,i2-1,i3-1); }
-  
-  for ( i=0;i<nbv;i++) {
-    float x,y;
-    f_in >> x >> y;
-    vertices[i].r.x =x;
-    vertices[i].r.y=y;
-    vertices[i].m=M1;}
-  
-  for ( i=0;i<nbt;i++) {
-    long i;
-    f_in >> i;
-    triangles[i].color=i;}
-  
-  for ( i=0;i<nbv;i++) {
-    long i;
-    f_in >> i;
-    vertices[i].ReferenceNumber=i;}
-}
-
-//////////////////////////////////
-
-void  Triangles::Read_nopo(MeshIstream & ff)
-{
-	long int verbosity=0;
-
-
- if (verbosity>1)
-    cout << "  -- ReadMesh .nopo file " <<  ff.CurrentFile  << endl;
- IFortranUnFormattedFile f_in(ff);
- 
- 
- Int4  l,i,j;
- l=f_in.Record(); 
- l=f_in.Record();
- f_in >> i;
- assert(i==32);
- Int4 niveau,netat,ntacm;
-
- char titre[80+1],  date[2*4+1], nomcre[6*4+1], typesd[5];
- f_in.read4(titre,20);
- f_in.read4(date,2);
- f_in.read4(nomcre,6);
- f_in.read4(typesd,1);
-
-
-   f_in >> niveau>>netat>>ntacm;
- if (strcmp("NOPO",typesd))
-   {
-     cout << " where in record  " << f_in.where() << " " << strcmp("NOPO",typesd) << endl;
-     cerr << " not a  nopo file but `" << typesd <<"`"<< " len = " << strlen(typesd) << endl;
-     cerr << (int) typesd[0] << (int) typesd[1] << (int) typesd[2] << (int) typesd[3] << (int) typesd[4] << endl;
-     cout << " nomcre :" << nomcre << endl;
-     cout << " date   :" << date << endl;
-     cout << " titre  :" << titre<< endl;
-     MeshError(112);
-   }
- if(verbosity>2)
-   cout << "    nb de tableau associe : " << ntacm << " niveau =" << niveau << endl;
- 
- for (i=0;i<ntacm;i++)
-   f_in.Record();
-
- f_in.Record();
- f_in >> l;
- assert(l == 27);
- Int4 nop2[27];
- for (i=0;i<27;i++)
-   f_in >> nop2[i];
- Int4 ndim = nop2[0];
- Int4 ncopnp = nop2[3];
- Int4 ne = nop2[4];
- Int4 ntria = nop2[7];
- Int4 nquad = nop2[8];
- Int4 np = nop2[21];
- // Int4 nef = nop2[13];
- Metric M1(1);
- if(verbosity>2) 
-   cout << "    ndim = " << ndim << " ncopnp= " << ncopnp << " ne = " << ne 
-	<< "    ntri = " << ntria << " nquad = " << nquad << " np = " << np << endl;
- nbv = np;
- nbt = 2*nquad + ntria;
- if (ne != nquad+ntria || ndim !=2 || ncopnp != 1 )
-   {
-     cerr << " not only tria & quad in nopo mesh on dim != 2 ou ncopnp != 1 " << endl;
-     MeshError(113);
-   }
- if( nop2[24]>=0)  f_in.Record();
-  NbOfQuad = nquad;
-  nbvx = nbv;
-  nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-  triangles =new Triangle[nbtx];
-  assert(triangles);
-  vertices=new Vertex[nbvx];
-  ordre=new (Vertex* [nbvx]);
-
-
- f_in >> l;
-
-  if(verbosity>9)
-    cout << " Read cnop4 nb of float  " << l << endl;
-  
-  assert(l==2*np);
-  for (i=0;i<np;i++)
-    {  float x,y;
-    f_in >>  x>> y;
-    vertices[i].r.x=x;
-    vertices[i].r.y=y;
-    vertices[i].m=M1;
-    vertices[i].DirOfSearch =NoDirOfSearch;
-
-    }
-  f_in.Record();
-  // lecture de nop5 bonjour les degats
-  f_in >> l;
-  if(verbosity>9)
-    cout << " Read nop5  nb of int4 " << l << endl;
- Int4 k=0;
- Int4 nbe4 =  3*ntria + 4*nquad;
- // cout << " nbv = " << nbv << " nbe4 " << nbe4 << endl;
- SetOfEdges4 * edge4= new SetOfEdges4(nbe4,nbv); 
- Int4 * refe = new Int4[nbe4];
- Int4 kr =0;
- for (i=0;i<ne;i++)
-   {
-     // Int4 ng[4]={0,0,0,0};
-     Int4 np[4],rv[4],re[4];
-     Int4 ncge,nmae,ndsde,npo;
-     f_in >> ncge >> nmae >> ndsde >> npo ;
-     //cout << " element " << i << " " << ncge << " "  
-     // << nmae <<" " <<  npo << endl;
-     if (ncge != 3 && ncge != 4)
-       {
-	 cerr << " read nopo type element[" << i << "] =" 
-	      << ncge << " not 3 or 4 " << endl;
-	 MeshError(115);
-       }
-     if (npo != 3 && npo != 4)
-       {
-	 cerr << " read nopo element[" << i << "] npo = "  
-	      << npo << " not 3 or 4 " << endl;
-	 MeshError(115);
-       }
-     
-     for( j=0;j<npo;j++)
-       {f_in >>np[j];np[j]--;}
-     
-     if (ncopnp !=1) 
-       {
-	 f_in >> npo;
-	 if (npo != 3 || npo != 4)
-	   {
-	     cerr << " read nopo type element[" << i << "]= "  
-		  << ncge << " not 3 or 4 " << endl;
-	     MeshError(115);
-	   }
-	 
-	 for(j=0;j<npo;j++)
-	   {f_in >>np[j];np[j]--;}
-	 
-       }
-     if (nmae>0) 
-       {
-	 Int4  ining; // no ref 
-	 
-	 f_in>>ining;
-	 if (ining==1)
-	   MeshError(116);
-	 if (ining==2)
-	   for (j=0;j<npo;j++)
-	     f_in >>re[j];
-	 for (j=0;j<npo;j++)
-	   f_in >>rv[j];
-	 
-	 
-	 // set the ref on vertex and the shift np of -1 to start at 0
-	 for (j=0;j<npo;j++)
-	   vertices[np[j]].ReferenceNumber = rv[j];
-	 
-	 if (ining==2)
-	   for (j=0;j<npo;j++)
-	     if (re[j])
-	       {
-		 kr++;
-		 Int4 i0 = np[j];
-		 Int4 i1 = np[(j+1)%npo];
-		 // cout << kr << " ref  edge " << i0 << " " << i1 << " " << re[j] << endl;
-		 refe[edge4->addtrie(i0,i1)]=re[j];
-	       }
-       }
-     
-     if (npo==3) 
-       { // triangles 
-	 triangles[k]  = Triangle(this,np[0],np[1],np[2]); 
-	 triangles[k].color = ndsde;
-	 k++;
-       }
-     else if (npo==4)
-       { // quad 
-	 Triangle & t1 = triangles[k++];
-	 Triangle & t2 = triangles[k++];
-	 t1 = Triangle(this,np[0],np[1],np[2]);
-	 t2 = Triangle(this,np[2],np[3],np[0]);
-	 t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created 
-	 t2.SetHidden(OppositeEdge[1]); 
-	 t1.color = ndsde;
-	 t2.color = ndsde;
-       }
-     else
-       {
-	 cerr << " read nopo type element =" << npo << " not 3 or 4 " << endl;
-	 MeshError(114);
-       }
-     
-     
-   }
- // cout << k << " == " << nbt << endl;
- assert(k==nbt);
-
- nbe = edge4->nb();
- if (nbe)
-   {
-     if (verbosity>7)
-     cout << " Nb of ref edges = " << nbe << endl;
-     if (edges)
-       delete [] edges;
-     edges = new Edge[nbe];
-     for (i=0;i<nbe;i++)
-       {
-	 edges[i].v[0] = vertices + edge4->i(i);
-	 edges[i].v[1] = vertices + edge4->j(i);
-	 edges[i].ref = refe[i];
-	 //	 cout << i << " " <<  edge4->i(i) << " " <<  edge4->j(i) << endl;
-       }
-      if (verbosity>7)
-	cout << " Number of reference edge in the  mesh = " << nbe << endl;
-   }
- delete [] refe;
- delete edge4;
-}
-  void  Triangles::Read_ftq(MeshIstream & f_in)
-{
-	long int verbosity=0;
-
-  //  
-  if (verbosity>1)
-    cout << "  -- ReadMesh .ftq file " <<  f_in.CurrentFile  << endl;
-  
-  Int4 i,ne,nt,nq;
-  f_in.cm() >> nbv >> ne >> nt >> nq ;
-  if (verbosity>3)
-    cout << "    nbv = " << nbv  << " nbtra = " << nt << " nbquad = " << nq << endl;
-  nbt = nt+2*nq;
-  
-
-  nbvx = nbv;
-  nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-  triangles =new Triangle[nbtx];
-  assert(triangles);
-  vertices=new Vertex[nbvx];
-  ordre=new (Vertex* [nbvx]);
-  Int4 k=0;
-  
-  for ( i=0;i<ne;i++) 
-    {
-      long ii,i1,i2,i3,i4,ref;
-      f_in >>  ii;
-	if (ii==3) 
-	  { // triangles 
-	    f_in >> i1>>  i2 >>   i3 >> ref ;
-	    triangles[k]  = Triangle(this,i1-1,i2-1,i3-1); 
-	    triangles[k++].color = ref;
-	  }
-	else if (ii==4)
-	  { // quad 
-	    f_in >> i1>>  i2 >>   i3 >> i4 >> ref ;
-	    Triangle & t1 = triangles[k++];
-	    Triangle & t2 = triangles[k++];
-	    t1 = Triangle(this,i1-1,i2-1,i3-1);
-	    t1.color=ref;
-	    t2 = Triangle(this,i3-1,i4-1,i1-1);
-	    t2.color=ref;
-	    t1.SetHidden(OppositeEdge[1]); // two time  because the adj was not created 
-	    t2.SetHidden(OppositeEdge[1]); 	  
-	  }
-	else
-	  {
-	    cout << " read ftq type element =" << ii << " not 3 or 4 " << endl;
-	    MeshError(111);
-	  }
-    }
-  assert(k==nbt);
-  Metric M1(1);
-  for ( i=0;i<nbv;i++)
-    {
-      f_in >> vertices[i].r.x >>   vertices[i].r.y >> vertices[i].ReferenceNumber;
-       vertices[i].DirOfSearch =NoDirOfSearch;
-       vertices[i].m = M1;
-
-    }
-}
-
-///////////////////////////////////////////////
-
-void  Triangles::Read_msh(MeshIstream &f_in)
-{
-	long int verbosity=0;
-
-    Metric M1(1.);
-  if (verbosity>1)
-    cout << "  -- ReadMesh .msh file " <<  f_in.CurrentFile  << endl;
-   	
-     Int4 i;
-     f_in.cm() >> nbv >> nbt ;
-     while (f_in.in.peek()==' ')
-         f_in.in.get();
-     if(isdigit(f_in.in.peek())) 
-       f_in >> nbe;
-     if (verbosity>3)
-       cout << "    nbv = " << nbv  << " nbt = " << nbt << " nbe = " << nbe << endl;
-     nbvx = nbv;
-     nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-     triangles =new Triangle[nbtx];
-     assert(triangles);
-     vertices=new Vertex[nbvx];
-     ordre=new (Vertex* [nbvx]);
-      edges = new Edge[nbe];
-     for ( i=0;i<nbv;i++)
-	{
-	 f_in >> vertices[i].r.x >>   vertices[i].r.y >> vertices[i].ReferenceNumber;
-	    vertices[i].on=0;
-	    vertices[i].m=M1;
-         //if(vertices[i].ReferenceNumber>NbRef)	NbRef=vertices[i].ReferenceNumber;  
-  	}
-     for (     i=0;i<nbt;i++)
-       {
-	 Int4 i1,i2,i3,r;
-	 f_in >>  i1 >>  i2 >>   i3 >> r;
-	 triangles[i]  = Triangle(this,i1-1,i2-1,i3-1);
-	 triangles[i].color = r;	 
-       }
-     for (i=0;i<nbe;i++)
-       {
-	 Int4 i1,i2,r;
-	 f_in >>  i1 >>  i2 >> r;
-	      edges[i].v[0]= vertices +i1-1;
-	      edges[i].v[1]= vertices +i2-1;
-	      edges[i].adj[0]=0;
-	      edges[i].adj[1]=0;
-	      edges[i].ref = r;
-	      edges[i].on=0;
-       }
-      
-}
-
-//////////////////////////////////////////////////
-
-void  Triangles::Read_amdba(MeshIstream &f_in )
-{
-	long int verbosity=0;
-
-  Metric M1(1);
-  if (verbosity>1)
-    cout << "  -- ReadMesh .amdba file " <<  f_in.CurrentFile  << endl;
-   	
-     Int4 i;
-     f_in.cm() >> nbv >> nbt ;
-     //    if (verbosity>3)
-       cout << "    nbv = " << nbv  << " nbt = " << nbt << endl;
-     f_in.eol() ;// 
-     nbvx = nbv;
-     nbtx =  2*nbv-2; // for filling The Holes and quadrilaterals 
-     triangles =new Triangle[nbtx];
-     assert(triangles);
-     vertices=new Vertex[nbvx];
-     ordre=new (Vertex* [nbvx]);
-     Int4 j;
-      for ( i=0;i<nbv;i++)
-	{
-	  f_in >> j ;
-	  assert( j >0 && j <= nbv);
-	  j--;
-	  f_in >> vertices[j].r.x >>   vertices[j].r.y >> vertices[j].ReferenceNumber;
-	   vertices[j].m=M1;
-	   vertices[j].DirOfSearch=NoDirOfSearch;
-	}
-     
-      for (     i=0;i<nbt;i++)
-       {
-	 Int4 i1,i2,i3,ref;
-	   f_in >> j ;
-	   assert( j >0 && j <= nbt);
-	   j--;
-	   f_in >> i1 >>  i2 >>   i3 >> ref;
-	 triangles[j]  = Triangle(this,i1-1,i2-1,i3-1);
-	 triangles[j].color =ref;
-       }
-      f_in.eol() ;// 
-      
-  // cerr << " a faire " << endl;
-  //MeshError(888);
-}
-
-
-Triangles::Triangles(const char * filename,Real8 cutoffradian) 
-: Gh(*(new Geometry())), BTh(*this)
-{ 
-
-  //  Int4 beginquad=0,begintria=0;
-  // Int4 endquad=0;endtria=0;
-  //int type_file=0;
-
-  int lll = strlen(filename);
-  int  am_fmt = !strcmp(filename + lll - 7,".am_fmt");
-  int  amdba = !strcmp(filename + lll - 6,".amdba");
-  int  am = !strcmp(filename + lll - 3,".am");
-  int  nopo = !strcmp(filename + lll - 5,".nopo");
-  int  msh = !strcmp(filename + lll - 4,".msh");
-  int  ftq = !strcmp(filename + lll - 4,".ftq");
-
- // cout << " Lecture type  :" << filename + lll - 7 <<":" <<am_fmt<<  endl;
-
-  char * cname = new char[lll+1];
-  strcpy(cname,filename);
-  Int4 inbvx =0;
-  PreInit(inbvx,cname);
-  OnDisk = 1;
-  //  allocGeometry = &Gh; // after Preinit ; 
-
-  MeshIstream f_in (filename);
-  
-  if (f_in.IsString("MeshVersionFormatted"))
-    {
-      int version ;
-      f_in >> version ;
-      Read(f_in,version,cutoffradian);
-    }
-  else {     
-    if (am_fmt) Read_am_fmt(f_in);
-    else if (am) Read_am(f_in);
-    else if (amdba) Read_amdba(f_in);
-    else if (msh) Read_msh(f_in);
-    else if (nopo) Read_nopo(f_in);
-    else if (ftq) Read_ftq(f_in);
-    else 
-      { 
-	cerr << " Unkown type mesh " << filename << endl;
-	MeshError(2);
-      }
-      ConsGeometry(cutoffradian);
-      Gh.AfterRead();    
-  }
-
-  SetIntCoor();
-  FillHoleInMesh();
-  // Make the quad ---
-   
-}
-
 Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
 
@@ -1148,14 +276,4 @@
   SetIntCoor();
   FillHoleInMesh();
-}
-
-void Geometry::ReadGeometry(const char * filename) {
-	long int verbosity=0;
-
-  OnDisk = 1;
-  if(verbosity>1)
-      cout << "  -- ReadGeometry " << filename << endl;
-  MeshIstream f_in (filename);
-  ReadGeometry(f_in,filename);
 }
 
Index: /issm/trunk/src/c/Bamgx/MeshWrite.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/MeshWrite.cpp	(revision 2794)
+++ /issm/trunk/src/c/Bamgx/MeshWrite.cpp	(revision 2795)
@@ -14,519 +14,4 @@
 
 namespace bamg {
-
-void Triangles::Write(const char * filename,const TypeFileMesh typein )
-{
-	long int verbosity=0;
-
-  TypeFileMesh type = typein;
-  const char * gsuffix=".gmsh";
-  int ls=0;
-  int lll = strlen(filename);
-  if (type==AutoMesh)
-    {
-      type = BDMesh;
-      if      (!strcmp(filename + lll - (ls=7),".am_fmt")) type = am_fmtMesh;
-      else if (!strcmp(filename + lll - (ls=6),".amdba"))  type = amdbaMesh;
-      else if (!strcmp(filename + lll - (ls=3),".am"))     type = amMesh;
-      else if (!strcmp(filename + lll - (ls=5),".nopo"))   type = NOPOMesh;
-      else if (!strcmp(filename + lll - (ls=4),".msh"))    type = mshMesh;
-      else if (!strcmp(filename + lll - (ls=4),".ftq"))    type = ftqMesh;
-      else if (!strcmp(filename + lll - (ls=7),".AM_FMT")) type = am_fmtMesh;
-      else if (!strcmp(filename + lll - (ls=6),".AMDBA"))  type = amdbaMesh;
-      else if (!strcmp(filename + lll - (ls=3),".AM"))     type = amMesh;
-      else if (!strcmp(filename + lll - (ls=5),".NOPO"))   type = NOPOMesh;
-      else if (!strcmp(filename + lll - (ls=4),".MSH"))    type = mshMesh;
-      else if (!strcmp(filename + lll - (ls=4),".FTQ"))    type = ftqMesh;
-      else ls=0;
-    } 
-  if (verbosity>1)
-    {
-      cout << "  -- Writing the file " << filename << " of type " ;
-     switch (type) 
-       {
-       case BDMesh     :  cout << " BD Mesh "  ; break;
-       case NOPOMesh   :  cout << " NOPO "     ; break;
-       case amMesh     :  cout << " am "       ; break;
-       case am_fmtMesh :  cout << " am_fmt "   ; break;
-       case amdbaMesh  :  cout << " amdba "    ; break;
-       case ftqMesh    :  cout << " ftq "      ; break;
-       case mshMesh    :  cout << " msh "      ; break;
-	default: 
-	  cerr << endl 
-	       <<  " Unkown type mesh file " << (int) type << " for Writing " << filename <<endl;
-	  MeshError(1);
-       }     
-     Int4 NbOfTria =  nbt-2*NbOfQuad-NbOutT ;
-     if (NbOfTria)      cout << " NbOfTria = " << NbOfTria;
-     if (NbOfQuad)      cout << " NbOfQuad = " << NbOfQuad;
-     if (nbe)   	cout << " NbOfRefEdge = " << nbe ;
-     cout    << endl;
-
-    }
-  ofstream f(filename /*,ios::trunc*/);
- f.precision(12);
-
-   if (f)
-     switch (type) 
-       {
-       case BDMesh     : 
-          {
-             if ( ! Gh.OnDisk)
-              {
-                 delete [] Gh.name;
-                 Gh.name = new char[lll+1+strlen(gsuffix)];
-                 strcpy(Gh.name,filename);
-                 if (Gh.name[lll-ls-1]=='.') strcpy(Gh.name+lll-ls, gsuffix+1);
-                 else strcpy(Gh.name+lll-ls,gsuffix);
-                cout << " write geo in " << Gh.name << endl;
-                  ofstream f(Gh.name) ;
-                  f << Gh ;
-                  Gh.OnDisk=true;
-              }
-	     f << *this     ;
-	     break;
-           }
-       case NOPOMesh   :  Write_nopo(f)  ; break;
-       case amMesh     :  Write_am(f)    ; break;
-       case am_fmtMesh :  Write_am_fmt(f); break;
-       case amdbaMesh  :  Write_amdba(f) ; break;
-       case ftqMesh    :  Write_ftq(f)   ; break;
-       case mshMesh    :  Write_msh(f)   ; break;
-	default: 
-	  cerr << " Unkown type mesh file " << (int) type << " for Writing " << filename <<endl;
-	  MeshError(1);
-       }
-   else
-     {
-       cerr << " Error openning file " << filename << endl;
-       MeshError(1);
-     }
-   if(verbosity>5)
-   cout << "end write" << endl;
-      
-}
-void Triangles::Write_nop5(OFortranUnFormattedFile * f,
-			   Int4 &lnop5,Int4 &nef,Int4 &lgpdn,Int4 ndsr) const 
-{
-	long int verbosity=0;
-
-  ndsr =0;
-  Int4 * reft = new Int4[nbt];
-  //Int4 nbInT = ;
-  ConsRefTriangle(reft);
-  Int4 no5l[20];
-
-  Int4 i5 = 0;
-  Int4 i,j,k=0,l5;
-  //  Int4 ining=0;
-  Int4 imax,imin;
-
-  lgpdn = 0;
-  nef=0;
-  // construction of a liste linked  of edge
-  Edge ** head  = new Edge *[nbv];
-  Edge  ** link = new Edge * [nbe];
-  for (i=0;i<nbv;i++)
-    head[i]=0; // empty liste
-  
-  for (i=0;i<nbe;i++)
-    { 
-      j = Min(Number(edges[i][0]),Number(edges[i][1]));
-      link[i]=head[j];
-      head[j]=edges +i;
-    }
-  for ( i=0;i<nbt;i++)
-    {
-      no5l[0]    = 0;
-      Int4 kining=0;
-      Int4 ining=0;
-      Int4 nmae =0;
-      Int4 np=0;
-      l5 = 2;
-      Triangle & t = triangles[i];
-      Triangle * ta; // 
-      Vertex *v0,*v1,*v2,*v3;
-      if (reft[i]<0) continue;
-      ta = t.Quadrangle(v0,v1,v2,v3);
-      if (!ta)
-	{ // a triangles
-	  no5l[l5++] = Max(subdomains[reft[i]].ref,(Int4) 1);
-	  np = 3;
-	  no5l[l5++] = np;
-	  no5l[0]    = np;
-	  no5l[l5++] = Number(triangles[i][0]) +1;
-	  no5l[l5++] = Number(triangles[i][1]) +1;
-	  no5l[l5++] = Number(triangles[i][2]) +1;
-	  imax = Max3(no5l[l5-1],no5l[l5-2],no5l[l5-3]);
-	  imin = Min3(no5l[l5-1],no5l[l5-2],no5l[l5-3]);
-	  lgpdn = Max(lgpdn,imax-imin);
-          kining=l5++;
-	  // ref of 3 edges 
-	  for (j=0;j<3;j++)
-	    {
-	      no5l[l5] = 0;
-	      int i0 = (int) j;
-	      int i1 = (i0+1) %3;
-	      Int4 j1= no5l[4+i0];
-	      Int4 j2= no5l[4+i1];
-	      Int4 ji = Min(j1,j2)-1;
-	      Int4 ja = j1+j2-2;
-	      Edge * e=head[ji];
-	      while (e)
-		if(Number((*e)[0])+Number((*e)[1]) == ja) 
-		   {
-		     no5l[l5] = e->ref;
-		     break;
-		   }
-		else
-		  e = link[Number(e)];
-	      l5++;
-	    }
-	  if ( no5l[l5-1] || no5l[l5-2] || no5l[l5-3]  )
-	    ining=2;
-	  else 
-	    l5 -= 3;
-	  
-	  no5l[l5++] = triangles[i][0].ref();
-	  no5l[l5++] = triangles[i][1].ref();
-	  no5l[l5++] = triangles[i][2].ref();
-	  if (ining ||  no5l[l5-1] || no5l[l5-2] || no5l[l5-3]  )
-            ining= ining ? ining : 3;
-	  else
-	    l5 -= 3,ining=0;
-
-	  
-	}
-      else if ( &t<ta)
-	{ 
-	  k++;
-	  no5l[l5++] = Max(subdomains[reft[i]].ref,(Int4) 1);
-	  np =4;
-	  no5l[l5++] = np;
-	  no5l[0]    = np;
-	  
-	  no5l[l5++] = Number(v0) +1;
-	  no5l[l5++] = Number(v1) +1;
-	  no5l[l5++] = Number(v2) +1;
-	  no5l[l5++] = Number(v3) +1;
-	  
-	  imax = Max(Max(no5l[l5-1],no5l[l5-2]),Max(no5l[l5-3],no5l[l5-4]));
-	  imin = Min(Min(no5l[l5-1],no5l[l5-2]),Min(no5l[l5-3],no5l[l5-4]));
-	  lgpdn = Max(lgpdn,imax-imin);
-
-
-          kining=l5++;
-	  // ref of the 4 edges 
-	  // ref of 3 edges 
-	  for (j=0;j<4;j++)
-	    {
-	      no5l[l5] = 0;
-	      int i0 = (int) j;
-	      int i1 = (i0+1) %4;
-	      Int4 j1= no5l[4+i0];
-	      Int4 j2= no5l[4+i1];
-	      Int4 ji = Min(j1,j2)-1;
-	      Int4 ja = j1+j2-2;
-	      Edge *e=head[ji];
-	      while (e)
-		if(Number((*e)[0])+Number((*e)[1]) == ja) 
-		   {
-		     no5l[l5] = e->ref;
-		     break;
-		   }
-		else
-		  e = link[Number(e)];
-	      l5++;
-	    }
-	  if ( no5l[l5-1] || no5l[l5-2] || no5l[l5-3] || no5l[l5-4] )
-	    ining=2;
-	  else 
-	    l5 -= 4;
-
-	  no5l[l5++] = v0->ref();
-	  no5l[l5++] = v1->ref();
-	  no5l[l5++] = v2->ref();
-	  no5l[l5++] = v3->ref();
-	  if (ining || no5l[l5-1] || no5l[l5-2] || no5l[l5-3] || no5l[l5-4] )
-            ining= ining ? ining : 3;
-          else
-	    l5 -= 4;
-
-	}
-      else l5=0;
-
-     if (l5)
-       {
-	 if (ining)
-	   {
-	     nef ++;
-	     nmae = l5-kining;
-	     no5l[kining] = ining;
-	   }
-	 else l5--;
-	 no5l[1]=nmae;
-	 // for all ref  
-	 for (j=kining+1;j<l5;j++)
-	   {
-	     no5l[j] = Abs(no5l[j]);
-	     ndsr = Max(ndsr,no5l[j]);
-	   }
-	 
-	 if (f && i < 10 && verbosity > 10)
-	   { 
-	     cout << " e[ " << i << " @" << i5 << "]=";
-	     for (j=0;j<l5;j++)
-	       cout << " " << no5l[j]; 
-	     cout << endl;
-	   }
-	 
-	 if (f)
-	   for (j=0;j<l5;j++)
-	     *f << no5l[j]; 
-	 i5 += l5;
-       }
-    }
-  if(verbosity>10)
-  cout << "   fin write nopo 5 i5=" << i5 << " " << i5*4 << endl;
-  lnop5=i5; 
-  lgpdn++; // add 1
-  delete [] reft;
-  delete [] head;
-  delete [] link;
-  
-}
-
-void Triangles::Write_nopo(ostream &ff) const
-
-{
- Int4  nef=0;
- Int4 lpgdn=0;
- Int4 ndsr=0;
- Int4 i;
- Int4 ndsd=1;
- Int4 lnop5=0;
- 
- OFortranUnFormattedFile f(ff);
- 
- for (i=0;i<NbSubDomains ;i++)
-   ndsd=Max(ndsd,subdomains[i].ref);
- 
- // to compute the lnop5,nef,lpgdn,ndsr parameter 
- Write_nop5(0,lnop5,nef,lpgdn,ndsr);
- 
- f.Record();
- 
- f <<  Int4(13)<<Int4(6)<<Int4(32)<<Int4(0)<<Int4(27)<<Int4(0) ;
- f << Int4(nbv+nbv) ;
- f << lnop5;
- f << Int4(1 )<<Int4(1)<<Int4(1 )<<Int4(1)<<Int4(2)<<Int4(1);
-
- f.Record(33*sizeof(Int4)); 
-
- f << Int4(32) ;
- 
- //char *c=identity;
- time_t timer =time(0);
- char buf[10];
- strftime(buf ,10,"%y/%m/%d",localtime(&timer));
- f.write4(identity,20);
- f.write4(buf,2);
- f.write4("created with BAMG",6);
- f.write4("NOPO",1);
-
- 
- f << Int4(0) << Int4(1) << Int4(0) ;
- f.Record();
- Int4 nbquad= NbOfQuad;
- Int4 nbtria= nbt-NbOutT - 2*NbOfQuad;
-
- cout << " lnop5      = " << lnop5 << endl;
- cout << " nbquad     = " << nbquad << endl;
- cout << " nbtrai     = " << nbtria << endl;
- cout << " lpgdn      = " << lpgdn << endl;
- cout << " nef        = " << nef  << endl;
- cout << " np         = " << nbv  << endl;
- cout << " ndsr       = " << ndsr << endl;
- f << Int4(27)  
-   << Int4(2)  << ndsr     << ndsd    << Int4(1) << nbtria+nbquad
-   << Int4(0)  << Int4(0)  << nbtria  << nbquad  << Int4(0)
-   << Int4(0)  << Int4(0)  << Int4(0) << nef     << Int4(nbv)
-   << Int4(0)  << Int4(0)  << Int4(0) << Int4(0) << Int4(0)
-   << Int4(0)  << nbv      << Int4(2) << lpgdn   << Int4(0)
-   << lnop5    << Int4(1) ;
-  f.Record();
-  f << (Int4) 2*nbv;
-  for (i=0;i<nbv;i++)
-    f << (float) vertices[i].r.x <<  (float) vertices[i].r.y;
-  f.Record();
-  f << lnop5;
-  Write_nop5(&f,lnop5,nef,lpgdn,ndsr);
-  // cout << "fin write nopo" << endl;
-}
-
-void Triangles::Write_am_fmt(ostream &f) const 
-{
-  Int4 i,j;
-  assert(this && nbt);
-  Int4 * reft = new Int4[nbt];
-  Int4 nbInT =    ConsRefTriangle(reft);
-  f.precision(12);
-  f << nbv << " " << nbInT << endl;
-  for (i=0;i<nbt;i++)
-    if(reft[i]>=0)
-      {
-	f << Number(triangles[i][0]) +1 << " " ;
-	f << Number(triangles[i][1]) +1 << " " ;
-	f << Number(triangles[i][2]) +1 << " " ;
-	f << endl;
-      }
-  for (i=0;i<nbv;i++)
-      f << vertices[i].r.x << " " << vertices[i].r.y << endl;
-   for (j=i=0;i<nbt;i++) 
-     if (reft[i]>=0)
-       f << subdomains[reft[i]].ref  << (j++%10 == 9 ?  '\n' : ' ');
-   f << endl;
-   for (i=0;i<nbv;i++)
-     f << vertices[i].ref()  << (i%10 == 9 ?  '\n' : ' ');
-   f << endl;
-   delete [] reft;
-
-
-}
-
-void Triangles::Write_am(ostream &ff) const 
-{
-  OFortranUnFormattedFile f(ff);  
-  Int4 i,j;
-  assert(this && nbt);
-  Int4 * reft = new Int4[nbt];
-  Int4 nbInT =    ConsRefTriangle(reft);
-  f.Record();
-  f << nbv << nbInT ;
-  f.Record();
-  for (i=0;i<nbt;i++)
-    if(reft[i]>=0)
-      {
-	f << Number(triangles[i][0]) +1 ;
-	f << Number(triangles[i][1]) +1 ;
-	f << Number(triangles[i][2]) +1 ;
-      }
-  for (i=0;i<nbv;i++)
-    {
-      float x= vertices[i].r.x;
-      float y= vertices[i].r.y;
-      f << x << y ;
-    }
-  for (j=i=0;i<nbt;i++) 
-    if (reft[i]>=0)
-      f << subdomains[reft[i]].ref;
-  for (i=0;i<nbv;i++)
-    f << vertices[i].ref() ;
-  delete [] reft;
-}
-
-void Triangles::Write_ftq(ostream &f) const 
-{
-
-  Int4 i;
-  assert(this && nbt);
-  Int4 * reft = new Int4[nbt];
-  Int4 nbInT =    ConsRefTriangle(reft);
-  f.precision(12);
-  Int4 nele = nbInT-NbOfQuad;
-  Int4 ntri =  nbInT-2*NbOfQuad;
-  Int4 nqua =  NbOfQuad;
-
-  f << nbv << " " << nele << " " << ntri <<  " " << nqua << endl;
-  Int4 k=0;
-  for( i=0;i<nbt;i++)
-    { 
-      Triangle & t = triangles[i];
-      Triangle * ta; // 
-      Vertex *v0,*v1,*v2,*v3;
-      if (reft[i]<0) continue;
-      ta = t.Quadrangle(v0,v1,v2,v3);
-      if (!ta)
-	{ // a triangles
-	  f << "3 " 
-	    << Number(triangles[i][0]) +1 << " " 
-	    << Number(triangles[i][1]) +1 << " " 
-	    << Number(triangles[i][2]) +1 << " " 
-	    << subdomains[reft[i]].ref << endl;
-	  k++;
-	}
-      if ( &t<ta)
-	{ 
-	  k++;
-	  f << "4 " << Number(v0)+1 << " " << Number(v1)+1  << " "  
-	    << Number(v2)+1 << " "  << Number(v3)+1 << " "  
-	    << subdomains[reft[i]].ref << endl;
-	}
-    }
-  assert(k == nele);
-  
-  for (i=0;i<nbv;i++)
-    f << vertices[i].r.x << " " << vertices[i].r.y 
-      << " " <<  vertices[i].ref() << endl;
-  delete [] reft;
-  
-  
-}
-void Triangles::Write_msh(ostream &f) const 
-{
-  Int4 i;
-  assert(this && nbt);
-  Int4 * reft = new Int4[nbt];
-  Int4 nbInT =    ConsRefTriangle(reft);
-  f.precision(12);
-  f << nbv << " " << nbInT << " " << nbe <<  endl;
-
-  for (i=0;i<nbv;i++)
-    f << vertices[i].r.x << " " << vertices[i].r.y << " " 
-      << vertices[i].ref() <<   endl;
-
-  for (i=0;i<nbt;i++)
-    if(reft[i]>=0)
-      f << Number(triangles[i][0]) +1 << " " 
-	<< Number(triangles[i][1]) +1 << " " 
-	<< Number(triangles[i][2]) +1 << " " 
-	<< subdomains[reft[i]].ref << endl;
-  
-
-  for (i=0;i<nbe;i++)
-    f << Number(edges[i][0]) +1 << " "  << Number(edges[i][1]) +1 
-      << " " << edges[i].ref << endl;
-      
-   delete [] reft;
-
-}
-
-void Triangles::Write_amdba(ostream &f) const 
-{
-  assert(this && nbt);
-
-  Int4 i,j;
-  Int4 * reft = new Int4[nbt];
-  Int4 nbInT =    ConsRefTriangle(reft);
-  f << nbv << " " << nbInT << endl;
-  cout.precision(12);
-  for (i=0;i<nbv;i++)
-    f << i+1 << " " 
-      << vertices[i].r.x 
-      << " " << vertices[i].r.y 
-      << " " << vertices[i].ref() << endl;
-  j=1;
-  for (i=0;i<nbt;i++)
-    if(reft[i]>=0)
-	f << j++ << " " 
-	  << Number(triangles[i][0]) +1 << " " 
-	  << Number(triangles[i][1]) +1 << " " 
-	  << Number(triangles[i][2]) +1 << " " 
-	  << subdomains[reft[i]].ref  << endl ;
-	f << endl;
-   delete [] reft;
-
-
-}
 
 void Triangles::WriteBamgMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
@@ -729,254 +214,4 @@
 	}
 }
-
-void Triangles::Write(const char * filename)
-{
-  ofstream f(filename);
-  if (f)
-    {
-       if (name) delete name;
-       name = new char[strlen(filename)+1];
-       strcpy(name,filename);
-       OnDisk =1;
-       f << *this;
-    }
-}
-void Triangles::WriteElements(ostream& f,Int4 * reft ,Int4 nbInT) const
-   { 
-     const Triangles & Th= *this;
-	 long int verbosity=0;
-
-     // do triangle and quad 
-     if(verbosity>9) 
-       cout  << " In Triangles::WriteElements " << endl
-	     << "   Nb of In triangles " << nbInT-Th.NbOfQuad*2 << endl
-	     << "   Nb of Quadrilaterals " <<  Th.NbOfQuad << endl
-	     << "   Nb of in+out+quad  triangles " << Th.nbt << " " << nbInT << endl;
-	 
-     Int4 k=nbInT-Th.NbOfQuad*2;
-     Int4 num =0;
-     if (k>0) {
-       f << "\nTriangles\n"<< k << endl;
-       for(Int4 i=0;i<Th.nbt;i++)
-	 { 
-	   Triangle & t = Th.triangles[i];
-	   if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) ))
-	    { k--;
-	      f << Th.Number(t[0])+1 << " " << Th.Number(t[1])+1 
-		   << " "  << Th.Number(t[2])+1  << " " << Th.subdomains[reft[i]].ref << endl;
-		   reft[i] = ++num;
-		 }
-	 }
-     } 
-     if (Th.NbOfQuad>0) {
-       f << "\nQuadrilaterals\n"<<Th.NbOfQuad << endl;
-       k = Th.NbOfQuad;
-       for(Int4 i=0;i<Th.nbt;i++)
-	 { 
-	   Triangle & t = Th.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)
-	      { 
-		k--;
-		f << Th.Number(v0)+1 << " " << Th.Number(v1)+1  << " "  
-		  << Th.Number(v2)+1 << " "  << Th.Number(v3)+1 << " "  
-		  << Th.subdomains[reft[i]].ref << endl;
-		  reft[i] = ++num;
-		  reft[Number(ta)] = num;
-	      }
-	 }
-       assert(k==0);
-     }
-     // warning reft is now the element number 
-   }
-
-ostream& operator <<(ostream& f, const   Triangles & Th) 
- {
-  //  Th.FindSubDomain();
-   // warning just on say the class is on the disk
-  //  ((Triangles *) &Th)->OnDisk = 1;
-
-   Int4 * reft = new Int4[Th.nbt];
-   Int4 nbInT =    Th.ConsRefTriangle(reft);
-   {
-     f << "MeshVersionFormatted 0" <<endl;
-     f << "\nDimension\n"  << 2 << endl;
-     f << "\nIdentifier\n" ;
-     WriteStr(f,Th.identity);
-     f << "\n\nGeometry\n" ;
-     if( Th.Gh.OnDisk)
-       WriteStr(f,Th.Gh.name),     f <<endl;
-     else
-       { // empty file name -> geom in same file
-	 f << "\"\"" << endl << endl;
-	 f << "# BEGIN of the include geometry file because geometry is not on the disk"
-	   << Th.Gh << endl;
-	 f << "End" << endl 
-	   << "# END of the include geometrie file because geometry is not on the disk"
-	   << endl ;
-       }
-   }
-   { 
-     f.precision(12);
-     f << "\nVertices\n" << Th.nbv <<endl;
-     for (Int4 i=0;i<Th.nbv;i++)
-       {
-	 Vertex & v =  Th.vertices[i];
-	 f << v.r.x << " " << v.r.y << " " << v.ref() << endl;
-       }
-   }
-  Int4 ie; 
-   {
-     f << "\nEdges\n"<< Th.nbe << endl;
-     for(ie=0;ie<Th.nbe;ie++)
-       { 
-	 Edge & e = Th.edges[ie];
-	 f << Th.Number(e[0])+1 << " " << Th.Number(e[1])+1;
-	f << " " << e.ref <<endl;
-       }
-     if(Th.NbCrackedEdges)
-       {
-	 f << "\nCrackedEdges\n"<< Th.NbCrackedEdges << endl;
-	 for( ie=0;ie<Th.NbCrackedEdges;ie++)
-	   { 
-	     Edge & e1 = *Th.CrackedEdges[ie].a.edge;
-	     Edge & e2 = *Th.CrackedEdges[ie].b.edge;
-	     f << Th.Number(e1)+1 << " " << Th.Number(e2)+1 <<endl;;
-	   }
-       }
-   }
-
-   Th.WriteElements(f,reft,nbInT);
-   {
-     f << "\nSubDomainFromMesh\n" << Th.NbSubDomains<< endl ;
-     for (Int4 i=0;i<Th.NbSubDomains;i++)
-       f << 3 << " " << reft[Th.Number(Th.subdomains[i].head)] << " " << 1 << " " 
-	 <<  Th.subdomains[i].ref << endl;
-     
-   }
-   if (Th.Gh.NbSubDomains)
-     {
-        f << "\nSubDomainFromGeom\n" << Th.Gh.NbSubDomains << endl ;
-       for (Int4 i=0;i<Th.NbSubDomains;i++)
-	 {  
-	 f << 2 << " " << Th.Number(Th.subdomains[i].edge)+1 << " " 
-	   <<  Th.subdomains[i].sens  << " " <<  Th.Gh.subdomains[i].ref << endl;
-	 } 
-   }
-   {
-     f << "\nVertexOnGeometricVertex\n"<<  Th.NbVerticesOnGeomVertex << endl;
-     for (Int4 i0=0;i0<Th.NbVerticesOnGeomVertex;i0++)
-       {
-	 VertexOnGeom & v =Th.VerticesOnGeomVertex[i0];
-	 assert(v.OnGeomVertex()) ;
-	 f << " " << Th.Number(( Vertex *)v)+1  
-	   << " " << Th.Gh.Number(( GeometricalVertex * )v)+1 
-	   << endl;
-       }
-   }
-   { 
-     f << "\nVertexOnGeometricEdge\n"<<  Th.NbVerticesOnGeomEdge << endl;
-     for (Int4 i0=0;i0<Th.NbVerticesOnGeomEdge;i0++)
-       {
-	 const VertexOnGeom & v =Th.VerticesOnGeomEdge[i0];
-	 assert(v.OnGeomEdge()) ;   
-	 f << " " << Th.Number((Vertex * )v)+1  ;
-	 f << " " << Th.Gh.Number((const  GeometricalEdge * )v)+1  ;
-	 f << " " << (Real8 ) v << endl;
-       }
-   }
-   {
-     Int4 i0,k=0;
-
-     for (i0=0;i0<Th.nbe;i0++)
-       if ( Th.edges[i0].on ) k++;
-     
-     f << "\nEdgeOnGeometricEdge\n"<< k << endl;
-      for (i0=0;i0<Th.nbe;i0++)
-		 if ( Th.edges[i0].on ) 
-		  f << (i0+1) << " "  << (1+Th.Gh.Number(Th.edges[i0].on)) <<  endl;
-		if (Th.NbCrackedEdges)
-	{
-	  f << "\nCrackedEdges\n"<< Th.NbCrackedEdges << endl;	  
-	  for(i0=0;i0< Th.NbCrackedEdges; i0++) 
-	    {
-	      f << Th.Number(Th.CrackedEdges[i0].a.edge) << " " ;
-	      f  << Th.Number(Th.CrackedEdges[i0].b.edge) << endl;
-	    }
-	}
-   }  
-   if (&Th.BTh != &Th && Th.BTh.OnDisk && Th.BTh.name) 
-     {
-       int *mark=new int[Th.nbv];
-       Int4 i;
-       for (i=0;i<Th.nbv;i++)
-	 mark[i]=-1;
-       f << "\nMeshSupportOfVertices\n" <<endl;
-       WriteStr(f,Th.BTh.name);
-       f <<endl;
-       f << "\nIdentityOfMeshSupport" << endl;
-       WriteStr(f,Th.BTh.identity);
-       f<<endl;
-
-       f << "\nVertexOnSupportVertex" << endl;
-       f<< Th.NbVertexOnBThVertex << endl;
-       for(i=0;i<Th.NbVertexOnBThVertex;i++) {
-	 const VertexOnVertex & vov = Th.VertexOnBThVertex[i];
-	 Int4 iv = Th.Number(vov.v);
-	 mark[iv] =0;
-	 f << iv+1<< " " << Th.BTh.Number(vov.bv)+1 << endl;}
-
-       f << "\nVertexOnSupportEdge" << endl;
-       f << Th.NbVertexOnBThEdge << endl;
-       for(i=0;i<Th.NbVertexOnBThEdge;i++) {
-	 const VertexOnEdge & voe = Th.VertexOnBThEdge[i];
-	 Int4 iv = Th.Number(voe.v);
-	 //	 assert(mark[iv] == -1]);
-	 mark[iv] = 1;
-	 f << iv+1 << " " << Th.BTh.Number(voe.be)+1 << " " << voe.abcisse <<  endl;}
-       
-       f << "\nVertexOnSupportTriangle" << endl;   
-       Int4 k = Th.nbv -  Th.NbVertexOnBThEdge - Th.NbVertexOnBThVertex;
-       f << k << endl;
-       //       Int4 kkk=0;
-       CurrentTh=&Th.BTh;
-       for (i=0;i<Th.nbv;i++) 
-	 if (mark[i] == -1) {
-	   k--;
-	   Icoor2 dete[3];
-	   I2 I = Th.BTh.toI2(Th.vertices[i].r);
-	   Triangle * tb = Th.BTh.FindTriangleContening(I,dete);
-	   if (tb->link) // a true triangle
-	     {
-	       Real8 aa= (Real8) dete[1]/ tb->det, bb= (Real8) dete[2] / tb->det;
-	       f << i+1 << " " << Th.BTh.Number(tb)+1 << " " << aa << " " << bb << endl ;
-	     }
-	   else 
-	     {
-	       double aa,bb,det[3];
-	       TriangleAdjacent ta=CloseBoundaryEdgeV2(I,tb,aa,bb);
-	       int k = ta;
-	       det[VerticesOfTriangularEdge[k][1]] =aa;
-	       det[VerticesOfTriangularEdge[k][0]] = bb;
-	       det[OppositeVertex[k]] = 1- aa -bb;
-	       Triangle * tb = ta;
-	       f << i+1 << Th.BTh.Number(tb)+1 << " " << det[1] << " " << det[2] <<endl;
-	     }
-	 }
-       assert(!k);
-       delete [] mark;
-	 
-
-     }
-   f << "\nEnd" << endl;
-   //  Th.ConsLinkTriangle();
-   delete [] reft;
-   return f;
-   
-}
-
-
 
 void Geometry::Write(const char * filename)