// -*- Mode : c++ -*-//// SUMMARY  :      // USAGE    :        // ORG      : // AUTHOR   : Frederic Hecht// E-MAIL   : hecht@ann.jussieu.fr///*  This file is part of 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 <stdio.h>#include <string.h>#include <math.h>#include <time.h>#include "Meshio.h"#include "Mesh2.h"#include "QuadTree.h"#include "SetOfE4.h"#ifdef __MWERKS__#pragma optimization_level 2//#pragma inline_depth 1#endif#ifdef DRAWING1  extern bool withrgraphique ;#endifnamespace bamg {static const  Direction NoDirOfSearch=Direction();void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian){  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){  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() ;//