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){ #ifdef DRAWING1   if (!withrgraphique) {initgraphique();withrgraphique=true;}   #endif  //  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 ---   }void Geometry::ReadGeometry(const char * filename){  OnDisk = 1;  if(verbosity>1)      cout << " -- ReadGeometry " << filename << endl;  MeshIstream f_in (filename);  ReadGeometry(f_in,filename);}void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)  {  if(verbosity>1)    cout << " -- ReadGeometry " << filename << endl;  assert(empty());  nbiv=nbv=nbvx=0;  nbe=nbt=nbtx=0;  NbOfCurves=0; // BeginOfCurve=0;  name=new char [strlen(filename)+1];  strcpy(name,filename);  Real8 Hmin = HUGE_VAL;// the infinie value //  Real8 MaximalAngleOfCorner = 10*Pi/180; ;   Int4 hvertices =0;  Int4 i;  Int4 Version,dim=0;  int field=0;  int showfield=0;  int NbErr=0;  while (f_in.cm())     {       field=0;      // warning ruse for on allocate fiedname at each time       char fieldname[256];      f_in.cm() >> fieldname ;      f_in.err();      if(f_in.eof()) break;//      cout <<  fieldname <<  " line " << LineNumber  <<endl;      if (!strcmp(fieldname,"MeshVersionFormatted") )        f_in  >> Version ;      else if (!strcmp(fieldname,"End"))	break;      else if (!strcmp(fieldname,"end"))	break;      else if (!strcmp(fieldname,"Dimension"))        {          f_in   >>  dim ;	  if(verbosity>5) 	    cout << "     Geom Record Dimension dim = " << dim << endl;                  assert(dim ==2);         }       else if (!strcmp(fieldname,"hVertices"))         { 	   if (nbv <=0) {	     cerr<<"Error: the field Vertex is not found before hVertices " << filename<<endl;	     NbErr++;}       	   if(verbosity>5) 	    cout << "     Geom Record hVertices nbv=" << nbv <<  endl;	   hvertices =1;           for (i=0;i< nbv;i++) 	     {	       Real4 h;	       f_in  >>  h ; 	       vertices[i].m = Metric(h);	     }	 }       else if (!strcmp(fieldname,"MetricVertices"))         { hvertices =1;	   if (nbv <=0) {	     cerr<<"Error: the field Vertex is not found before MetricVertices " << filename<<endl;	     NbErr++;}                  if(verbosity>5) 	     cout << "     Geom Record MetricVertices nbv =" << nbv <<  endl;           for (i=0;i< nbv;i++) 	     {	       Real4 a11,a21,a22;	       f_in  >>  a11 >> a21 >> a22  ; 	       vertices[i].m = Metric(a11,a21,a22);	     }	 }       else if (!strcmp(fieldname,"h1h2VpVertices"))         { hvertices =1;	   if (nbv <=0) {	     cerr<<"Error: the field Vertex is not found before h1h2VpVertices " << filename<<endl;	     NbErr++;}                  if(verbosity>5) 	     cout << "     Geom Record h1h2VpVertices nbv=" << nbv << endl;           for (i=0;i< nbv;i++) 	     {	       Real4 h1,h2,v1,v2;	       f_in  >> h1 >> h2 >>v1 >>v2 ; 	       vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));	     }	 }      else if (!strcmp(fieldname,"Vertices"))        {           assert(dim ==2);          f_in   >>  nbv ;	  //          if(LineError) break;          nbvx = nbv;                    vertices = new GeometricalVertex[nbvx];	  if(verbosity>5) 	    cout << "     Geom Record Vertices nbv = " << nbv << "vertices = " << vertices<<endl;          assert(nbvx >= nbv);          nbiv = nbv;          for (i=0;i<nbv;i++) {            f_in  >> vertices[i].r.x  ;            // if(LineError) break;            f_in  >> vertices[i].r.y ;	    // if(LineError) break;            f_in >>   vertices[i].ReferenceNumber   ;            vertices[i].DirOfSearch=NoDirOfSearch;	    //            vertices[i].m.h = 0;            vertices[i].color =0;            vertices[i].Set();}	  // if(LineError) break;	    pmin =  vertices[0].r;	    pmax =  vertices[0].r;	    // recherche des extrema des vertices pmin,pmax	    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));	    assert(coefIcoor >0);	    if (verbosity>5) {	      cout << "     Geom: min="<< pmin << "max ="<< pmax << " hmin = " << MinimalHmin() <<  endl;}        }      else if(!strcmp(fieldname,"MaximalAngleOfCorner")||!strcmp(fieldname,"AngleOfCornerBound"))        {                   f_in >> MaximalAngleOfCorner;              	  if(verbosity>5) 	    cout << "     Geom Record MaximalAngleOfCorner " << MaximalAngleOfCorner <<endl;          MaximalAngleOfCorner *= Pi/180;        }      else if (!strcmp(fieldname,"Edges"))        {	  if (nbv <=0) {	    cerr<<"Error: the field edges is not found before MetricVertices " << filename<<endl;	    NbErr++;}   	  else 	    {	      int i1,i2;	      R2 zero2(0,0);	      f_in   >>  nbe ;	      	      edges = new GeometricalEdge[nbe];	      if(verbosity>5) 		cout << "     Record Edges: Nb of Edge " << nbe <<endl;	      assert(edges);	      assert (nbv >0); 	      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  ;		  		  i1--;i2--; // for C index		  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;		  		  if(verbosity>3) 		    cout << "     Geom Hmin " << Hmin << endl;		}	      	    }	}      else if (!strcmp(fieldname,"EdgesTangence") ||!strcmp(fieldname,"TangentAtEdges")  )        {           int n,i,j,k;          R2 tg;          f_in  >> n ;          	  if(verbosity>5) 	    cout << "     Record TangentAtEdges: Nb of Edge " << n <<endl;                    for (k=0;k<n;k++)            {	      f_in  >>  i  >> j ;	      f_in >> tg.x  >> tg.y ;	      assert( i <= nbe );	      assert( i > 0 );	      assert ( j == 1 || j==2 );	      i--;j--;// for C index	      edges[i].tg[j] = tg;            }        }      else if (!strcmp(fieldname,"Corners"))        {           int i,j,n;          f_in  >> n ;	  if(verbosity>5) 	    cout << "     Record Corner: Nb of Corner " << n <<endl;                    for (i=0;i<n;i++) {                 f_in  >>  j ;            assert( j <= nbv );            assert( j > 0 );            j--;            vertices[j].SetCorner();            vertices[j].SetRequired();  }        }      else if (!strcmp(fieldname,"RequiredVertices"))        {           int i,j,n;          f_in  >> n ;          for (i=0;i<n;i++) {                 f_in  >>  j ;            assert( j <= nbv );            assert( j > 0 );            j--;            vertices[j].SetRequired();  }      }      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,"SubDomain") || !strcmp(fieldname,"SubDomainFromGeom"))      { 	f_in   >>  NbSubDomains ;	if (NbSubDomains>0) 	  {	    subdomains = new GeometricalSubDomain[  NbSubDomains];	    Int4 i0,i1,i2,i3;	    for (i=0;i<NbSubDomains;i++) 	      {		f_in  >> i0  >>i1 		      >> i2  >>i3 ; 				assert(i0 == 2);		assert(i1<=nbe && i1>0);		subdomains[i].edge=edges + (i1-1);		subdomains[i].sens = (int) i2;		subdomains[i].ref = i3;	      }	  }      }      else	{ // unkown field	  field = ++showfield;	  if(showfield==1) // just to show one time 	    if (verbosity>3)	      cout << "    Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;	}      showfield=field; // just to show one time     } // while !eof()  // generation  de la geometrie   // 1 construction des aretes   // construire des aretes en chaque sommets     if (nbv <=0) {    cerr<<"Error: the field Vertex is not found in " << filename<<endl;    NbErr++;}  if(nbe <=0) {    cerr <<"Error: the field Edges is not found in "<< filename<<endl      ;NbErr++;}  if(NbErr) MeshError(1); }}  // end of namespace bamg